MODELING AND SIMULATION OF A FLUTTERING BIOINSPIRED SUBMERSIBLE
By
Aren Michael Hellum

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Mechanical Engineering
2011

ABSTRACT
MODELING AND SIMULATION OF A FLUTTERING BIOINSPIRED
SUBMERSIBLE
By
Aren Michael Hellum
A bioinspired submersible propelled by a ﬂuttering, ﬂuid-conveying tail was designed and
analyzed. A simple model of the submersible was created by treating the rigid hull of
the submersible as a rigid body boundary condition on the ﬂuid-conveying tail. Curves of
neutral stability were determined in the internal/external velocity parameter space for several
values of rigid body mass, and the thrust produced by these neutrally-stable waveforms was
determined using Lighthill’s methods. The power required to produce these waveforms
and their eﬃciency were also determined. The eﬃciency was found to be greater than
50%, comparable to a small marine propeller, but below the 90% quoted in ﬁsh propulsion
literature.
A more general model of the submersible which removed many of the linearizing assumptions made in the simpliﬁed model was also developed and solved numerically. This model is
comprised of the equations of motion for a non-inertial frame ﬁxed on the rigid head of the
submersible and the equations of motion for a ﬂuid-conveying tail within that non-inertial
frame. Simulations were made using several ﬂexible tail geometries, as well as for a rigid
tube and a rigid tail which was dimensionally identical to one of the geometries. The forward speed of the ﬂuttering ﬂexible tail conﬁguration was found to be higher than than that
of the rigid tube conﬁguration for one geometry, and higher than that of the rigid tail for
most geometries. A prototype of this submersible was constructed, and qualitative features
and ﬁndings of the general model were veriﬁed. The general model was also extended to

incorporate a time-variable velocity of the conveyed ﬂuid, and a functional description of
this velocity was found which caused the submersible to turn without the incorporation of
additional actuators on the tail or hull-mounted ﬁns.

Copyright by
AREN MICHAEL HELLUM
2011

TABLE OF CONTENTS

List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
List of Figures

viii

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

1 Motivation

ix
1

2 Dynamics of a Pipe Conveying Fluid with a Non-Uniform Velocity
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Background - Uniform Velocity Proﬁle . . . . . . . . . . . . . . . . .
2.3 Non-Uniform Velocity Proﬁle . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Triple plug ﬂow model . . . . . . . . . . . . . . . . . . . . . .
2.3.2 N-plug ﬂow model . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Analysis of a Cantilever Pipe . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Solution of the diﬀerential equation . . . . . . . . . . . . . . .
2.4.2 Determination of momentum ﬂux correction factor . . . . . .
2.5 Model Comparison:
Uniform and Non-Uniform Flow . . . . . . . . . . . . . . . . . . . .
2.5.1 Turbulent ﬂow . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2 Laminar ﬂow . . . . . . . . . . . . . . . . . . . . . . . . . . .

Proﬁle
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .

5
5
6
8
8
11
13
13
14

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

16
16
20

3 Simpliﬁed Dynamics of a Fluid-Conveying, Fluid-Immersed Pipe
to a Rigid Body
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Fluid-conveying pipes . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Slender body swimming . . . . . . . . . . . . . . . . . . . . .
3.3 Fluid-conveying pipe aﬃxed to a rigid body . . . . . . . . . . . . . .
3.3.1 Simpliﬁed boundary conditions . . . . . . . . . . . . . . . . .
3.3.2 Method of analysis . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Stability, thrust and eﬃciency . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Neutral stability . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Thrust characteristics . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Hydrodynamic eﬃciency . . . . . . . . . . . . . . . . . . . . .

Aﬃxed
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

23
23
24
24
28
32
32
35
38
38
40
42

4 General Dynamics of a Fluid-Conveying, Fluid-Immersed Pipe Aﬃxed to
a Rigid Body
46
4.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
v

4.2
4.3
4.4

4.5

4.6

4.7

Introduction . . . . . . . . . . . . . . . . . . . .
Assumptions . . . . . . . . . . . . . . . . . . . .
Dynamic Model - Newtonian Formulation . . .
4.4.1 Force and Moment Equations . . . . . .
4.4.2 Added Mass . . . . . . . . . . . . . . . .
4.4.3 Expansion of a General Volume Term . .
4.4.4 Individual Volume Terms . . . . . . . . .
4.4.4.1 Rigid Head: Ωh . . . . . . . . .
4.4.4.2 Flexible Tail: Ωt . . . . . . . .
4.4.4.3 Head Internal Fluid: Ω1 . . . .
4.4.4.4 Tail Internal Fluid: Ω2 . . . . .
4.4.4.5 Head External Fluid: Ω3 . . .
4.4.4.6 Tail External Fluid: Ω4 . . . .
4.4.5 Individual Surface Terms . . . . . . . . .
4.4.5.1 Jet Inlet: Γ1 . . . . . . . . . .
4.4.5.2 Jet Outlet: Γ2 . . . . . . . . .
4.4.5.3 External Flow Head Inlet: Γ3 .
4.4.5.4 External Flow Head Outlet: Γ4
4.4.5.5 External Flow Tail Inlet: Γ5 . .
4.4.5.6 External Flow Tail Outlet: Γ6 .
4.4.6 Combined Scalar Equations . . . . . . .
4.4.7 Pressure and Drag Forces . . . . . . . .
Dynamics of the Flexible Tail . . . . . . . . . .
4.5.1 Equations of Motion . . . . . . . . . . .
4.5.2 Discretization . . . . . . . . . . . . . . .
Numerical Simulations . . . . . . . . . . . . . .
4.6.1 Simulation Parameters and Procedure .
4.6.2 Estimation of Drag Coeﬃcients . . . . .
4.6.3 Simulation Results . . . . . . . . . . . .
Validation . . . . . . . . . . . . . . . . . . . . .

5 Experimental Studies and Control
5.1 Introduction . . . . . . . . . . . . . .
5.2 Experiments . . . . . . . . . . . . . .
5.2.1 SPI Performance . . . . . . .
5.3 Eﬀect of Time-Varying Internal Fluid
5.3.1 Additional Terms . . . . . . .
5.3.2 Numerical Simulations . . . .

. . . . .
. . . . .
. . . . .
Velocity
. . . . .
. . . . .

.
.
.
.
.
.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

51
52
55
55
60
64
66
66
66
67
68
69
69
70
70
71
72
73
73
74
75
79
81
81
86
88
88
88
91
96

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

99
99
101
103
105
105
108

6 Concluding Remarks and
Future Work
113
6.1 Dynamics of a Pipe Conveying Fluid with a Non-Uniform Velocity Proﬁle . . 113
6.2 Simpliﬁed Dynamics of a Fluid-Conveying, Fluid-Immersed Pipe Aﬃxed to a
Rigid Body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
vi

6.3
6.4

General Dynamics of a Fluid-Conveying,
Fluid-Immersed Pipe Aﬃxed to a Rigid Body . . . . . . . . . . . . . . . . . 116
Experimental Studies and Control . . . . . . . . . . . . . . . . . . . . . . . . 118

Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

121

LIST OF TABLES

4.1

Names of diﬀerent ﬂuid control volumes and control surfaces. Note that there
is no need for a “head internal outlet” and “tail internal inlet”, since these
surfaces are the same. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

4.2

Parameters used in simulations . . . . . . . . . . . . . . . . . . . . . . . . .

90

4.3

Diﬀerent tail geometries and their forward speeds and steady-state oscillation
amplitudes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

5.1

Mean speeds of the three conﬁgurations in experiments . . . . . . . . . . . . 104

viii

LIST OF FIGURES

2.1

A ﬂuid-conveying pipe and a magniﬁed view of a small length element . . . .

7

2.2

Free-body diagram of (a) a ﬂuid element and (b) its corresponding pipe element

7

2.3

A cross-sectional view of the three ﬂuid volumes of a triple plug ﬂow

. . . .

9

2.4

Plot of the momentum ﬂux correction factor µ as a function of the Reynolds
number. The data points obtained through simulation are shown by circles. .

15

(a) Locus of the ﬁrst three roots of Det (Z) = 0 for uniform (dashed line) and
non-uniform (solid line) ﬂow models with β = 0.3 (b) A magniﬁed image of
the second root loci in (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

Plot of (a) ucr and (b) ωcr for uniform (dashed line) and non-uniform (solid
line) turbulent ﬂow for diﬀerent values of β. . . . . . . . . . . . . . . . . . .

19

Left: Modes two and three of the uniform proﬁle case prior to mode switch,
β = 0.385. Right: Closeup of boxed region. The marked points are the values
of each locus at u = 8.498, and are at points (23.75, 5.64) and (24.46, 6.25) for
the second and third mode, respectively. . . . . . . . . . . . . . . . . . . . .

19

Left: Modes two and three of the uniform proﬁle case after mode switching
has occurred, β = 0.387. Right: Closeup of boxed region. The marked points
are the values of each locus at u = 8.502, and are at points (23.60, 6.08) and
(24.59, 5.92) for the second and third mode, respectively. . . . . . . . . . . .

20

Closeup of second mode, where β = 0.290, 0.297, 0.305 in the direction of
the plotted arrow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

2.10 Plot of (a) ucr and (b) ωcr for uniform (thin line) and laminar (thick line)
ﬂow for diﬀerent values of β. . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

2.5

2.6

2.7

2.8

2.9

3.1

A cantilevered ﬂuid-conveying pipe, with a magniﬁed view of a small length
element. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ix

25

3.2

The cross-section of a ﬁnned tube with ﬂuid-conduit diameter D and tail span
S. The dotted circle describes the area responsible for the added external mass
Me , which is equal to 0.25ρf πS 2 , where ρf is the density of the external ﬂuid. 26

3.3

The proposed submersible, comprised of a rigid body and a ﬂuid-conveying
ﬂexible tail. Note that the tail is assumed to have a ﬁnned-tube cross section,
as shown in Fig.3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

3.4

Free-body diagrams of the rigid body and tail of the submersible in Fig.3.3. .

33

3.5

Argand diagram for the ﬁrst three modes of oscillation for βi = 0.01, βe = 0.9,
λ = 1/2, ψB = 1/12, µ = 2.25 and ue = 1.0. . . . . . . . . . . . . . . . . . .

36

3.6

Illustration of the procedure used to ﬁnd the neutral stability curve. . . . . .

37

3.7

Neutral stability curves for diﬀerent values of µ. . . . . . . . . . . . . . . . .

39

3.8

Plot showing the relationship between ωcr and ue for diﬀerent values of µ. .

40

3.9

The darkened region of each neutral stability curve (left of the dotted lines)
depicts the region of negative thrust. . . . . . . . . . . . . . . . . . . . . . .

41

3.10 Eﬃciency for neutrally stable points in the ui -ue space and for diﬀerent values
of µ. Note that the curve for µ = 2.5 is unique in that it has an interim region
of zero eﬃciency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

3.11 Eﬃciency as a function of ui and ue for the mass fraction µ = 2.5. . . . . . .

44

4.1

4.2

4.3

A submersible with a ﬂuid-conveying ﬂexible tail. The submersible is propelled
by the combined action of the ﬂuid jet and ﬂuttering tail. In this ﬁgure, the
deﬂection of the tail has been exaggerated. . . . . . . . . . . . . . . . . . . .

53

The cross-section of (a) the rigid head of the submersible at some point along
its length, and (b) the ﬂexible tail of the submersible, which is a ﬁnned tube.
In both ﬁgures, the dotted circle describes the area responsible for the added
ﬂuid mass in the y direction. . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

Depiction of control volumes Ωh , Ωt , and Ωi , i = 1, 2, 3, 4; and control surfaces
Γj , j = 1, 2, 3, 4, 5, 6, along with their outward normals. . . . . . . . . . . . .

56

x

4.4

Free-body diagrams of (a) a diﬀerential element of Ω2 (tail internal ﬂuid), (b)
a diﬀerential element of Ωt (ﬂexible tail), and (c) a diﬀerential element of Ω4
(tail external ﬂuid). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

4.5

A ﬂuid-conveying beam in the form of a ﬁnned-tube with ﬁn height h. . . . .

89

4.6

The subplot on the left shows the forward speed of the submersible with ﬂexible tail, rigid tail, and rigid tube conﬁgurations. On the right, the displacement of the tail tip y(L) and orientation θ of the rigid head of the submersible
are given for the ﬂexible tail conﬁguration. . . . . . . . . . . . . . . . . . . .

92

Sequence of images of the submersible with the ﬂexible tail. The scaling for
the y axis has been adjusted to make the tail oscillations more visible. . . . .

94

The subplot on the left shows the forward speed of the submersible with
L = 0.30, 0.34 and 0.38 where h = 0.077. On the right, the displacement of
the tail tip y(L) is shown for these tail conﬁgurations. . . . . . . . . . . . . .

95

The subplot on the left shows the forward speed of the submersible with
h = 0.067, 0.072 and 0.077 where L = 0.34. On the right, the deﬂection of
the tail tip y(L) is shown for these tail conﬁgurations. . . . . . . . . . . . . .

96

4.10 Displacement of the tail tip y(L) in mm for a cantilever, U = Ucr = 2.372
m/s. Neutral stability is evident, since the oscillations of the tail are neither
growing nor shrinking over time. . . . . . . . . . . . . . . . . . . . . . . . . .

97

4.7

4.8

4.9

5.1

Exterior of the current prototype SPI, with marked features of interest. For
interpretation of the references to color in this and all other ﬁgures, the reader
is referred to the electronic version of this dissertation. . . . . . . . . . . . . 101

5.2

High-speed images of the submersible with the ﬂexible tail acquired over a
period of 1 sec, soon after it started from rest in the MSU swimming pool.
The average speed of the submersible in these images was ≈ 0.4 BL/s, less
than the maximum speed of ≈ 1.15 BL/s. Note the large deﬂections of the
tail which occur during acceleration. . . . . . . . . . . . . . . . . . . . . . . 103

5.3

Mean speeds of the conﬁgurations listed in Table 5.1. The error bars express
the 95% conﬁdence intervals. . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

5.4

Comparison of the constant-velocity U = 4.0 case to non-constant velocity
case where α = 40, k = 0.75 . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
xi

5.5

Comparison of three values of k: 0.25, 0.75 and 1.0. The scaling parameter
α = 40 for all cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

5.6

An image sequence of the submersible during a turn. Note that the dashed
line is collinear with the non-inertial x axis. . . . . . . . . . . . . . . . . . . 112

xii

Chapter 1
Motivation
Fish-like propulsion has been a matter of academic interest since Gray’s pioneering work
in the 1930’s; the oft-cited “Gray’s Paradox” is derived from a 1936 paper [14] in which
the speed attained by a dolphin was calculated to require approximately seven times the
power available to the animal. Though later workers resolved this apparent paradox (see
the work by Fish [12]), interest in the eﬃciency of ﬁsh-like motion has been sustained to the
present day. An early model of oscillatory propulsion Taylor [42] was based on calculation
of the resistive force applied by the surrounding ﬂuid; this model, which neglects inertial
eﬀects, is more appropriate to the study of low-Reynolds number propulsion, such as that
of spermatozoa. These inertial eﬀects were accounted for by Lighthill in his seminal 1960
paper [23] which used slender-body theory to approximate the eﬀect of the pressure ﬁeld
surrounding the ﬁsh. Lighthill found that a travelling waveform with higher amplitude at
the tail than the head is required to produce eﬃcient thrust. Later works have extended
these initial eﬀorts to account for a planform of variable height [44], deﬂections of arbitrary
amplitude [26], and to account for the wing-like “lunate” tail which is a feature of the fastest
1

carangiform1 swimmers [11], along with a host of other works.
The continued drive within the academic community to “chase the tail” is in large part
due to impressive estimates of eﬃciency. For example, Lighthill’s slender body methods
estimate that 90% of the power expended by a swimming ﬁsh is available to propel the ﬁsh,
and that the remainder is wasted, serving only to raise the kinetic energy in the wake. This
value of 90% “Froude” eﬃciency lines up well with a ﬁgure of 87% found by Anderson and
coworkers [2] found on a simpliﬁed experimental setup. Some simpliﬁcations were required
to obtain these ﬁgures; it is diﬃcult to precisely measure the eﬃciency of a self-propelled
swimming body, since the sources of drag and thrust cannot be separated. An interesting
paper by Schultz and Webb [37] describes a variety of approaches which have been employed
to determine drag and thrust, but they conclude that forward speed is the most robust way
of determining whether a tail waveform is superior to another for a given platform.
This interest has manifested itself in a number of robotic platforms designed to harness
the eﬃciency of ﬁsh-like motion. The MIT RoboTuna [5, 43] is a well known “biomimetic”
platform. That device employs a tail composed of several links covered by a Neoprene
sheath. These links were then individually actuated in an approximation of a ﬁshes motion.
A diﬀerent type of linked system was built by McMasters and coworkers [27] which used a
novel geared mechanism to produce a ﬁsh-like motion using only a single servomotor, while a
third approach was taken by Alvarado and Youcef-Toumi [1], who described a system which
used a wholly ﬂexible tail, excited with a single actuator at the base. These platforms are
simply some of the more complex and recent devices; experimental devices used to measure
1: “Carangiform” swimming is typiﬁed by large amplitude motion concentrated in the
posterior half of a swimming ﬁshes body. The term is used in contrast to “anguilliform”
motion, in which large amplitude deﬂections take place over the length of the creature’s body.
Lighthill’s excellent review paper [24] describes these and other propulsion mechanisms in
more detail.
2

ﬁsh propulsion were built as early as 1933 [13]. Other devices have also been built to
investigate single aspects of ﬁsh swimming, such as the caudal ﬁn work of Anderson and
coworkers [2], and studies of pectoral ﬁn swimming undertaken by Lauder and coworkers
[22].
A fundamental diﬀerence between the approaches taken by Triantafyllou [43] and Alvarado and Youcef-Tuomi [1] is that the RoboTuna’s need to independently actuate its
linked tail means that the system must essentially be driven against its natural dynamics.
In contrast, the ﬂexible tail described in [1] can be induced to oscillate near its natural
frequency. This approach can lead to reduced transmission losses and lower motor bandwidth requirements, and the mechanism so produced can easily be made very robust, since
it consists of a single ﬂexible element. The present work proposes a sort of ﬂexible tail that
is fundamentally diﬀerent from that described by Alvarado. Speciﬁcally, the desired traveling waveform is produced not by an actuator at the base, but instead by ﬂutter instability
induced in a slender tail by conveying ﬂuid down its center. This instability is sometimes
referred to as the “garden hose instability”, after the tendency of such a hose to whip about
when the faucet is opened too far.
The ﬂuid-conveying pipes and ﬂutter instability are well-represented in the academic
literature - the lists of papers cited in Paidoussis’ review books on the topic [30, 31] are a
testament to this. However, the literature is largely grounded in a relatively small number of
applications, such as pipeline vibrations [3]. Therefore, even a simpliﬁed model of a ﬂuttering
ﬂuid-conveying submersible requires expansion of the literature.
This work is organized as follows. Chapter 2 derives the equations of motion for a ﬂuidconveying pipe, and presents a simple way to account for a non-uniform velocity proﬁle found
3

in real pipe ﬂow. Chapter 3 presents a simple, analytically tractable model of a ﬂuttering
ﬂuid-conveying submersible, consisting of a rigid body aﬃxed to one end of a ﬂuid-conveying
tail. This rigid body boundary condition, while not particularly unusual, had never been
presented in the ﬂuid-conveying pipe literature. Chapter 3 also presents some estimates of
thrust and eﬃciency for a ﬂuid-conveying tail’s waveform based on the results of Lighthill
[23] and Wu [44].
Chapter 4 derives a considerably more complex model of the submersible, a model which
is able to account for the general motion of a ﬂuttering ﬂuid-conveying submersible within
a plane. As with the rigid body boundary condition, it is interesting that despite the
well-developed state of the literature, the equations of motion for a ﬂuid-conveying pipe
in a non-inertial reference frame have not been presented. Chapter 5 discusses previous
designs of ﬁsh-like submersibles and presents preliminary experimental data from a platform
designed at Michigan State University. The equations of motion derived in Chapter 4 are
also extended to account for a time-varying velocity of the conveyed ﬂuid and are solved with
an eye toward maneuvering the submersible. Conclusions and thoughts for the direction of
future research on this topic are provided in Chapter 6.

4

Chapter 2
Dynamics of a Pipe Conveying Fluid
with a Non-Uniform Velocity Proﬁle

2.1

Introduction

Analytical treatments of the ﬂuid-conveying pipe problem typically invoke the “plug ﬂow”
assumption, under which the velocity proﬁle is uniform across the cross-section of the pipe.
At high Reynolds number, this assumption is approximately correct, since the velocity proﬁle
is nearly uniform over the central region of the cross-section, with only a thin, highly sheared
annular region near the pipe wall. The plug ﬂow assumption is less valid at low Reynolds
numbers, since a larger fraction of the ﬂuid at a given cross section has signiﬁcantly diﬀerent
momentum than the average ﬂuid element at that cross section. The deviation from the plug
ﬂow assumption is particularly apparent for laminar ﬂow. This chapter presents a method
by which the boundaries of ﬂutter instability of pipes conveying low Reynolds number ﬂuid
may be determined.
5

This chapter is organized as follows. The equation of motion of a pipe conveying ﬂuid
with a uniform velocity proﬁle (plug ﬂow) is ﬁrst presented in Section 2.2 as background
material. The derivation of plug ﬂow is extended to a non-uniform velocity proﬁle in Section
2.3. A ﬂow model composed of three concentric plugs of ﬂuid derived ﬁrst and extended to
N plugs; large values of N can be used to model an arbitrary velocity proﬁle. The resulting
equation for the non-uniform ﬂow diﬀers from the standard equation for uniform ﬂow by a
single constant; it is highly tractable, requiring no more eﬀort to solve than the standard
equation. This constant, a momentum ﬂux correction factor similar to that used in many
ﬂuid mechanics texts [33], is determined in Section 2.4. The equations of motion for a
cantilever pipe are solved and the results for the uniform and non-uniform ﬂow equations
are compared in Section 2.5 for both laminar and turbulent ﬂow.

2.2

Background - Uniform Velocity Proﬁle

Consider the ﬂuid-conveying pipe in Fig.2.1 where U denotes the ﬂuid velocity (constant)
relative to the pipe. Assuming an Euler-Bernoulli beam model of the pipe, its equation of
motion [30] can be written as follows

EI

∂2y
∂2y
∂2y
∂4y
+ MU 2
+ 2MU
= 0,
+ (M + m)
∂x∂t
∂x4
∂x2
∂t2

(2.1)

where E and I are the Young’s modulus of elasticity and area moment of inertia of the
cross-section of the pipe, and M and m are the ﬂuid mass and pipe mass per unit length.
The above equation can be obtained by showing that the rate of change of linear momentum
6

y

r

U

y(x,t)

U

R

x

x

dx
dx

Figure 2.1: A ﬂuid-conveying pipe and a magniﬁed view of a small length element
of a ﬂuid element is given by the expression [30]
dL
∂
∂ 2
y dx ˆ
j
=M
+U
dt
∂t
∂x

and from the equations of motion of a ﬂuid element and its corresponding pipe element
Fig.2.2 in the x and y directions

q S dx

∂p
A(p + 1 dx)
2 ∂x

F dx

y

1 ∂p
A(p − 2 dx)
∂x

x

q S dx
dx

∂Q
Q+ 1
2 ∂x dx
M + 1 ∂M dx
q S dx
2 ∂x
F dx q S dx
T + 1 ∂T dx
2 ∂x

M − 1 ∂M dx
2 ∂x
1 ∂Q dx
Q− 2
∂x
1
T − 2 ∂T dx
∂x

(a)

(b)

Figure 2.2: Free-body diagram of (a) a ﬂuid element and (b) its corresponding pipe element

Plug ﬂuid element:

∂p
∂y
− qS − F
= 0,
∂x
∂x
∂ ∂y
∂y
∂ 2
∂
F − A (p ) − qS
y,
=M
+U
∂x ∂x
∂x
∂t
∂x
−A

7

(2.2)
(2.3)

Pipe element:
∂T
∂y
∂2y
= 0,
+ qS + F
−Q
∂x
∂x
∂x2
∂Q
∂
∂y
∂y
∂2 y
−F +
(T ) + qS
=m
.
∂x
∂x ∂x
∂x
∂t2

(2.4)
(2.5)

In Eqs.(2.2) through (2.5), A denotes the cross-sectional area of the pipe, S denotes its
internal surface area per unit length, p denotes the ﬂuid pressure, q denotes the shear stress
in the ﬂuid, and F denotes the force per unit length normal to the wall. For the pipe element,
Q, M and T denote the traverse shear force, bending moment, and tension, respectively.

2.3
2.3.1

Non-Uniform Velocity Proﬁle
Triple plug ﬂow model

The triple plug ﬂow model assumes three concentric volumes of ﬂuid being conveyed through
the pipe; the cross-sectional area of these volumes are shown in Fig.2.3. The ﬂow velocity is
diﬀerent for the three volumes but is assumed to be constant within each volume. The triple
plug ﬂow is not physically realizable but its analysis provides the framework for investigation
of a general velocity proﬁle. Note that the ﬂuid volume at the center (marked 1 in Fig.2.3)
has a single ﬂuid-ﬂuid interface; the volume in the middle (marked 2 in Fig.2.3) has two
ﬂuid-ﬂuid interfaces; and the outermost ﬂuid volume (marked 3 in Fig.2.3) has one ﬂuidﬂuid and one ﬂuid-pipe interface. The analysis for a general velocity proﬁle will require us
to introduce more volumes with two ﬂuid-ﬂuid interfaces, similar to volume 2. It is for this
reason that three “plugs” are required; using fewer does not give rise to volumes with two
8

ﬂuid-ﬂuid interfaces, and using more creates redundant elements.

y

z
1
2
3

Figure 2.3: A cross-sectional view of the three ﬂuid volumes of a triple plug ﬂow

To extend the analysis of plug ﬂow in section 2.2 to triple plug ﬂow, we ﬁrst denote
the cross-sectional areas of volumes 1, 2, and 3 as A1 , A2 and A3 respectively; their ﬂow
velocities as U1 , U2 and U3 respectively; and their mass per unit length as M1 , M2 and
M3 respectively. F12 , F23 and F3p denote the radial force per unit length between ﬂuid
volumes 1 and 2, ﬂuid volumes 2 and 3, and ﬂuid volume 3 and the pipe, respectively. The
shear force at these interfaces are denoted by q12 , q23 and q3p respectively, and the surface
area per unit length of these interfaces are denoted by S12 , S23 and S3p respectively. The
dynamics of the ﬂuid volumes 1, 2 and 3 can now be replicated from Eqs.(2.2) and (2.3),
and that of the pipe from Eqs.(2.4) and (2.5), as follows
Volume 1:

∂y
∂p
− q12 S12 − F12
= 0,
−A1
∂x
∂x
∂2y
∂y
∂
∂ 2
F12 − pA1
− q12 S12
y,
= M1
+ U2
∂x
∂t
∂x
∂x2
9

(2.6)
(2.7)

Volume 2:

∂p
∂y
−A2
− q23 S23 + q12 S12 − (F23 − F12 )
= 0,
∂x
∂x
∂
∂ 2
∂2y
∂y
y,
= M2
+ U2
(F23 − F12 ) − pA2
− (q23 S23 − q12 S12 )
∂x
∂t
∂x
∂x2

(2.8)
(2.9)

Volume 3:

∂p
∂y
−A3
− q3p S3p + q23 S23 − (F3p − F23 )
= 0,
∂x
∂x
∂y
∂2y
∂
∂ 2
(F3p − F23 ) − pA3
− (q3p S3p − q23 S23 )
y
= M3
+ U3
∂x
∂t
∂x
∂x2

(2.10)
(2.11)

Pipe:
∂T
∂2 y
∂y
= 0,
+ q3p S3p + F3p
−Q
∂x
∂x
∂x2
∂y
∂Q
∂2y
∂2y
+ q3p S3p
.
− F3p + T
=m
∂x
∂x
∂x2
∂t2

(2.12)
(2.13)

The summation of equations in the x direction, namely Eqs.(2.6), (2.8), (2.10) and (2.12),
gives
∂p ∂T
−(A1 + A2 + A3 )
+
=0
∂x ∂x

Note that (A1 + A2 + A3 ) = A, where A is the inner cross-sectional area of the pipe. The
above equation simpliﬁes to
∂
(T − pA) = 0
∂x
if A = A(x), the same result as that which is obtained for plug ﬂow. The summation of
10

equations in the y direction, Eqs.(2.7), (2.9), (2.11) and (2.13), gives
3
∂
∂ 2 y ∂Q
∂2y
∂2y
∂ 2
+
=
Mn
−pA
y+m
+T
+ Un
∂x
∂t
∂x
∂x2
∂x2
∂t2
n=1

Following the derivation of plug ﬂow [30] and substituting (T − pA) = 0 and Q = −EI

∂3y
∂x3

into the above equation, the following expression is obtained:

EI

∂4y



3



∂2y





3

∂2y



3



∂2y
2
+
Mn Un 
+2 
Mn Un 
= 0. (2.14)
+ m +
Mn 
∂x∂t
∂x4
∂x2
∂t2
n=1
n=1
n=1

2.3.2

N-plug ﬂow model

It is simple to envision the analysis of the preceding section with more volumes. The additional volumes, similar to Volume 2 of the triple plug ﬂow in that they possess two ﬂuidﬂuid interfaces, do not complicate the analysis since all interfacial terms are cancelled. Thus
Eq.(2.14) can be rewritten as

EI



N



∂2y





N

∂2y



N



∂2y
2
Mn Un 
Mn Un 
Mn 
+
+2
= 0, (2.15)
+ m +
∂xt
∂x4
∂x2
∂t2
n=1
n=1
n=1
∂4y

where N is any integer greater than three. For very large values of N, the volumes have
inﬁnitesimal thickness and the summations in the coeﬃcients of Eq.(2.15) can be replaced
with integrals. The average velocity of ﬂuid over a cross section is deﬁned as

1
¯
U=
A

A

U(A) dA =

for a cylindrical pipe of radius R.
11

R
2
U(r) rdr
R2 0

The coeﬃcients of the dynamical equation in Eq.(2.15) can now be expressed as



N



2
Mn Un  = ρf

U 2 (A) dA = 2πρf

R

¯
U 2 (r) rdr = µM U 2 ,

A
0
n=1


N
¯
¯

2
Mn Un  = 2ρf
U(A) dA = 2ρf AU = 2M U ,
A
n=1


N
m +
Mn  = m + M,
n=1

(2.16)

where ρf is the density of the conveyed ﬂuid, and µ is the non-dimensional momentum ﬂux
correction factor. For a cylindrical pipe with a known ﬂuid velocity proﬁle U(r), µ has the
expression
µ=

R U(r) 2
2
rdr.
¯
U
R2 0

(2.17)

Using the algebraic simpliﬁcations in Eq.(2.16), Eq.(2.15) can be rewritten as follows

EI

2
2
2
∂4y
¯ ∂ y + 2M U ∂ y + (m + M) ∂ y = 0.
¯
+ µM U 2
∂x∂t
∂x4
∂x2
∂t2

(2.18)

Note the similarity of Eq.(2.18) with that of plug ﬂow given by Eq.(2.1). All terms are
essentially identical with the exception of the additional constant µ, which is a function of
the velocity proﬁle. Note also that evaluation of Eq.(2.17) for a uniform velocity proﬁle
yields µ = 1, meaning that Eq.(2.1) is consistent with its plug ﬂow assumption under this
expanded derivation.
12

2.4

Analysis of a Cantilever Pipe

2.4.1

Solution of the diﬀerential equation

The behavior of a ﬂuid-conveying pipe for cantilever boundary conditions has been analyzed,
namely

y(0, t) = 0,

∂y(0, t)
(0, t) = 0,
∂x

∂ 3 y(L, t)
= 0,
∂x3

∂ 2 y(L, t)
= 0,
∂x2

(2.19)

The equation of motion is made non-dimensional with the following change of variables

Y =

y
,
L

X=

x
,
L

T = t Ω.

By introducing the non-dimensional velocity, mass fraction and frequency as follows

u=

M 1/2 ¯
U L,
EI

β=

M
,
m+M

ω=

M + m 1/2
ΩL2 ,
EI

and assuming a separable form for Y (X, T ) such that

Y (X, T ) = φ(X) e−iωT

the following non-dimensional equation of motion and boundary conditions are obtained:
d4 φ
dφ
d2 φ
+ µu2
+ 2β 1/2 uiω
− ω2φ = 0
4
2
dX
dX
dX
∂ 2 φ(1)
∂φ(0)
= 0,
(0, t) = 0,
φ(0) = 0,
∂X
∂X 2
13

∂ 3 φ(1)
=0
∂X 3

(2.20)

The solution of φ is assumed to be of the form φ(X) = AezX and this yields the characteristic
polynomial
z 4 + µu2 z 2 + 2β 1/2 uiωz − ω 2 = 0.

(2.21)

For speciﬁc values of µ, u and β, Eq.(2.21) provides four roots, zn , n = 1, 2, 3, 4, where
zn = zn (ω). The complete solution of φ(X) has the form

φ(X) = A1 ez1 X + A2 ez2 X + A3 ez3 X + A4 ez4 X

The solution of the equation above based on the boundary conditions in Eq.(2.20) results in
the complete solution
4
Y (X, T ) =
n=1

An ezn X eiωT =

4

An eRe[zn ]X ei(Im[zn ]X+Re[ω]T ) e−Im[ω]T
n=1
(i)
(ii)
(iii)

An inspection of the above equation indicates that Y (X, T ) is a product of three exponential
terms of which the ﬁrst term is bounded since X is bounded and the second term is oscillatory
since the exponent is imaginary. The third term can grow unbounded with time if Im[ω] < 0
and this represents unstable dynamics of the pipe. The exact mode and velocity at which
the ﬂuid-conveying pipe becomes unstable depends on the ﬂuid mass fraction β.

2.4.2

Determination of momentum ﬂux correction factor

For laminar ﬂow, the Poiseuille solution of Navier-Stokes equation holds [33] and the value
of µ can be analytically determined to be equal to 4/3 for a circular pipe. For turbulent
ﬂow, the value of µ approaches unity as the value of Reynolds number approaches inﬁnity.
14

The literature [7] commonly cites single values of µ for turbulent ﬂow based on the boundary
layer model employed, and the values cited are derived using assumptions better suited to
high-Reynolds number turbulence. Since the present work relaxes the assumption of highReynolds number turbulence, a model for µ such that µ = µ(Re) is desired. Numerical values
of µ for turbulent ﬂow were calculated based on Eq.(2.17) using velocity proﬁles generated
by the commercial software STAR-CCM. Even with this improved resolution, knowledge of µ
is required for more values of the Reynolds number than are feasible to simulate. To predict
the value of µ for laminar, transition, and turbulent ﬂow regimes, the following curve-ﬁt was
employed:



 4/3
Re ≤ 2200




µ(Re) =
 3.647 − 0.001052 × Re 2200 < Re < 2413




 1.04 + 167.2/Re
Re ≥ 2413.

(2.22)

1.4
laminar

µ (Re)

1.3
transition

1.2

1.1
4

turbulent

1.0
0

x 10

2

4

6

Re
Figure 2.4: Plot of the momentum ﬂux correction factor µ as a function of the Reynolds
number. The data points obtained through simulation are shown by circles.
It can be seen from Fig.2.4 that the data points obtained from simulation matches well with
the expression of µ in Eq.(2.22). The choice of 2200 < Re < 2413 to deﬁne the “transition
15

region” is somewhat arbitrary, as might be the choice of using one single value for µ to
represent a phenomenon as complex as turbulent transition. However, the authors note that
for a given ﬂow setup the range of velocities corresponding to transition ﬂow is small, making
the choice of the transition model relatively unimportant1 .
For a circular pipe, expressions can be determined which relate the non-dimensional
velocity u and the Reynolds number:

u=

M 1/2 ¯
UL,
EI

Re =

¯
UD
ν

⇒

u = νL Re

πρf
.
4EI

(2.23)

Clearly, a value of u does not uniquely determine the Reynolds number. The geometric
properties of the pipe (L, E, I) and density and kinematic viscosity of the working ﬂuid (ρf ,
ν) should be speciﬁed to determine the Reynolds number for a given u. Since the momentum
ﬂux correction factor µ is a function of the Reynolds number, it follows that the value of µ
cannot be determined from the value of u alone.

2.5

Model Comparison:
Uniform and Non-Uniform Flow

2.5.1

Turbulent ﬂow

The locus of the ﬁrst three roots of Det (Z) = 0 is shown in the Argand diagram in Fig.2.5(a)
for β = 0.308 for both uniform and non-uniform ﬂow models. Since the ﬂuid-conveying pipe
undergoes ﬂutter instability in the second mode, a separate Argand diagram of the loci of the
1: The laminar to turbulent transition region can be seen in the Argand diagram in
Fig.2.5.
16

4

40

ucr = 8.13
ωcr = 23.06

transition

first root

20

Im (ω)

Im (ω)

2
third root
u=0

0

second root

0

20

40

u=0

0

−2

60

ucr = 6.94
ωcr = 14.46

15

20

Re (ω)

Re (ω)

(a)

25

(b)

Figure 2.5: (a) Locus of the ﬁrst three roots of Det (Z) = 0 for uniform (dashed line) and
non-uniform (solid line) ﬂow models with β = 0.3 (b) A magniﬁed image of the second root
loci in (a).

second root is shown in Fig.2.5(b). The root loci for the uniform ﬂow model are a function
of u alone (µ is implicitly assumed to be unity) but they are a function of both u and µ for
the non-uniform ﬂow model. As noted in the previous section, the value of µ is not uniquely
deﬁned in terms of u. Certain dimensional coeﬃcients related to the working ﬂuid and pipe
geometry must be assumed to obtain this relationship. The non-uniform ﬂow model assumes
water to be the working ﬂuid (ρf = 1000 kg/m3 , ν = 1.0 × 10−6 m2 s) and the following
parameters for the pipe:

E = 1.7 MP a,

I = 6.48 × 10−10 m4 ,

L = 0.5 m

It is acknowledged that the need to specify dimensional parameters is a limitation but this
limitation is not signiﬁcant. An inspection of Fig.2.4 reveals that µ is weakly related to
17

the Reynolds number for turbulent ﬂow and therefore dependence of µ on the dimensional
parameters is not signiﬁcant.
It can be seen from Fig.2.5 that the root loci for the uniform and non-uniform ﬂow models
are quite diﬀerent though their µ values are quite similar1 . A close look at the second root
loci indicates that the uniform ﬂow model predicts ﬂutter instability of the pipe to occur
for a critical velocity of ucr = 8.13 with ωcr = 23.06 whereas the non-uniform ﬂow model
predicts signiﬁcantly lower values of ucr = 6.94 (15% lower) and ωcr = 14.46 (37% lower).
Clearly, the dynamics of the system are very sensitive to the value of µ in the neighborhood
of β = 0.3. The values of ucr and ωcr are plotted in Fig.2.6 for diﬀerent values of β. This
ﬁgure indicates that the non-uniform ﬂow model predicts signiﬁcantly lower values of ucr and
ωcr for β in the neighborhood of 0.7 as well. There is good agreement between the uniform
and non-uniform turbulent ﬂow models for values of β that are not in the neighborhood of
0.3 or 0.7.
The system’s behavior in the regions near β = 0.3 and β = 0.7 is somewhat diﬀerent than
what has appeared in the literature [15], since we are examining only the lowest velocity at
which the pipe becomes unstable. The source of the system’s sensitivity to changes in β is
not a mode-switching phenomenon [38]. A mode switch for the uniform ﬂow system does
occur near β = 0.386, and this phenomenon is depicted in Figs.2.7 and 2.8.
Rather, the discontinuous behavior is a consequence of the unusual shape of the second
mode’s locus in this region of the β parameter space. Fig.2.9 shows the shape of second
mode’s locus in this sensitive region of the β-space. The β = 0.297 curve is nearly tangent
to the real axis near ω = 16 and the β = 0.290, 0.305 curves ﬁrst cross the real axis at
1: The value of µ is implicitly assumed to be unity for the uniform ﬂow model whereas
it has values in the neighborhood of 1.05 for the non-uniform ﬂow model.
18

16

50
43.87

40

12.85

12
ωcr

ucr

10.52

30
26.90

8.28

8

21.88

20

6.80

14.15

4

0

0.3

0.5

0.7

10

1

0

0.3

0.5

β

1

β

(a)

0.7

(b)

Figure 2.6: Plot of (a) ucr and (b) ωcr for uniform (dashed line) and non-uniform (solid
line) turbulent ﬂow for diﬀerent values of β.

20
7
6.5

15

Im(ω)

Im(ω)

6
10

5.5
5

5

4.5
0
0

20

40

4
20

60

Re(ω)

22

24
Re(ω)

26

28

Figure 2.7: Left: Modes two and three of the uniform proﬁle case prior to mode switch,
β = 0.385. Right: Closeup of boxed region. The marked points are the values of each locus
at u = 8.498, and are at points (23.75, 5.64) and (24.46, 6.25) for the second and third mode,
respectively.

19

20

7

18
16

6.5
14

6
Im(ω)

Im(ω)

12
10
8
6

5.5
5

4

4.5

2
0
0

10

20

30

40

50

4
20

60

Re(ω)

22

24
Re(ω)

26

28

Figure 2.8: Left: Modes two and three of the uniform proﬁle case after mode switching has
occurred, β = 0.387. Right: Closeup of boxed region. The marked points are the values of
each locus at u = 8.502, and are at points (23.60, 6.08) and (24.59, 5.92) for the second and
third mode, respectively.
ω = 14.35 and ω = 22.54, respectively.

2.5.2

Laminar ﬂow

In prior work, the plug ﬂow model implicitly made the assumption that µ = 1. It is clear from
Fig.2.4 that this assumption is reasonable only for high Reynolds number. The non-unique
relationship between u and Re allows situations where this assumption is not reasonable.
For example, inspection of Eq.(2.23) reveals that a suﬃciently long pipe could have a large
value for u at low Reynolds number. A suﬃciently long ﬂuid-conveying pipe could therefore
undergo ﬂutter instability with laminar ﬂow, in which case µ = 4/3. In the case of turbulent
ﬂow, the near-unity value of µ yields values of ucr and ωcr which diﬀer signiﬁcantly from the
uniform case only in certain thin regions of the β parameter space. Figure 2.10 is analogous
20

4

3

Im(ω)

2

1

0

−1

14

16

18
20
Re(ω)

22

24

Figure 2.9: Closeup of second mode, where β = 0.290, 0.297, 0.305 in the direction of the
plotted arrow.
to Fig.2.6, but assumes laminar ﬂow and gives a very diﬀerent result. The values of ucr
for laminar ﬂow are quite diﬀerent from that of uniform ﬂow in all regions of the parameter
space. The values of ωcr for laminar ﬂow are similar to that of uniform ﬂow over much of the
range of β but there are certain regions where the predicted values are signiﬁcantly diﬀerent.
Unlike the thin regions near β = 0.3 and β = 0.7 for turbulent ﬂow (see Fig.2.6), the regions
of large separation for laminar ﬂow extend from approximately β = 0.29 to β = 0.40 and
from β = 0.69 to β = 0.93.
This dramatic diﬀerence between the solution found using the “standard” plug ﬂow assumption and the improved model presented in this communication is perhaps the best
argument for the standard use of this improved model, particularly at scales outside of
typical engineering practice.

21

40
ωcr

50

12
ucr

16

8

0.29
0.40

30
20

4
2

0

0.5

10

1

0.69

0

0.5

β

1

β

(a)

0.93

(b)

Figure 2.10: Plot of (a) ucr and (b) ωcr for uniform (thin line) and laminar (thick line) ﬂow
for diﬀerent values of β.

22

Chapter 3
Simpliﬁed Dynamics of a
Fluid-Conveying, Fluid-Immersed
Pipe Aﬃxed to a Rigid Body

3.1

Introduction

In the 1970s, Paidoussis and coworkers built and tested [28] a propulsor consisting of a ﬂuidconveying pipe aﬃxed to the underside of a large surface hull. That work has limited theoretical justiﬁcation, but the superﬁcial similarity between a ﬂuttering ﬂuid-conveying cantilever
and a ﬁsh is readily apparent — both motions take the form of a traveling wave which grows
in amplitude from tip to tail. It is therefore easy to appreciate that the well-documented efﬁciency of ﬁsh-like motion [14, 25] was the motivation for this type of mechanism. Although
the ﬂuttering motion of the tail was found to be a net gain above the thrust provided by the
jet exhausting into the water, the eﬃciency did not approach that of a propeller, and after
23

a patent [29], the idea appears to have been dropped.
Interim developments since these initial eﬀorts justify a reexamination of the ﬂuttering
ﬂuid-conveying tail as a propulsor. Improvements in battery and motor technologies now
permit the prime mover and power source to be packaged into the neutrally buoyant hull of a
small submersible, whereas Paidoussis’ surface hull was very large relative to the propulsor.
Reducing the size of the hull changes the system’s dynamics, since the ﬂuid-conveying tail
may produce suﬃcient transverse force and moments to alter the motions of the hull. That
is, the boundary conditions of the tail are likely to be signiﬁcantly diﬀerent than those of
a cantilever. The rigid body boundary condition developed in this chapter and used in its
investigations is a close linear approximation to the boundary condition acting on the base
of a ﬂuttering, ﬂuid-conveying submersible’s tail.
This chapter is organized as follows. Section 3.2 provides the background for the dynamic equations of ﬂuid conveying, ﬂuid immersed pipes and the mechanics of slender body
swimming. Section 3.3 presents the rigid body boundary conditions and makes them nondimensional. The analytical method used for solution is also describe in this section. Section
3.4 investigates the ﬂow requirements for controlled ﬂutter and calculations of thrust and
eﬃciency of the waveforms produced by the ﬂuttering tail.

3.2
3.2.1

Background
Fluid-conveying pipes

The equations of motion and boundary conditions for a cantilever pipe conveying ﬂuid with
constant velocity Ui , immersed in an inviscid ﬂuid ﬂowing with constant velocity Ue , and
24

ignoring gravitational, viscous, pressurization and tensile eﬀects, are as follows [31]:

EI

2
2
2
∂4y
2
2 ∂ y + 2(MU + M U ) ∂ y + (m + M + M ) ∂ y = 0, (3.1)
+ (MUi + Me Ue )
e e
e
i
∂x∂t
∂x4
∂x2
∂t2

y(0, t) = 0,

∂y
(0, t) = 0,
∂x

∂3y
(L, t) = 0
∂x3

∂2y
(L, t) = 0,
∂x2

y
Ui

r

Ui

y(x,t)
x

R

x
dx
dx

Figure 3.1: A cantilevered ﬂuid-conveying pipe, with a magniﬁed view of a small length
element.
where y(x, t) is the displacement of the pipe, as shown in Fig.3.1, E, I and L denote the
Young’s modulus, area moment of inertia, and length of the pipe, respectively, and m, M
and Me represent the mass per unit length of the beam, the internal (conveyed) ﬂuid, and
the external ﬂuid. The masses per unit length of the beam and internal ﬂuid can be easily
computed but the mass per unit length of the external ﬂuid requires approximation. One
method for computation of Me is to use the added mass coeﬃcient [9]. For thin cross sections,
such as that of a ﬂat plate, the added mass is equal to the mass of water within the cylinder
which circumscribes the plate cross-section. A “ﬁnned tube”, well suited to providing both a
ﬂuid conduit and a tail of adequate span, is depicted in Fig.3.2 with the area responsible for
the added external mass marked. Viscous terms characteristic of the external ﬂow [17] and
internal ﬂow [20] are of course needed to compute the overall eﬃciency of the submersible,
25

D
S

Figure 3.2: The cross-section of a ﬁnned tube with ﬂuid-conduit diameter D and tail span
S. The dotted circle describes the area responsible for the added external mass Me , which
is equal to 0.25ρf πS 2 , where ρf is the density of the external ﬂuid.
and aﬀect beam stability as well. These eﬀects are neglected in this section to permit an
analytical solution of the equations, and will be considered in the following section along
with the additional accelerations and nonlinear rotations of the rigid body.
Equation (3.1) may be made non-dimensional via the following change of variables

X=

x
,
L

Y =

y
,
L

T =

t
L2

1/2
EI
.
m + M + Me

(3.2)

The non-dimensional velocities ui and ue and the mass fractions βi and βe are deﬁned as
follows
ui =
ue =

M 1/2
Ui L,
EI
Me 1/2
Ue L,
EI

M
βi =
,
m + M + Me
Me
,
βe =
m + M + Me

then Eq.(3.1) can be written in its non-dimensional form
∂2Y
∂4Y
+ (u2 + u2 )
+ 2(ui
e
i
∂X 4
∂X 2

βi + ue

βe )

∂2Y
∂2Y
=0
+
∂X∂T
∂T 2

A separable form for y(x, t) is assumed, i.e., y(x, t) = f (x)eiΩt , such that the non-dimensional
26

variable Y takes the form

Y (X, T ) = φ(X)eiωt,

ω=

m + M + Me 1/2
ΩL2 .
EI

(3.3)

Separation yields the ordinary diﬀerential equation and boundary conditions
∂4φ
∂2φ
+ (u2 + u2 )
+ 2(ui
e
i
∂X 4
∂X 2

φ(0) = 0,

∂φ(0)
= 0,
∂X

βi + ue

βe ) iω

∂ 2 φ(1)
= 0,
∂X 2

∂φ
− ω 2 φ = 0,
∂X

(3.4)

∂ 3 φ(1)
= 0,
∂X 3

The solution of φ is assumed to be of the form φ(X) = AezX . For speciﬁc values of ui , ue , βi
and βe , the characteristic polynomial of Eq.(3.4) provides four roots zn , where zn = zn (ω),
n = 1, 2, 3, 4. The solution of φ(X) therefore takes the form

φ(X) = A1 ez1 X + A2 ez2 X + A3 ez3 X + A4 ez4 X .

(3.5)

Substitution of Eq.3.5 into Eq.3.4 yields the identity


1
1
1
 1


 z
z2
z3
z4
1


 2 z
 z e 1 z 2 ez2 z 2 ez3 z 2 ez4
2
3
4
 1

3
3
3
3
z1 ez1 z2 ez2 z3 ez3 z4 ez4
Z









  A1   0 
  

  

 A   0 
 2   
 =  .

  

  A3   0 
  

  

0
A4

(3.6)

A non-trivial solution (for ω) of Eq.3.6 is obtained by numerical evaluation of the roots of
Det(Z) = 0. This equation has inﬁnite roots in ω. Substitution of Eq.(3.5) into Eq.(3.3)
27

yields
4
Y (X, T ) =
n=1

An ezn X eiωT =

4

An eℜ[zn ]X ei(ℜ[zn ]X+ℜ[ω]T ) e−ℜ[ω]T .
n=1
(i)
(ii)
(iii)

(3.7)

It can be seen that Y (X, T ) is a product of three exponential terms of which the ﬁrst term
is bounded (since X is bounded), and the second term is periodic since the exponent is
imaginary. The third term can grow unbounded with time if ℜ[ω] < 0 and this represents
the onset of ﬂutter instability. The mode and velocity at which the pipe becomes unstable
depends on the ﬂuid mass fractions βi and βe . The coeﬃcients An , n = 1, 2, 3, 4, can
be computed from the nullspace of the matrix Z in Eq.3.6, once ω and zn , n = 1, 2, 3, 4,
have been determined. These coeﬃcients are needed to estimate the force exerted by the
ﬂuid-conveying tube on the surrounding ﬂuid.

3.2.2

Slender body swimming

Fish-like propulsion has been a topic of interest in the academic community for more than 60
years and several robotic platforms have been built [43, 27] to exploit the phenomenon. The
mechanism proposed in this communication (see Fig.3.3) is composed of a ﬂuttering ﬂuidconveying tail providing thrust by both jet and tail action; the tail has the cross-sectional
proﬁle of a ﬁnned tube like the one shown in Fig.3.2. A similar mechanism was constructed
in the 1970s [28] and was found to produce positive thrust only if the phase velocity of the
tail displacement was greater than the forward speed of the vessel.
Thrust production via a high phase velocity traveling wave was described ﬁrst in a paper
by Lighthill [23], which used slender body analysis to approximate the thrust produced by
28

ℓ

Ue

Ui

y
L

x

Figure 3.3: The proposed submersible, comprised of a rigid body and a ﬂuid-conveying
ﬂexible tail. Note that the tail is assumed to have a ﬁnned-tube cross section, as shown in
Fig.3.2.

an idealized ﬁsh. Lighthill found that a traveling waveform, for example

y(x, t) = f (x) cos(kx + Ωt)

(3.8)

can produce positive thrust if (Ω/k) > Ue , where Ue is the speed of the body relative to the
external ﬂuid. The quantity Ω/k is known as the phase velocity. The following equation is
the dimensional form of Eq.(3.7) for a single waveform

yn (x, t) = An eℜ[Zn ]x ei(ℜ[Zn ]x+ℜ[Ω]t) e−iℜ[Ω]t ,

(3.9)

where Zn and Ω are the dimensional wavenumber and frequency. The waveform described
by Eq.(3.9) will result in positive thrust if

ℜ[Ω]
> Ue .
ℜ[Zn ]
29

The above equation can be made non-dimensional, taking the form

ℜ[ω]
L2

1/2 L
EI
> ue
m + M + Me
ℜ[zn ]

EI 1/2 1
Me
L

⇒

ℜ[ω]
ue
>
.
ℜ[zn ]
β 1/2

(3.10)

Equation (3.10) gives us a condition under which a waveform will generate positive thrust,
where ℜ[ω]/ℜ[zn ] is the non-dimensional phase velocity of the waveform.
Paidoussis’ described his icthyoid propulsor [28] as having a single phase velocity, which
was measured by direct observation. While a single phase velocity is relatively simple to determine experimentally, determination of positive thrust is not straightforward in the context
of Eq.(3.7) since it has four traveling waveforms of diﬀerent, spatially variable amplitudes
and phase velocities. In the context of Eq.(3.7), it is easier to estimate thrust by the method
laid out by Lighthill [23] and Wu [44]. In those papers, a slender1 ﬁsh is considered and the
time-averaged thrust τ is given by the relation

1
τ = Me
2

∂y 2
∂y 2
− Ue
∂t
∂x

−
x=L

∂y 2
∂y 2
− Ue
∂t
∂x

,

(3.11)

x=0

where y = y(x, t) denotes the displacement of the slender body from its neutral position, y
˙
and y ′ denote partial derivatives of y with respect to time and x, and overbar refers to a longterm time average. In the above equation, increasing Me increases the thrust generated; this
leads to the design choice of the ﬁnned-tube in Fig.3.2. The form of Eq.(3.11) also makes
it clear that a higher forward speed Ue requires a higher velocity y to sustain it. Both
˙
Lighthill (Lighthill, 1960) and Wu (Wu, 1971) derived the expression in Eq.(3.11) without
the assumption of harmonic motion. If harmonic motion is assumed, the time average over a
1: A deﬁnition of “slender body” can be found in §2 of Wu [44].
30

single cycle is suﬃcient. Note that Eq.(3.11) diﬀers slightly from that found in (Wu, 1971)1 .
Wu’s assumption that no mass is aﬀected at x = 0 is relaxed, since we have assumed that
the tail has a geometry which is uniform along its length rather than a tapered ﬁsh which
has zero area at the tip.

A similar expression (Lighthill, 1960) was derived for the average power P required to
provide the displacements y(x, t)

P = Ue Me

∂y
∂t

∂y
∂y
+ Ue
∂t
∂x

−
x=L

∂y
∂t

∂y
∂y
+ Ue
∂t
∂x

,

(3.12)

x=0

which includes the power lost in the vortex wake. Equations (3.11) and (3.12) indicate that a
large magnitude of y produces greater thrust but also requires higher energy input. This can
˙
also be veriﬁed from the expression of the Froude eﬃciency (Lighthill, 1960) of the motion
of the slender body

η=

τ Ue
.
P

(3.13)

The expression for eﬃciency in Eq.(3.13) does not account for power lost to internal ﬂuid
shearing in the pipe, external drag, or the thrust produced by the ﬂuid expelled from the
ﬂuid-conveying tail. The actual eﬃciency of the submersible will therefore be somewhat
lower after these eﬀects are accounted for.

1: Equation (47) in that work.
31

3.3

Fluid-conveying pipe aﬃxed to a rigid body

3.3.1

Simpliﬁed boundary conditions

The complete planar dynamics of a submersible with a ﬂuid-conveying ﬂuttering tail, as
shown in Fig.3.3, will be presented in the next chapter. However, a simpliﬁed model that
permits better understanding of the important parameters governing the motion is desirable.
To this end, we make the following assumptions:

A1. The rigid body is symmetric about the plane containing the neutral surface of the
undeformed beam.

A2. The rotation of the rigid body, denoted by θ in Fig.3.4, is small. This is in addition
to the Euler-Bernoulli beam assumption that the slope of the tail is small everywhere
along its length.

A3. The submersible has zero acceleration in the x direction. This allows us to ignore the
“forces” that arise from a non-inertial reference frame.

A4. The added mass coeﬃcient associated with the rigid body is zero. This assumption
simpliﬁes the analysis by allowing us to concentrate on the geometry of the rigid body
in the xy plane. This is equivalent to the assumption that the rigid body is a planar
lamina in the xy plane.
32

ℓ

V
M

M

θ

F

F

y
V

MB , JB

x

Figure 3.4: Free-body diagrams of the rigid body and tail of the submersible in Fig.3.3.
The above assumptions allow us to write the boundary conditions for the ﬂuid-conveying
tail1 at x = 0 as follows:

EI
EI

∂2y
∂x2

∂3y
∂2y
−ℓ
∂t2
∂x∂t2

∂3y
+ MB
∂x3

− (JB + MB ℓ2 )

∂3y
∂x∂t2

+ MB ℓ

= 0,

(3.14a)

= 0,

(3.14b)

x=0

∂2y
∂t2 x=0

where MB is the mass of the rigid body, and JB is the mass moment of inertia of the rigid
body about its center of mass, and ℓ is the distance of the center of mass of the rigid body
from the base of the tail. The boundary conditions in Eq.(3.14) can be derived from the
free-body diagram of the rigid body in Fig.3.4 as follows

−V = MB

∂2y
∂3y
,
−ℓ
∂t2
∂x∂t2 x=0

M − ℓV = JB

∂3y
∂x∂t2

.
x=0

It should be pointed out that the variable θ in Fig.3.4 denotes the orientation of the rigid
body, which is equal to the slope of the tail at x = 0, i.e., θ = [∂y/∂x]x=0 . Since the rigid
body is symmetric about the plane containing the neutral surface of the undeformed beam,
1: The result given here is similar to that in [8] with the exception that they have been
derived here at x = 0 rather than x = L.
33

the force F does not aﬀect the boundary conditions.

Using Eq.(3.2), we can obtain the non-dimensional form of Eq.(3.14):
∂3Y
∂2Y
∂3Y
+ µ(
−λ
)
= 0,
∂X 3
∂T 2
∂X∂T 2 X=0
∂3Y
∂2Y
∂2Y
= 0,
− µ (ψB + λ2 )
−λ
∂X 2
∂X∂T 2
∂T 2
X=0

(3.15a)
(3.15b)

where
µ=

MB
,
(m + M + Me )L

λ=

ℓ
,
L

JB
,
ψB =
MB L2

are non-dimensional parameters of the rigid body. Physically, µ is the ratio of its mass to the
mass of the rest of the system, λ is a non-dimensional distance, and ψB is the square of the
non-dimensional radius of gyration. Note that Eq.(3.15) provides the boundary conditions
for a free end as µ → 0, and that of a clamped end as µ → ∞. The free-end conditions can
be shown easily whereas the clamped-end conditions can be shown by allowing µ to approach
∞ in Eq.(3.15), which yields:
∂3Y
∂2Y
= 0,
−λ
∂T 2
∂X∂T 2 X=0
∂3Y
∂2Y
(ψB + λ2 )
= 0.
−λ
∂X∂T 2
∂T 2 X=0

(3.16a)
(3.16b)

Since ψB = 0, it can be readily shown from the above that ∂ 2 Y /∂T 2 = ∂ 3 Y /∂X∂T 2 = 0
at X = 0. Assuming zero initial velocities, i.e., ∂Y (0, 0)/∂T = ∂ 2 Y (0, 0)/∂X∂T = 0, we
obtain the clamped end conditions ∂Y (0, T )/∂T = ∂ 2 Y (0, T )/∂X∂T = 0.
34

3.3.2

Method of analysis

The simpliﬁed boundary conditions for a ﬂuid-conveying pipe aﬃxed to a rigid body are
investigated in the same manner as that of a cantilever pipe, which was discussed in section
2.1. Using the boundary conditions in Eq.(3.15), we get the following relation that is similar
to Eq.(3.6):


η2
η3
η4
 η1


 ζ
ζ2
ζ3
ζ4
1


 2 z
 z e 1 z 2 ez2 z 2 ez3 z 2 ez4
2
3
4
 1

3
3
3
3
z1 ez1 z2 ez2 z3 ez3 z4 ez4









  A1   0 

  

  
 A   0 
 2   

 =  ,

  
  A3   0 

  

  
A4
0

(3.17)

where ηn , ζn , n = 1, 2, 3, 4, are deﬁned by the relations
2
ηn = zn + µω 2 (ψB + λ2 )zn − λ ,
3
ζn = zn − µω 2 (1 − λzn ).
Equation (3.17) leads to a solution of the form given by Eq.(3.7). Since the proposed application of the swimming submersible requires the oscillations of the ﬂuttering tail not to
grow with time, the points of neutral stability, i.e. ℜ[ω] = 0, are sought in the ui -ue space.
For a given µ, λ, ψB , βi , βe , and a given forward velocity ue , the value of ui for which
ℜ[ω] = 0 can be found through analysis of an Argand diagram, an example of which is
provided in Fig.3.5. This diagram is constructed [15] by determining the natural frequencies
of a set of modes at ui = 0, then gradually incrementing ui to determine ℜ[ω] and ℜ[ω]
for higher values of ui . The natural frequencies at the onset of ﬂutter instability, ωcr , are
the locations where the resultant curves cross the imaginary axis. The velocities at which
35

this occurs are referred to as critical velocities, ucr . This procedure yields a single neutrally
stable point in the ui -ue space. Clearly, while it is possible to search the entire ui -ue space
for these neutrally stable points, this would require construction of a large number of Argand
diagrams to obtain an acceptable resolution in ue , which is prohibitively time-expensive. To
mitigate this problem, an automated method similar in character to that discussed above is
proposed.

20
first
root

10
I (ω)

third root
second
root

0

ωcr = 14.64
ucr = 5.24

0

20

40

60

ℜ (ω)
Figure 3.5: Argand diagram for the ﬁrst three modes of oscillation for βi = 0.01, βe = 0.9,
λ = 1/2, ψB = 1/12, µ = 2.25 and ue = 1.0.
In the automated method, the value of ui required for neutral stability at ue = 0 is ﬁrst
computed by interpolating values of ui for which ℜ[ω] ≈ 0. This is repeated once more for
ue = ǫ, where ǫ is a small number. Subsequent points may be found in the following manner,
which is explained with the help of Fig.3.6. For a set of two consecutive neutrally stable
points already determined, such as the points marked 1 and 2 in Fig.3.6, the next point is
guessed to exist at point 3′ at a distance ǫ1 from point 2 along the vector v drawn from
36

neutral stability curve

ǫ2

ue

3′
ǫ1

3
2
v

1
ui
Figure 3.6: Illustration of the procedure used to ﬁnd the neutral stability curve.

point 1 to point 2. Several points at which to compute ω are then chosen near point 3′ along
a vector perpendicular to v. These points are separated from each other by the distance ǫ2 .
The value for ωcr may now be found by interpolation, and computation of the critical values
of ui an ue follows trivially. This procedure, in which the direction of iteration is nearly
perpendicular to the curve, was necessary to navigate some of the sharp turns in the neutral
stability curves given in the next section. In general, lower values of ǫ1 and ǫ2 are needed
for curves with sharper turns, and in the current work, ǫ1 = 0.02, ǫ2 = 0.002 were found
to suﬃce. The method was found to be robust as well; while we present results for a single
rigid body of variable mass here, changes in other parameters can be easily accomodated.
37

3.4
3.4.1

Stability, thrust and eﬃciency
Neutral stability

We obtained neutral stability curves in the ui -ue space for a ﬂuid-conveying pipe aﬃxed
to a rigid body, using the simpliﬁed boundary conditions discussed in the previous section.
These neutral stability curves were obtained for various values of µ. The values of the other
parameters were chosen as follows:

βi = 0.01,

βe = 0.9,

λ = 1/2,

ψB = 1/12.

The values of βi and βe used here reﬂect the ﬁnned-tube geometry, shown in Fig.3.2. The
nature of this geometry is such that a large amount of external ﬂuid is associated with the
oscillations of the tail and this explains the relatively large value of βe . Likewise, the small
value of βi reﬂects the small area through which internal ﬂuid is conveyed. The values of λ
and ψB correspond to a uniform cylinder with length equal to that of the tail.
The neutral stability curves are shown in Fig.3.7. The area inside each curve (towards
the origin) represents the region where the tail does not ﬂutter whereas the area outside
represents the region where the tail ﬂutters. The diﬀerence between the low-µ and high-µ
curves is striking. For values of ue < 2.9, the value of ui required to create ﬂutter decreases
with increase in the value of µ. The conﬂuence of all the curves at ue ≈ 2.9 is also interesting
though no theoretical reason for this conﬂuence is apparent. It should be mentioned that the
curves’ tendency to pass close to one another is not a unique scenario. Recall that the rigid
body can be described by the parameters ψJ , λ, and µ, of which we have only investigated
the latter in detail in this work. Changing the values of ψJ and λ does not, in our experience,
38

15
2.5

10 2.25

5

25

ue

2.0

∞

1.0

0.33

5
ui ≈ 5.7, ue ≈ 2.9

0

0

2

4

6

8

ui
Figure 3.7: Neutral stability curves for diﬀerent values of µ.
remove this tendency of the curves to intersect at a point, though the point in ui -ue space
is diﬀerent.
For values of ue > 2.9, higher values of µ are neutrally stable at a higher value of ue for
a given value of ui . This implies that for a given ﬂow rate provided by the prime mover, a
tail aﬃxed to a hull of larger mass will tend to ﬂutter at higher forward speed. The curves
also indicate that higher values of µ are neutrally stable at a higher value of ui for a given
value of ue . This implies that for a given external ﬂow, a tail aﬃxed to a hull of larger mass
will require higher internal ﬂow to ﬂutter.
Figure 3.8 plots the relationship between ue and ωcr . It can be seen from this ﬁgure
that for high mass ratios, there exist certain regions of the ue parameter space where ωcr
is sensitive to ue . This sensitivity may be used to determine a “sweet spot”, where a
higher oscillation frequency may be obtained with a minimal change in ue . This increase
in oscillation frequency and concomitant increase in y will increase the thrust produced by
˙
39

50
40
25

30

∞

ωcr

5 2.5
2.25

20
10

2.0
0.33

0

0

5

10

15

ue
Figure 3.8: Plot showing the relationship between ωcr and ue for diﬀerent values of µ.
the beam, per Eq.(3.11), and can be used for designing high-acceleration maneuvers for the
submersible. Thrust production and eﬃciency in the context of Eqs.(3.11) and (3.12) is
discussed in the next two sections.

3.4.2

Thrust characteristics

To determine the thrust produced by the ﬂuttering tail of the proposed submersible, we
non-dimensionalize Eq.(3.11) and compute the average over one cycle as follows:

τ∗ =

τ L2
ωcr 2π/ωcr
=
EI
4π 0

2
˙
βe Y 2 − u2 Y ′
e

X=1

2
˙
− βe Y 2 − u2 Y ′
e

X=0

dT. (3.18)

The function Y (X, T ) is found by the method presented in section 2.1, and takes the form of
Eq.(3.7). Equation (3.7) has both real and imaginary parts; only the real part is physically
manifest and contributes to the thrust. Assuming neutral stability (ℜ[ω] = 0), the real part
40

of Eq.(3.7) is
4
Y (X, T ) =

eℜ[zn ]X ℜ[An ] cos(ℜ[zn ]X + ℜ[ω]T ) − ℜ[An ] sin(ℜ[zn ]X + ℜ[ω]T )

n=1
The coeﬃcients An in the above equation are found by computing the nullspace of the matrix
˙
in Eq.(3.17). The terms Y and Y ′ in Eq.(3.18) can then be obtained through diﬀerentiation.

15
2.5
2.25

10

5

25

ue

2.0

∞

1.0

0.33

5
ui ≈ 5.7, ue ≈ 2.9

0

0

2

4

6

8

ui
Figure 3.9: The darkened region of each neutral stability curve (left of the dotted lines)
depicts the region of negative thrust.
Figure 3.9 reproduces the curves of neutral stability in Fig.3.7 with dark lines depicting
the region on each curve where the thrust is negative. Notice that no value of µ allows
thrust-producing ﬂutter instability at ui = 0. This matches with our physical intuition that
a ﬂapping ﬂag cannot generate thrust. At low values of ui , the positive hydrodynamic work
is predominantly contributed by the external ﬂuid and reduces the energy of that ﬂuid; this
phenomenon can be put to use in power generation [41]. It is also interesting to note that
systems with lower values of µ can produce thrust at lower values of ui than systems with
41

higher values of µ, though the high-µ systems have higher forward speed ue .
It is important to remember, however, that merely having positive thrust from the tail
does not guarantee that a given ui -ue point can be reached. The system’s drag and the
thrust of the ﬂuid jet will also govern the submersible’s top speed. Since the drag of the
system will, for a neutrally-buoyant vessel, be strongly related to the displacement and mass,
we will reserve these concerns for a later work more closely tied to the physical realization
of the submersible.
Finally, it should be mentioned that the calculations presented here are purely concerned
with the thrust produced by the traveling waveform in the ﬂuttering tail. As such, the
thrust produced by the ﬂuid jet is not considered. However, it has been shown through
experiments [28] that this loss of thrust can be overcome such that the combined jet and tail
action produces a net thrust that is higher than that of a ﬁxed jet.

3.4.3

Hydrodynamic eﬃciency

Similar to the expression for thrust, Eq.(3.12) can be made non-dimensional and the average
over one cycle computed, to give the average non-dimensional power P ∗ :
1/2 3
∗ = P Me L =
P
(EI)3/2
ωcr 2π/ωcr
1/2 ˙
1/2 ˙
˙
˙
βe ue Y 2 − u2 βe Y ′ Y
− βe ue Y 2 − u2 βe Y ′ Y
e
e
2π 0
X=1
X=0

dT

From Eq.(3.13), the expression for Froude eﬃciency can be written as

η=

τ Ue
τ ∗ ue
=
.
P∗
P
42

(3.19)

0.6
∞
1.0

0.4

25

η

5

0.2
2.25

2.5

0.33 2.0

0.0

0

5

10

15

ue
Figure 3.10: Eﬃciency for neutrally stable points in the ui -ue space and for diﬀerent values
of µ. Note that the curve for µ = 2.5 is unique in that it has an interim region of zero
eﬃciency.

Eﬃciencies computed using Eq.(3.19) are plotted with respect to ue for various values of µ
in Fig.3.10. Each curve plots the eﬃciency for neutrally stable points in the ui -ue space.
Per the discussion in [44], Eq.(3.19) has meaning only when the thrust is positive, and
therefore, Fig.3.10 plots the eﬃciency only for neutrally stable points with positive thrust,
the non-darkened regions of Fig.3.9.
It is interesting to note that the maximum eﬃciency is insensitive to the value of µ.
This will provide ﬂexibility in submersible design since the mass of the hull can be chosen
based on other factors such as power source, drag and buoyancy, rather than hydrodynamic
eﬃciency. For each value of µ, the eﬃciency remains near its peak value for a wide range
of ue . This trend resembles the eﬃciency curves for tail-swimming ﬁsh which are known to
maintain eﬃciency over a broad range of swimming speeds. Incidentally, this broad peak in
eﬃciency is not characteristic of typical marine propellers, which tend to be most eﬃcient
43

over a narrow range of velocities.
It should be mentioned that the values of eﬃciency cited here consider only the power
used to generate tail motion in the external pressure ﬁeld and the thrust provided by the tail.
Other factors, such as pipe losses and viscous drag will aﬀect the eﬃciency of the vehicle
when considered as a whole. However, the thrust provided by the ﬂuid jet has also not been
considered; a full accounting of the system’s eﬃciency would require further investigation.

0.6

η

0.4
8
0.2
0.0
15

4 ui
10

5

ue

0
0

Figure 3.11: Eﬃciency as a function of ui and ue for the mass fraction µ = 2.5.
A combined perspective on Figs.3.9 and 3.10 is presented in Fig.3.11, which depicts the
neutral stability curve for µ = 2.5 in three dimensions, with the eﬃciency shown out of the
plane of the page. It is easy to see from this ﬁgure that the eﬃciency of tail motion for a
propulsor of this type is dependent on both ue and ui . Since a value of zero on the above
chart represents a net zero or negative thrust, it should be pointed out that high values of ue
require high values of ui in order for the tail to oscillate with a thrust-generating waveform.
44

This result could be reasonably predicted; the limit case of a cantilever perturbed only by
an axially ﬂowing external ﬂuid (ui = 0) is not expected to add momentum in the axial
direction to the surrounding ﬂuid. Nonetheless, it is good to see this expectation conﬁrmed
by analysis.

45

Chapter 4
General Dynamics of a
Fluid-Conveying, Fluid-Immersed
Pipe Aﬃxed to a Rigid Body

4.1

Nomenclature

A large number of symbols have been employed in this chapter. A nomenclature has therefore
been provided for the reader’s convenience.
h

Height of the ﬂexible tail in the z direction, [m]

ˆˆ k
i, j, ˆ

Unit vectors along the x, y and z axes, respectively

xx
kh , k

xy yx yy
,k ,k
h
h
h

Added mass coeﬃcients associated with accelerations of the rigid head of
the submersible in the x and y directions, [−]

xx xy yx yy
kt , kt , kt , kt

Added mass coeﬃcients associated with accelerations of the ﬂexible tail
of the submersible in the x and y directions, [−]
46

θθ θθ
kh , kt

Added mass coeﬃcients associated with angular acceleration of the rigid
head and ﬂexible tail of the submersible, [−]

mh , mt

Mass of the rigid head and ﬂexible tail of the submersible, [kg]

mi

Mass of the ﬂuid within control volume Ωi , [kg]

n
ˆ

The outward-pointing unit vector normal to a surface

nj
ˆ

The outward-pointing unit vector normal to a control surface Γj

p2 (x), p4 (x) Pressure in ﬂuid control volumes Ω2 and Ω4 , [N/m2 ]
q

Shear stress at the wall of the ﬂuid-conveying tube [N/m2 ]

r

Position vector of a point in the xy reference frame, [m]

˙
rxy

Time derivative of r as seen by an observer in the xy frame, [m/s]

¨
rxy

˙
Time derivative of rxy as seen by an observer in the xy frame, [m/s2 ]

rc

Position vector of the center-of-mass of the system comprised of the submersible and ﬂuid control volumes in the xy reference frame, [m]

rcx , rcy

x and y components of rc , [m]

˙
rcxy

Time derivative of rc as seen by an observer in the xy frame, [m/s]

rcx , rcy
˙
˙

˙
x and y components of rcxy , [m/s]

¨
rcxy

˙
Time derivative of rcxy as seen by an observer in the xy frame, [m/s2 ]

rcx , rcy
¨
¨

¨
x and y components of rcxy , [m/s2 ]

xy-frame

Body-ﬁxed reference frame as shown in Fig.4.1

x, y
¯ ¯

x and y coordinates of the center-of-mass of control volume Ω, [m]

xh , yh
¯ ¯

x and y coordinates of the center-of-mass of the rigid head of the submersible, [m]

47

xt , yt
¯ ¯

x and y coordinates of the center-of-mass of the ﬂexible tail of the submersible, [m]

xi , yi
¯ ¯

x and y coordinates of the center-of-mass of the i-th ﬂuid control volume,
[m]

y(x, t)

Displacement of a point on the ﬂexible tail at a distance of x from the
origin of the xy-frame, [m]

yL

y(L, t), [m]

y(x, t)
˙

Time derivative of y(x, t) as seen by an observer in the xy-frame, [m/s]

yL
˙

y(L, t), [m/s]
˙

y ′ (x, t) First spatial derivative of y(x, t)
′
yL

y ′ (L, t)

Aj

Area of the control surface Γj , [m2 ]

ˆ ˆ ˆ
I, J, K

Unit vectors along the X, Y and Z axes, respectively

Fext

External force acting on the submersible, [N]

Fd

External force on the submersible per unit length due to drag, generated
by ﬂuid in the external control volumes, [N/m]

Fp

Pressure on the submersible due to non-ambient pressure at the entrance
of the internal tube, [N/m2 ]

FI

Force magnitude per unit length normal to the wall of the ﬂuid-conveying
tube of the ﬂexible tail due to internal ﬂow, [N/m]

FE

Force magnitude per unit length normal to the outer surface of the ﬂexible
tail due to external ﬂow, [N/m]

48

Jz

Mass moment of inertia of control volume Ω about the z-axis passing
through the origin of the xy-frame, [kg.m2 ]

Jzh , Jzt

Mass moment of inertia of the rigid head and ﬂexible tail of the submersible, respectively.

Both are computed about the z-axis passing

through the origin of the xy-frame, [kg.m2 ]
Jzi

Mass moment of inertia of the i-th ﬂuid control volume about the z-axis
passing through the origin of the xy-frame, [kg.m2 ]

L

Length of the ﬂexible tail in the undeformed conﬁguration, [m]

Lh

Length of the rigid head of the submersible along the x axis, [m]

M

Bending moment at a cross-section of the ﬂexible tail, [N.m]

Mext

Moment of Fext about the origin of the inertial reference frame, [N.m]

Ph , Pt

Perimeters of the head and tail, [m]

Q

Shear force magnitude at a cross-section of the ﬂexible tail, [N]

R

Position vector of a point in the XY frame, [m]

R0

Position vector of the origin of the xy reference frame in the XY frame,
[m]

R0x , R0y

x and y components of R0 , [m]

Rc

Position vector of the center-of-mass of the system comprised of the submersible and ﬂuid control volumes in the XY frame, [m]

S

Internal surface area per unit length of the ﬂuid conveying tube within
the ﬂexible tail, [m2 ]

T

Tension at a cross-section of the ﬂexible tail, [N]

49

U

Speed of the ﬂuid ﬂowing through the ﬂexible tail of the submersible,
[m/s]

Ue

Speed of the external ﬂuid in the x direction as seen by an observer in the
xy-frame, [m/s]

Uh

Speed of the ﬂuid ﬂowing through the rigid head of the submersible, [m/s]

V

Velocity of a particle in the XY reference frame, [m/s]

Vr

Velocity of a ﬂuid particle relative to a control surface, [m/s]

Vs

Velocity of a point on a control surface, [m/s]

XY -frame Inertial reference frame as shown in Fig.4.1
X0 , Y 0

X and Y components of R0 , [m]

ρ

Volumetric density, [kg/m3 ]

ρh , ρt

Volumetric density of the rigid head and ﬂexible tail of the submersible,
[kg/m3 ]

ρf

Volumetric density of the ﬂuid, [kg/m3 ]

θ

orientation of xy reference frame with respect to the XY frame, [rad]

Γ

Surface, [m2 ]

Γj

j-th control surface, [m2 ]

Ω

Volume, [m3 ]

Ωh , Ωt

Control volume of the rigid head and ﬂexible tail of the submersible, [m3 ]

Ωi

i-th ﬂuid control volume, [m3 ]

50

4.2

Introduction

In this chapter, we derive and solve the equations of motion for a submersible propelled by a
ﬂuttering ﬂuid-conveying tail while relaxing some of the assumptions imposed in the previous
chapter. In the previous chapter, it was assumed the rotations and transverse accelerations of
the rigid body were small, meaning that the equation of motion of the ﬂuid-conveying ﬂexible
tail could be stated within an inertial reference frame. In this more general formulation, the
submersible is assumed to have distinctly non-zero acceleration and rotation, which makes it
more convenient to write the equations of motion in a non-inertial frame of reference. Note
that the equation of motion for a ﬂuid-conveying pipe in a non-inertial frame has yet to
appear in the literature as of this writing.
It is perhaps interesting that the equation of motion for a ﬂuid-conveying pipe in a noninertial frame appears here for the ﬁrst time, given that an earlier attempt at this type of
propulsor was made by Paidoussis and coworkers in the mid-1970’s [28, 29] 1 . There are
two likely reasons for this. First, those works are primarily experimental in nature and
provide very little in the way of theory to begin with; as will be shown in this chapter, the
equations of motion in a non-inertial frame are not trivial, and do not lend themselves to
being “thrown in” to an experimental work. Also, while no dimensions have been published
for the surface ship-type hull described by [28], the ﬁgures provided in that work and the
related patent [29] imply that the hull is large relative to the tail. A suﬃciently large hull
is of course insensitive to the time-dependent transverse loads and moments imposed by a
ﬂuttering tail, and the testing method employed in those works was not sensitive to the
transient in forward acceleration.
1: An earlier cited work is an undergraduate thesis [36], which is diﬃcult to ﬁnd.
51

This chapter is organized as follows. Section 4.3 states the assumptions made for derivation of the equations of motion of the ﬂuid-conveying ﬂexible tail submersible. In section 4.4,
the derivation method is laid out, and the terms for the momentum and external forces are
described for the various components of the submersible. Section 4.5 describes the equations
of motion of a ﬂuid-conveying tail or pipe within a non-inertial reference frame. The simulation procedure used to solve the equations presented in sections 4.4 and 4.5 is presented
in section 4.6, along with the results of those simulations.
This chapter makes use of many variable names, only some of which will be recognizable
from prior sections - a nomenclature for this chapter is therefore provided in Appendix A.

4.3

Assumptions

Consider the ﬂuid-conveying submersible shown in Fig.4.1. It is propelled by the combined
action of a ﬂuid jet exiting from its ﬂexible tail and the ﬂuttering action of the tail. XY
denotes the inertial reference frame, xy is a body-ﬁxed reference frame, and θ is the orientation of the xy frame with respect to the inertial frame. The equations of motion of the
submersible are derived using the following simplifying assumptions:
A1. The submersible is comprised of a head, which is a rigid body, and a ﬂexible tail
clamped to the head. The submersible is neutrally buoyant and its motion is conﬁned
to the XY plane. The xy-frame is ﬁxed to the head at the base of the tail such that
the neutral axis of the undeformed tail is congruent with the x-axis. The head is
symmetric with respect to the xz and xy planes and is uniform along its length. In
its undeformed conﬁguration, the ﬂexible tail has a rectangular projection in the xz
plane. A corollary to this assumption is that both the head and tail have three planes
52

(x, y)

R

y

r

U
rc

Rc
θ

x

Y
Lh

R0

L

X
Figure 4.1: A submersible with a ﬂuid-conveying ﬂexible tail. The submersible is propelled
by the combined action of the ﬂuid jet and ﬂuttering tail. In this ﬁgure, the deﬂection of
the tail has been exaggerated.
of symmetry. The projection of the head and the undeformed tail in the xy plane is
slender.
A2. An internal cylindrical ﬂuid-conveying tube runs along the length of the rigid head and
the ﬂexible tail. The axis of the cylindrical tube is congruent with the neutral axis of
the ﬂexible tail and with the x-axis along the length of the head.
A3. The ﬂexible tail behaves like an Euler-Bernoulli beam. This implies that all points
along the length of the tail move in the y direction only, i.e., without foreshortening.
Furthermore, the displacement and slope of the tail with respect to the x-axis are
small.
A4. The ﬂow in the internal cylindrical tube is steady and has a uniform tube-axial velocity
proﬁle. Secondary ﬂows resulting from the curvature of the tube have been neglected.
53

A5. The ﬂuid conveyed through the internal cylindrical tube is the same as the ﬂuid surrounding the submersible. The submersible is far from any free or solid surfaces.

cross-section
of
rigid body

z

finned
tube
y

(a)

(b)

Figure 4.2: The cross-section of (a) the rigid head of the submersible at some point along its
length, and (b) the ﬂexible tail of the submersible, which is a ﬁnned tube. In both ﬁgures,
the dotted circle describes the area responsible for the added ﬂuid mass in the y direction.

The following assumptions are made in regards to the external ﬂuid:

A6. Geometry: The ﬂuid control volumes around the head and tail of the submersible,
which are deﬁned in the next section, have identical motion as that of the submersible.
The control volume around the head has identical axes of symmetry as those of the
head. The control volume around the tail is symmetric about the neutral axis of the
tail.

A7. Inertial Eﬀects: Each of the external control volumes contain the ﬂuid associated with
“added mass” [9]; note that the added-mass volumes may be diﬀerent for the x, y and
θ coordinates. A cross-section of the added mass volume for the y coordinate is shown
in Fig.4.2 for both the rigid head and the ﬂexible tail. The ﬂuid within these volumes
is inviscid.
54

A8. Drag Eﬀects: The “drag force” experienced by the head and tail of the submersible is
reﬂected within the surrounding control volumes of ﬂuid as a momentum deﬁcit within
that ﬂuid. Instead of approximating this momentum deﬁcit, the drag force on the body
will be computed directly, via separate analysis.

A9. Inlet Pressure: There is a non-ambient pressure at the entrance of the internal tube
which results from drawing surrounding ﬂuid into the head of the submersible. This
pressure will be computed via analysis of the captured streamtube.

The internal ﬂuid aﬀects the dynamics of the submersible in the following way:

A10. Inertial Eﬀects: The ﬂuid control volumes inside the head and tail of the submersible
are subjected to the same acceleration as that of the head and tail of the submersible,
respectively.

4.4
4.4.1

Dynamic Model - Newtonian Formulation
Force and Moment Equations

The equations of motion of the submersible are derived using the standard control volume
approach. The ﬂuid control volumes in and around the rigid head and ﬂexible tail of the
submersible are shown in Fig.4.3 along with their control surfaces. The names of the control
volumes and control surfaces are listed in Table 4.1. The equations presented in this section
are concerned with the motion of the xy-frame assuming a general function y(x, t) for the
motion of the tail. The tail’s dynamics are developed separately in a later section. The
55

translational motion of the xy-frame can be obtained from Newton’s second law of motion

d
d
Fﬂuid =
ρ V dΩ
ρh V dΩ +
dt Ω (t)
dt Ωt t
h

Section A-A

Γ3

Section B-B

(4.1)

Γ1

Section C-C

Section D-D

Γ5

Γ4

Γ6

Γ2

identical areas

Ωh

Ω1

Ωt

Ω2

Ω4

A

n2
ˆ
n6
ˆ

Ω3

B

n4
ˆ

D

D

C

n5
ˆ

n1
ˆ
n3
ˆ

C
B

A

Figure 4.3: Depiction of control volumes Ωh , Ωt , and Ωi , i = 1, 2, 3, 4; and control surfaces
Γj , j = 1, 2, 3, 4, 5, 6, along with their outward normals.
The ﬂuid forces on the left hand side of Eq.(4.1) are comprised of inertial forces associated
with the ﬂuid control volumes (assumptions A7 and A10), drag forces (assumption A8),
and the force generated by the non-ambient pressure at the entrance of the internal tube
56

Table 4.1: Names of diﬀerent ﬂuid control volumes and control surfaces. Note that there
is no need for a “head internal outlet” and “tail internal inlet”, since these surfaces are the
same.
Ωh
Ωt
Ω1
Ω2
Ω3
Ω4

Control Volumes
Rigid Head
Flexible Tail
Head Internal Fluid
Tail Internal Fluid
Head External Fluid
Tail External Fluid

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6

Control Surfaces
Jet Inlet
Jet Outlet
External Flow Head Inlet
External Flow Head Outlet
External Flow Tail Inlet
External Flow Tail Outlet

(assumption A9). The inertial forces are comprised of volume and surface terms obtained
from Reynolds’ transport theorem [33] and therefore the external forces can be written as
4

Fﬂuid = −
+

6





d

ρf V (Vr · nj )dΓ
ˆ
ρ V dΩ −
dt Ω (t) f
Γj (t)
i
i=1
j=1

Ωh ,Ωt

Fd dx +

Γ1 (t)

Fp dΓ

(4.2)

Substitution of the above equation into Eq.(4.1) yields

Ωh ,Ωt

Fd dx +

Γ1 (t)



4



d 
ρt V dΩ +
ρf V dΩ
ρ V dΩ +
dt Ω (t) h
Ωi (t)
Ωt
h
i=1


6

ρf V (Vr · nj )dΓ
ˆ
+
(4.3)
Γj (t)
j=1

Fp dΓ =

The volume integral terms of Eq.(4.3) may be simpliﬁed using Liebniz’ rule, written in
Eq.(4.4) for a general control volume Ω(t) with control surface Γ(t). Note that Liebniz’ rule
applies only to the volume integral term on the left-hand side of Eq.(4.4); the surface integral
which arises from the application of Liebniz rule can be combined with the surface integral
57

since Vs + Vr = V .
d
ρV dΩ +
ρV (Vr · n)dΓ
ˆ
dt Ω(t)
Γ(t)
=

∂
(ρV )dΩ +
ρV (Vs · n)dΓ +
ˆ
ρV (Vr · n)dΓ,
ˆ
Γ(t)
Γ(t)
Ω(t) ∂t

=

∂
(ρV )dΩ +
ρV (V · n)dΓ
ˆ
Γ(t)
Ω(t) ∂t

(4.4)

Using Liebniz’ rule as written above, Eq.(4.3) can be rewritten as follows

Ωh ,Ωt

Fd dx +

Γ1 (t)

4
∂
∂
∂
(ρh V )dΩ +
(ρt V )dΩ +
(ρf V )dΩ
Ωh ∂t
Ωi ∂t
Ωt ∂t
i=1
4
+
ρf V (V · nj )dΓ
ˆ
(4.5)
Γj
j=1

Fp dΓ =

An equation for the rotational motion of the submersible may be obtained by equating
the sum of external moments about the inertial origin to the rate of change of angular
momentum

d
Mext dΩ =
R × ρh V
dt Ω
Ωh ,Ωt
h

dΩ +

d
R × ρt V
dt Ωt

dΩ

(4.6)

˙
where V = R. Similar to the external forces, the external moments on the left hand side of
Eq.(4.6) are comprised of inertial moments, drag moments, and pressure moments. The inertial moments are comprised of volume and surface terms obtained from Reynolds’ transport
58

theorem [33] and therefore the external moments can be written as
4

6





d

Mext dΩ = −
R × ρf V dΩ −
R × ρf V (Vr · nj )dΓ
ˆ
dt Ω (t)
Ωh ,Ωt
Γj (t)
i
i=1
j=1
+

Ωh ,Ωt

R × Fd dx +

Γ1 (t)

R × Fp dΓ

where, as with Eq.(4.3), Liebniz’ rule may be applied to move the time derivative inside the
integral:
4

6

Mext dΩ = −
Ωh ,Ωt
i=1
+

Ωh ,Ωt





∂

R × ρf V dΩ −
R × ρf V (V · nj )dΓ
ˆ
Ωi (t)∂t
Γj (t)
j=1
R × Fd dx +

Γ1 (t)

R × Fp dΓ

(4.7)

Similar to the linear momentum equations, the equation above may be combined with
Eqs.(4.6) to form:

Ωh ,Ωt
4
i=1

R × Fd dx +

Ωi (t)

R×

Γ1 (t)

R × Fp dΓ =

∂
ρ V dΩ −
∂t f

6
j=1




Ωh

R×

Γj (t)

∂
ρ V
∂t h

dΩ +

Ωt


R × ρf V (V · nj )dΓ
ˆ

R×

∂
ρV
∂t t

dΩ

(4.8)

using the simpliﬁcation

∂
R × ρV
∂t

=R×

∂
ρV
∂t

+ V × ρV = R ×

∂
ρV
∂t

It is useful to rewrite Eq.(4.8) by taking moments and deﬁning the angular momentum about
the center of mass of the submersible and its attached control volumes. Using Eqs.(4.5) and
59

(4.8) yields the expressions

(R − Rc ) × Fd dx +
(R − Rc ) × Fp dΓ =
Γ1 (t)
Ωh ,Ωt
∂
∂
(R − Rc ) × ρh V dΩ +
(R − Rc ) × ρt V dΩ
Ωh ∂t
Ωt ∂t


6
4
∂

(R − Rc ) × ρf V dΩ +
(R − Rc ) × ρf V (V · nj )dΓ
ˆ
Γj (t)
Ωi (t) ∂t
j=1
i=1

(4.9)

4.4.2

Added Mass

When a body is accelerated through a ﬂuid, a volume of ﬂuid surrounding the body is
accelerated together with the body. The time rate of change of the momentum of this
surrounding ﬂuid results from a force to that ﬂuid by the body, and an equal and opposite
force is applied by that ﬂuid to the body. This force is known as the added mass eﬀect, and is
modeled as the product of the acceleration of the body and a directionally dependent “added
mass” related to the geometry of the body [9]. The directional dependence is reﬂected by
rewriting the terms associated with external ﬂow in Eqs.(4.5) and (4.9). Let αh and αt be
deﬁned as follows:

∂
(ρf V )dΩ +
ρf V (V · nj )dΓ
ˆ
Γ3 ,Γ4
Ω3 ∂t
∂
αt =
(ρf V )dΩ +
ρf V (V · nj )dΓ
ˆ
Ω4 ∂t
Γ5,6

αh =

The following replacements may now be made in Eq.(4.5):
60

∂
(ρf V )dΩ +
ρf V (V · nj )dΓ
ˆ
Γ3 ,Γ4
Ω3 ∂t
xy
yx
yy
xx
i i
k h αh · ˆ ˆ + k
αh · ˆ ˆ + k
i j
αh · ˆ ˆ + k
j i
αh · ˆ ˆ
j j
h
h
h

Replace:
With:

∂
(ρf V )dΩ +
ρf V (V · nj )dΓ
ˆ
Γ5 ,Γ6
Ω4 ∂t
yy
yx
xy
xx
j j
j i
k t αt · ˆ ˆ + k t αt · ˆ ˆ + k t αt · ˆ ˆ + k t αt · ˆ ˆ
i j
i i

Replace:
With:

where the coeﬃcients k xx and k xy scale the momentum imparted to the external ﬂuid in
the x and y directions by acceleration in the x direction. Likewise, k yx , and k yy scale
the momentum imparted to the external ﬂuid in the x and y directions by accelerations
in the y direction. The volume Ω3 is equal in size to that of the head [9] and is assumed
to surround the head with uniform thickness. The boundary of volume Ω4 is a cylinder
which circumscribes the ﬁnned tube tail [9]. As per assumption A1, both the head and tail
have three planes of symmetry and assumption A7 states that the velocity ﬁeld of the ﬂuid
surrounding the tail is adequately modeled by potential ﬂow theory. Thus the replacements
above simplify to [9]:

Replace:
With:

∂
(ρf V )dΩ +
ρf V (V · nj )dΓ
ˆ
Γ3 ,Γ4
Ω3 ∂t
yy
xx
i i
k h αh · ˆ ˆ + k
αh · ˆ ˆ
j j
h
61

(4.10)

Replace:
With:

∂
(ρf V )dΩ +
ρf V (V · nj )dΓ
ˆ
Γ5 ,Γ6
Ω4 ∂t
yy
xx
k t αt · ˆ ˆ + k t αt · ˆ ˆ
j j
i i

(4.11)

Using the substitutions in Eqs.(4.10) and (4.11), Eq.(4.5) can now be rewritten as follows

∂V
∂V
∂V
ρh
ρt
dΩ +
dΩ +
ρf
dΩ
∂t
∂t
Γ1 (t)
Ωh
Ωh ,Ωt
Ωt ∂t
Ω1 ,Ω2
∂V
∂V
xx
xx
dΩ + kt
dΩ · ˆ ˆ
i i
ρf
+
ρf V V · nj dΓ +
ˆ
ρf
kh
∂t
∂t
Ω4
Γ1 ,Γ2
Ω3
∂V
∂V
yy
yy
+
k
ρ
ρf
dΩ + kt
dΩ · ˆ ˆ
j j
h Ω f ∂t
∂t
Ω4
3
Fd dx +

+

xx
kh

+

k

Fp dΓ =

xx
ρf V V · nj dΓ + kt
ˆ

ρf V V · nj dΓ
ˆ

·ˆ ˆ
i i

yy
yy
ρf V V · nj dΓ
ˆ
ρ V V · nj dΓ + kt
ˆ
h Γ ,Γ f
Γ5 ,Γ6
3 4

·ˆ ˆ
j j

Γ3 ,Γ4

Γ5 ,Γ6

(4.12)

A replacement similar to Eqs.(4.10) and (4.11) can be performed on analogous terms in
Eq.(4.9):

Replace:
With:

∂
(R − Rc ) × ρf V dΩ +
(R − Rc ) × ρf V (V · nj )dΓ
ˆ
(4.13)
Γ3 ,Γ4
Ω3 ∂t
∂
θθ
kh
(R − Rc ) × ρf V dΩ +
(R − Rc ) × ρf V (V · nj )dΓ
ˆ
Ω3 ∂t
Γ3 ,Γ4
62

Replace:
With:

∂
(R − Rc ) × ρf V dΩ +
ρf V (V · nj )dΓ
ˆ
(4.14)
Γ5 ,Γ6
Ω4 ∂t
∂
θθ
(R − Rc ) × ρf V dΩ +
(R − Rc ) × ρf V (V · nj )dΓ
ˆ
kt
Γ5 ,Γ6
Ω4 ∂t

Using the substitutions in Eqs.(4.13) and (4.14) and replacing (R − Rc ) with (r − rc ), which
can be veriﬁed from Fig.4.1, Eq.(4.9) can now be rewritten as follows

(r − rc ) × Fd dx +
(r − rc ) × Fp dΓ =
Γ1 (t)
Ωh ,Ωt
∂
∂
(r − rc ) × ρh V dΩ +
(r − rc ) × ρt V dΩ
Ωh ∂t
Ωt ∂t
∂
(r − rc ) × ρf V V · nj dΓ
ˆ
(r − rc ) × ρf V dΩ +
+
Γ1 ,Γ2
Ω1 ,Ω2 ∂t
∂
∂
θθ
θθ
+ kh
(r − rc ) × ρf V dΩ + kt
(r − rc ) × ρf V dΩ
Ω3 ∂t
Ω4 ∂t
θθ
+ kh

Γ3 ,Γ4

(r − rc ) × ρf V

θθ
V · nj dΓ + kt
ˆ

Γ5 ,Γ6

(r − rc ) × ρf V

(4.15)

V · nj dΓ
ˆ

To account for the eﬀect of added mass, rc in the above equation is deﬁned as follows:

rc = rcx ˆ + rcy ˆ
i
j,

(4.16)

θθ
θθ
mh xh + mt xt + m1 x1 + m2 x2 + kh m3 x3 + kt m4 x4
¯
¯
¯
¯
¯
¯
rcx =
θθ
mh + mt + m1 + m2 + k θθ m3 + kt m4
h
θθ
θθ
mh yh + mt yt + m1 y1 + m2 y2 + kh m3 y3 + kt m4 y4
¯
¯
¯
¯
¯
¯
rcy =
θθ
mh + mt + m1 + m2 + k θθ m3 + kt m4
h

where yh , y1 and y3 are zero due to Assumptions A1, A2 and A6.
¯ ¯
¯
63

4.4.3

Expansion of a General Volume Term

This section contains expansion of the volume terms in Eqs.(4.12) and (4.15), which are
the force and moment equations, respectively. From Fig.4.1, the position, velocity and
acceleration of an arbitrary particle is

R = R0 + r

˙
˙
˙
˙ˆ
V = R = R0 + rxy + (θk × r)

˙
¨
¨
¨
¨ˆ
˙ˆ ˙
˙ˆ
˙ˆ
V = R = R0 + rxy + (θk × r) + (2θk × rxy ) + (θk × (θk × r))

(4.17)

where the particle’s velocity and acceleration have been obtained through diﬀerentiation.
The volume terms in Eq.(4.12) are expanded via Eq.(4.17), and are integrated separately to
give rise to physically intuitive terms:

∂V
¨
¨
¨ˆ
˙ˆ ˙
˙ˆ
˙ˆ
dΩ =
ρ R0 + rxy + (θk × r) + (2θk × rxy ) + (θk × (θk × r)) dΩ
∂t
Ω
Ω
¨
¨
˙
¨ˆ
˙ˆ
˙
ρr dΩ
=
ρ(R0 + rxy ) dΩ + θ k ×
ρr dΩ + 2θk ×
ρrxy dΩ − θ2
Ω
Ω
Ω
Ω
¨
¨
˙
¨ xj ¯i)
˙ˆ
˙ xi ¯j)
=
ρ(R0 + rxy ) dΩ + mθ(¯ˆ − yˆ + 2θk ×
ρrxy dΩ − mθ2 (¯ˆ + yˆ
Ω
Ω
ρ

(4.18)
where

Ω

ρr dΩ = m(¯ˆ + yˆ
xi ¯j)

(4.19)

The volume terms in Eq.(4.15) are likewise expanded via Eq.(4.17) and integrated to give
64

rise to the following terms




∂
∂
∂
˙
˙ 
˙
(r − rc ) × ρ (R − Rc ) dΩ +
(r − rc ) × ρV dΩ =
(r − rc ) × ρRc dΩ


Ω ∂t
Ω ∂t
Ω ∂t
˙ rc )
˙
=(r−



=

=


¨ ¨ 
¨
(r − rc ) × ρ (R − Rc ) dΩ +
(r − rc ) × ρRc dΩ +


Ω
Ω
¨ ¨
=(r−rc )
=0

Ω

Ω

˙
˙ ˙
(r − rc ) × ρRc dΩ
=0

¨
¨
¨
¨
˙ˆ
˙ˆ
(r − rc ) × ρ (rxy − rcxy ) + θ k × (r − rc ) + 2θk × (rxy − rcxy )

˙ˆ
˙ˆ
+ θk × (θk × (r − rc )) dΩ
=

Ω

¨
¨ˆ
˙ˆ
(r − rc ) × ρrxy dΩ + θk Jz − m(¯2 + y 2 ) + 2θk
x
¯

Ω

˙
(r − rc ) · ρrxy dΩ

(4.20)

In the above equation, the following simpliﬁcations were employed:

Ω

Ω

¨
(r − rc ) × ρrcxy dΩ =

˙ˆ
˙ˆ
(r − rc ) × θ k × ρ(r − rc ) dΩ = θk

Ω

Ω

Ω

¨
ρ(r − rc ) dΩ × rcxy = 0

˙ˆ
ρ [(r − rc ) · (r − rc )] dΩ = θ k Jz − m(¯2 + y 2 )
x
¯

˙
˙
˙ˆ
(r − rc ) × ρ 2θk × (rxy − rcxy ) dΩ
˙ˆ
= 2θ k

Ω

Ω

˙
˙ˆ
ρ (r − rc ) · rxy dΩ − 2θk

Ω

˙
ρ (r − rc ) · rcxy dΩ

˙ˆ
˙ˆ
(r − rc ) × ρ θ k × (θk × (r − rc )) dΩ = 0

65

4.4.4

Individual Volume Terms

In this section, the expansions in Eqs.(4.18) and (4.20) were used to compute the volume
terms in Eqs.(4.12) and (4.15). In particular, the integrals for the volumes Ωh , Ωt , and Ωi ,
i = 1, 2, 3, 4 were evaluated.

4.4.4.1

Rigid Head: Ωh

All the particles in the rigid head of the submersible are ﬁxed with respect to the xy-frame
and therefore

¨
rxy = 0,

x = xh ,
¯ ¯

y = 0,
¯

˙
rxy = 0,

Jz = Jzh

The volume terms in Eq.(4.12), expanded via Eq.(4.18), take the form

∂V
¨
¨x j ˙ ¯ i
ρh
dΩ = mh R0 + θ¯hˆ − θ2 xhˆ
∂t
Ωh

(4.21)

The volume terms in Eq.(4.15), expanded via Eq.(4.20), take the form
∂
¨ Jzh − x2
(r − rc ) × ρh V dΩ = mh θ
¯h
mh
Ωh ∂t
4.4.4.2

ˆ
k

(4.22)

Flexible Tail: Ωt

For the particles in the ﬂexible tail of the submersible, the following relations hold:

¨
rxy = yˆ
¨j,

x = xt ,
¯ ¯

y=
¯

1 L
y dx,
L 0

˙
rxy = yˆ
˙ j,
66

Jz = mt

1 L 2
L2
+
y dx
3
L 0

The volume terms in Eq.(4.12), expanded via Eq.(4.18), take the form
∂V
1 L
1 L
¨
ρt
dΩ = mt R0 +
y dxˆ + θ xt ˆ −
¨ j ¨ ¯ j
y dx ˆ
i
∂t
L 0
L 0
Ωt
L
1 L
˙1
¯ i
y dx ˆ − θ2 xtˆ +
˙ i ˙
y dx ˆ
j
−2θ
L 0
L 0

(4.23)

The volume terms in Eq.(4.15), expanded via Eq.(4.20), take the form

∂
(r − rc ) × ρt V dΩ =
Ωt ∂t



2
2
L
L
L
1
1
¨ L + 1
(x − rcx )¨ dx + θ 
y
y 2 dx − x2 −
¯t
y dx 
mt 
L 0
3
L 0
L 0
˙
+2θ

4.4.4.3

1 L
ˆ
(y − rcy )y dx k
˙
L 0

(4.24)

Head Internal Fluid: Ω1

The ﬂuid within the head ﬂows with steady velocity Uh , per assumption A4. The conduit
through which it ﬂows is ﬁxed within the head of the submersible, and therefore

¨
rxy = 0,

x = x1 ,
¯ ¯

y = 0,
¯

˙
rxy = Uhˆ
i,

L2
Jz = m1 h
3

The volume terms in Eq.(4.12), expanded via Eq.(4.18), take the form

∂V
¨
¨x j
˙ j ˙ ¯ i
dΩ = m1 R0 + θ¯1 ˆ + 2θUh ˆ − θ2 x1ˆ
ρf
∂t
Ω1
67

(4.25)

The volume terms in Eq.(4.15), expanded via Eq.(4.20), take the form
L2
∂
h − x2 dx
¨
(r − rc ) × ρf V dΩ = m1 θ
¯1
∂t
3
Ω1

0
˙ 1
ˆ
+ 2θ
(x − rcx )Uh dx k
Lh −L
h
(4.26)

4.4.4.4

Tail Internal Fluid: Ω2

For the particles in this control volume, the following relations hold:

¨
rxy = (¨ + 2U y ′ + U 2 y ′′ )ˆ
y
˙
j,
˙
rxy = U ˆ + (y + Uy ′ ) ˆ
i
˙
j,

x = x2
¯ ¯
Jz = m2

1 L
y dx,
L 0
1 L 2
L2
+
y dx
3
L 0
y=
¯

The volume terms in Eq.(4.12), expanded via Eq.(4.18), take the form
1 L
∂V
1 L
¨
ρf
dΩ = m2 R0 +
(¨ + 2U y ′ + U 2 y ′′ ) dx ˆ + θ x2 ˆ −
y
˙
j ¨ ¯ j
y dx ˆ
i
∂t
L 0
L 0
Ω2
1 L
1 L
˙
˙
+2θ U ˆ −
j
¯ i
[y + Uy ′ ] dx ˆ − θ2 x2ˆ +
˙
i
y dx ˆ
j
(4.27)
L 0
L 0

The volume terms in Eq.(4.15), expanded via Eq.(4.20), take the form
1 L
∂
(r − rc ) × ρf V dΩ = m2
(x − rcx )(¨ + 2U y ′ + U 2 y ′′ ) dx
y
˙
∂t
L 0
Ω2


2
2
L
L
1
¨ L + 1
+θ
y 2 dx − x2 −
¯2
y dx 
3
L 0
L 0
˙
+2θ

1 L
ˆ
(x − rcx )U + (y − rcy )(y + Uy ′ ) dx k (4.28)
˙
L 0
68

4.4.4.5

Head External Fluid: Ω3

The ﬂuid in this control volume ﬂows past the rigid head of the submersible with an unsteady
velocity Ue . For particles in this control volume, the following relations therefore hold

¨
˙ i,
rxy = Ue ˆ

x = x3
¯ ¯

L2
θθ
Jz = kh m3 h
3

˙
rxy = Ue ˆ
i

y = 0,
¯

The volume terms in Eq.(4.12), expanded via Eq.(4.18), take the form

xx
ρf kh
Ω3

∂V ˆ ˆ
yy
·i i+k
h
∂t

∂V ˆ ˆ
· j j dΩ =
∂t

yy ¨
xx ¨
˙ ¯ i
¨x
˙
˙
m3 kh R0x + Ue − θ2 x3 ˆ + m3 k
R0y + θ¯3 + 2θUe ˆ
j
h

(4.29)

The volume terms in Eq.(4.15), expanded via Eq.(4.20), take the form
2
∂
θθ r U + θ Lh − x2
θθ
˙e ¨
(r − rc ) × ρf V dΩ = m3 kh
¯3
kh
cy
3
Ω3 ∂t

˙
ˆ
+ 2θ(x − rcx )Ue k
(4.30)

4.4.4.6

Tail External Fluid: Ω4

The ﬂuid in this control volume ﬂows past the ﬂexible tail of the submersible with an unsteady
velocity Ue . For particles in this control volume, the following relations therefore hold

2
¨
˙ i
˙
rxy = Ue ˆ + (¨ + 2Ue y ′ + Ue y ′′ + Ue y ′ )ˆ
y
˙
j,
˙
rxy = Ue ˆ + (y + Ue y ′ ) ˆ
i
˙
j,

θθ
Jz = kt m4
69

1 L
y dx,
L 0
L2
1 L 2
+
y dx
3
L 0

x = x4 ,
¯ ¯

y=
¯

The volume terms in Eq.(4.12), expanded via Eq.(4.18), take the form

Ω4

xx
ρf kt

∂V ˆ ˆ
yy
· i i + kt
∂t

xx R + U − θ 1
¨
˙e ¨
m4 kt
0x
L

∂V ˆ ˆ
· j j dΩ =
∂t

L

L
˙1
˙ ¯ i
y dx − 2θ
(y + Ue y ′ ) dx − θ2 x4 ˆ +
˙
L 0
0

L
1 L
yy ¨
2
¨x
˙
˙ 1
˙
m4 kt R0y +
(¨ + 2Ue y ′ + Ue y ′′ + Ue y ′) dx + θ¯4 + 2θUe − θ2
y
˙
y dx ˆ
j
L 0
L 0
(4.31)

The volume terms in Eq.(4.15), expanded via Eq.(4.20), take the form

θθ
m4 kt
θθ
m4 kt

2
¨ L
θ
3
˙
+2θ

∂
(r − rc ) × ρf V dΩ =
Ω4 ∂t
1 L
2
˙
˙
(x − rcx )(¨ + 2Ue y ′ + Ue y ′′ + Ue y ′) − (y − rcy )Ue
y
˙
L 0

2
1 L 2
1 L
+
y dx − x2 −
¯4
y dx 
L 0
L 0

dx +

1 L
ˆ
(x − rcx )Ue + (y − rcy )(y + Ue y ′ ) dx k
˙
L 0

4.4.5

Individual Surface Terms

4.4.5.1

(4.32)

Jet Inlet: Γ1

The jet inlet is located at the tip of the head as shown in Fig.4.3. At this location,

r = −Lhˆ
i,

n1 = −ˆ
ˆ
i,

˙
˙ˆ
Vs = R0 + θk × −Lhˆ
i,

˙
V = (Vs + Vr ) = R0 + Uhˆ − θLh ˆ
i ˙ j,
70

i,
Vr = Uhˆ

˙
V · n1 = −R0x − Uh
ˆ

The surface terms in Eq.(4.12) have the form

Γ1

˙
˙
ρf V [V · n1 ] dΓ = ρf A1 R0 + Uhˆ − θLhˆ (−R0x − Uh )
ˆ
i ˙ j

(4.33)

The surface terms in Eq.(4.15) have the form

Γ1

(r − rc ) × ρf V [V · n1 ] dΓ =
ˆ

˙
ˆ
˙
˙
˙
ρf A1 (−Lh − rcx )R0y + rcy (R0x + Uh ) + θ(L2 + Lh rcx ) [−R0x − Uh ] k
h

4.4.5.2

(4.34)

Jet Outlet: Γ2

At the outlet of the jet,

r = Lˆ + y(L) ˆ
i
j,
Vr = U n2 ,
ˆ

˙
Vs = R0 + y(L) ˆ + θ k × Lˆ + y(L)ˆ ,
˙
j ˙ˆ
i
j

n2 ≈ ˆ + y ′ (L) ˆ
ˆ
i
j,

˙
V = (Vs + Vr ) = R0 + y(L) ˆ + U ˆ + Uy ′ (L)ˆ + θ Lˆ − y(L) ˆ ,
˙
j
i
j ˙ j
i

˙
˙
˙
˙
V · n2 = R0x + U − θy(L) + y ′ (L) R0y + y(L) + θL + Uy ′ (L)
ˆ
˙

The surface terms in Eq.(4.12) have the form

Γ2

ρf V [V · n2 ] dΓ =
ˆ

˙
ρf A2 R0 + y(L) ˆ + Uˆ + Uy ′ (L)ˆ + θ Lˆ − y(L)ˆ
˙
j
i
j ˙
j
i

×

˙
˙
˙
˙
R0x + U − θy(L) + y ′ (L) R0y + y(L) + θL + Uy ′ (L)
˙
71

(4.35)

The surface terms in Eq.(4.15) have the form

Γ2

˙
˙
(r − rc ) × ρf V [V · n2 ] dΓ = ρf A2 (L − rcx ){R0y + y(L) + θL + Uy ′ (L)}
ˆ
˙
˙
˙
˙
−{y(L) − rcy }{R0x + U − θy(L)} + θ L2 + y 2 (L) − Lrcx − y(L)rcy

×

˙
˙
ˆ
˙
˙
{R0x + U − θy(L)} + y ′ (L){R0y + y(L) + θL + Uy ′ (L)} k
˙

4.4.5.3

(4.36)

External Flow Head Inlet: Γ3

For this surface, shown in Fig.4.3,

i,
r = −Lh ˆ

n3 = −ˆ
ˆ
i,

˙
˙ˆ
i,
Vs = R0 + θk × −Lhˆ

˙
V = (Vs + Vr ) = R0 + Ue ˆ − θLh ˆ
i ˙ j,

Vr = Ue ˆ
i,

˙
V · n3 = −R0x − Ue
ˆ

The surface terms in Eq.(4.12) have the form

yy
xx
i i
ρf kh V · ˆ ˆ + k
V · ˆ ˆ [V · n3 ] dΓ =
j j
ˆ
h
Γ3
yy ˙
xx ˙
˙
˙
ρf A3 kh R0x + Ue ˆ + k
i
R0y − θLh ˆ [−R0x − Ue ]
j
h

(4.37)

The surface terms in Eq.(4.15) have the form

θθ
kh

Γ3

(r − rc ) × ρf V [V · n3 ] dΓ =
ˆ

θθ
ρf A3 kh

˙
˙
˙
−Lh − rcx R0y + rcy R0x + Ue + θ L2 + Lh rcx
h

ˆ
˙
[−R0x − Ue ] k
(4.38)

72

4.4.5.4

External Flow Head Outlet: Γ4

This surface, shown in Fig.4.3, is located at the base of the ﬂexible tail. For this surface, the
following relations hold:

r = 0,

n4 = ˆ
ˆ
i,

˙
Vs = R0 ,

˙
V = (Vs + Vr ) = R0 + Ue ˆ
i,

Vr = Ue ˆ
i,

˙
V · n4 = R0x + Ue
ˆ

The surface terms in Eq.(4.12) have the form

yy
xx
i i
ρf kh V · ˆ ˆ + k
V · ˆ ˆ [V · n4 ] dΓ =
j j
ˆ
h
Γ4
yy ˙
xx ˙
ρf A4 kh R0x + Ue ˆ + k R0y ˆ [R0x + Ue ]
i
j ˙
h

(4.39)

The surface terms in Eq.(4.15) have the form

θθ
kh

Γ4

(r − rc ) × ρf V [V · n4 ] dΓ =
ˆ

θθ
˙
˙
ρf A4 kh −rcx R0y + rcy R0x + Ue

4.4.5.5

ˆ
˙
[R0x + Ue ] k

(4.40)

External Flow Tail Inlet: Γ5

This surface also lies at the base of the ﬂexible tail. Although its normal n5 is opposite to
ˆ
that of Γ4 , the net momentum ﬂux through these surfaces is not zero, since their areas and
73

added mass coeﬃcients are not equal. For this surface,

r = 0,

n5 = −ˆ
ˆ
i,

˙
Vs = R0 ,

˙
V = (Vs + Vr ) = R0 + Ue ˆ
i,

Vr = Ue ˆ
i,

˙
V · n5 = −R0x − Ue
ˆ

The surface terms in Eq.(4.12) have the form

Γ5

xx V · ˆ ˆ + k yy V · ˆ ˆ [V · n ] dΓ =
j j
ˆ5
ρf kt
i i
t

xx R + U ˆ + k yy R ˆ [−R − U ]
˙
˙
˙
ρf A5 kt
e i
e
0x
0x
t 0y j

(4.41)

The surface terms in Eq.(4.15) have the form

θθ
kt

Γ5

(r − rc ) × ρf V [V · n5 ] dΓ =
ˆ

θθ
˙
˙
ρf A5 kt −rcx R0y + rcy R0x + Ue

4.4.5.6

ˆ
˙
[−R0x − Ue ] k

(4.42)

External Flow Tail Outlet: Γ6

This surface lies at the tip of the tail and has the same normal as that of Γ2 . For this surface,

r = Lˆ + y(L) ˆ
i
j,
Vr = Ue n6 ,
ˆ

n6 ≈ ˆ + y ′(L) ˆ
ˆ
i
j,

˙
Vs = R0 + y(L) ˆ + θk × Lˆ
˙
j ˙ˆ
i,

˙
V = (Vs + Vr ) = R0 + y(L) ˆ + Ue ˆ + Ue y ′ (L)ˆ + θ Lˆ − y(L) ˆ ,
˙
j
i
j ˙ j
i

˙
˙
˙
˙
V · n6 = R0x + Ue − θy(L) + y ′(L) R0y + y(L) + θL + Ue y ′ (L)
ˆ
˙
74

The surface terms in Eq.(4.12) have the form

Γ6

xx V · ˆ ˆ + k yy V · ˆ ˆ [V · n ] dΓ =
j j
ρf kt
i i
ˆ6
t
xx R + U − θy(L) ˆ + k yy R + θL + y(L) + U y ′ (L) ˆ ×
˙
˙
˙
˙
ρf A6 kt
j
i
e
e ˙
0y
0x
t
˙
˙
˙
˙
R0x + Ue − θy(L) + y ′(L) R0y + y(L) + θL + Ue y ′(L)
˙

(4.43)

The surface terms in Eq.(4.15) have the form

θθ
kt

Γ6

θθ
˙
˙
(r − rc ) × ρf V [V · n6 ] dΓ = ρf A6 kt (L − rcx ) R0y + y(L) + θL + Ue y ′ (L)
ˆ
˙

− y(L) − rcy

˙
˙
˙
R0x + Ue − θy(L) + θ L2 + y(L)2 − Lrcx − y(L)rcy

˙
˙
˙
˙
R0x + Ue − θy(L) + y ′(L) R0y + y(L) + θL + Ue y ′(L)
˙

4.4.6

ˆ
k

×
(4.44)

Combined Scalar Equations

The expressions in sections 4.4.4 and 4.4.5 are substituted back into Eqs.(4.12) and (4.15) to
obtain the equations of motion for translation in x and y, and rotation in θ. These equations,
which are are given subsequently, use the following terms:

xx
xx
mT x = mh + mt + m1 + m2 + kh m3 + kt m4
yy
yy
mT y = mh + mt + m1 + m2 + k m3 + kt m4
h
75

mh =
mt =
mi =

Ωh
Ωt
Ωi

1
ρ r dΩ
xh ˆ + yh ˆ =
¯ i ¯ j
mh Ω h
h
1
xt ˆ + yt ˆ =
¯ i ¯ j
ρ r dΩ
mt Ωt t
1
ρ r dΩ, i = 1, 2, 3, 4
xi ˆ + yi ˆ =
¯ i ¯ j
mi Ω f
i

ρh dΩ,

ρt dΩ,
ρf dΩ,

Translation in x:

Ωh ,Ωt

Fd dx +

Γ1 (t)

Fp dΓ

·ˆ =
i

L
xx
xx
xx 1
¨
¨
˙
mT x R0x + kh m3 + kt m4 Ue − θ (mt + m2 + m4 kt )
y dx
L 0
L
y(L)
xx 1
xx
˙
− 2θ (mt + m2 + m4 kt )
y dx + (m2 U + m4 kt Ue )
˙
L 0
L
2
xx ¯
xx ¯
˙
˙
− θ2 mh xh + mt xt + m1 x1 + m2 x2 + m3 kh x3 + m4 kt x4 − ρf A1 R0x + Uh
¯
¯
¯
¯
˙
˙
˙
˙
˙
˙
+ ρf A2 R0x + U − θy(L) R0x + U − θy(L) + y ′ (L) R0y + y(L) + θL + Uy ′ (L)
˙
2
2
xx
xx
xx ˙
˙
− ρf A3 kh R0x + Ue + ρf (kh A4 − kt A5 ) R0x + Ue
xx ˙
˙
˙
˙
˙
˙
+ρf A6 kt R0x + Ue − θy(L) R0x + Ue − θy(L) + y ′(L) R0y + y(L) + θL + Ue y ′(L)
˙
(4.45)
76

Translation in y:

Ωh ,Ωt

Fd dx +

Γ1 (t)

Fp dΓ

·ˆ=
j

yy 1
¨
mT y R0y + (mt + m2 + m4 kt )
L

L
0

y dx + m2
¨

1 L
(2U y ′ + U 2 y ′′ ) dx
˙
L 0

L

yy 1
2
˙
(2Ue y ′ + Ue y ′′ + Ue y ′ ) dx
˙
+ m4 kt
L 0
yy
yy
¨
¯
+ θ mh xh + mt xt + m1 x1 + m2 x2 + m3 k x3 + m4 kt x4
¯
¯
¯
¯
¯
h
yy
yy
˙
+ 2θ m1 Uh + m2 U + (m3 k + m4 kt )Ue
h
L
xx 1
˙ j
˙
˙
˙
y dx + ρf A1 R0y − θLh ˆ (−R0x − Uh )
− θ2 (mt + m2 + m4 kt )
L 0
˙
˙
+ ρf A2 R0y + y(L) + θL + Uy ′ (L)
˙

˙
˙
R0x + U − θy(L)

˙
˙
+y ′(L) R0y + y(L) + θL + Uy ′ (L)
˙
yy ˙
yy
yy
˙
˙
˙
˙
R0y − θLh [−R0x − Ue ] + ρf (k A4 − kt A5 )R0y R0x + Ue
h
h
yy ˙
˙
˙
˙
˙
+ ρf A6 kt R0y + θL + y(L) + Ue y ′ (L) R0x + Ue − θy(L)

+ ρf A3 k

˙
˙
+y ′(L) R0y + y(L) + θL + Ue y ′(L)
˙

77

(4.46)

Rotation in θ:

Ωh ,Ωt

(r − rc ) × Fd dx +

Γ1

(r − rc ) × Fp dΓ =

θθ
(mt + m2 + m4 kt ) L
(x − rcx )¨ dx
y
L
0

θθ L
m4 kt
m2 L
2
˙
(x − rcx )(2U y ′ + U 2 y ′′ ) dx +
˙
(x − rcx )(2Ue y ′ + Ue y ′′ + Ue y ′ ) dx
˙
L 0
L
0
L2
m k θθ L
¨ J + (m + m k θθ ) h
˙ e m k θθ rcy − 4 t
+U
(y − rcy ) dx + θ
3 h
1
3 h 3
zh
L
0

+

2
2
L
L
2 dx − (m + m + m k θθ ) 1
θθ ) L + 1
y
y dx
+ (mt + m2 + m4 kt
t
2
4 t
3
L 0
L 0
θθ ¯
θθ ¯
− mh x2 + mt x2 + m1 x2 + m2 x2 + m3 kh x2 + m4 kt x2
¯h
¯t
¯1
¯2
3
4

θθ
θθ
m1 Uh + m3 kh Ue 0
m2 U + m4 kt Ue L
˙
(x − rcx ) dx +
(x − rcx ) dx
+ 2θ 
Lh
L
−Lh
0
θθ L
θθ L
mt + m2 + m4 kt
m2 U + m4 Ue kt
+
(y − rcy )y dx +
˙
(y − rcy )y ′ dx
L
L
0
0
˙
˙
˙
˙
+ ρf A1 (−Lh − rcx )R0y + rcy (R0x + Uh ) + θ(L2 + Lh rcx ) [−R0x − Uh ]
h
˙
˙
+ ρf A2 (L − rcx ) R0y + y(L) + θL + Uy ′ (L)
˙
˙
˙
˙
−{y(L) − rcy }{R0x + U − θy(L)} + θ L2 + y 2 (L) − Lrcx − y(L)rcy

×

˙
˙
˙
˙
R0x + U − θy(L) + y ′ (L) R0y + y(L) + θL + Uy ′ (L)
˙
θθ
+ ρf A3 kh

˙
˙
˙
−Lh − rcx R0y + rcy R0x + Ue + θ L2 + Lh rcx
h

θθ
θθ
˙
˙
+ ρf (A4 kh − A5 kt ) −rcx R0y + rcy R0x + Ue

˙
[−R0x − Ue ]

˙
[R0x + Ue ]

θθ
˙
˙
+ ρf A6 kt (L − rcx ) R0y + y(L) + θL + Ue y ′ (L)
˙
− y(L) − rcy

˙
˙
˙
R0x + Ue − θy(L) + θ L2 + y(L)2 − Lrcx − y(L)rcy

˙
˙
˙
˙
R0x + Ue − θy(L) + y ′ (L) R0y + y(L) + θL + Ue y ′ (L)
˙

78

ˆ
k

×
(4.47)

4.4.7

Pressure and Drag Forces

In this section, expressions for the pressure and drag forces, Fp and Fd that appear on the
left side of Eqs.(4.45), (4.46) and (4.47) are provided. The force exerted by the pressure
acting at the inlet is computed via analysis of a captured quasi-steady streamtube:

Γ1

Fp dA = −A1 ρf Uh (Uh − Ue ) ˆ
i

(4.48)

On the left side of the scalar equations in the previous section, the pressure Fp is integrated
over the surface of the inlet Γ1 . It is assumed to be uniform over this area.
An expression for the drag force is obtained by ﬁrst writing the components of the velocity
T
N
of the rigid head and ﬂexible tail in the normal and tangential directions (Vh , Vh , VtN ,
VtT ):
N
˙
˙
˙
˙
˙
VtN = R0y + y + θx − (R0x − θy)y ′
˙ ˙
Vh = R0y + θx
(4.49)
T
˙
˙
˙
˙
VtT = R0x − θy + (R0y + y + θx)y ′
˙ ˙
Vh = R0x
Using a model similar to that proposed by [42]1 , the distributed drag forces on the rigid head
N
T
N
T
and ﬂexible tail in the normal and tangential directions (Dh , Dh , Dt , Dt ) are written as
N
N T N
N
Dh = −ρf Ph (Ch Vh Vh + ch Vh )

N
N
Dt = −ρf h(Ct VtT VtN + ct VtN )

T
T T T
Dh = −ρf Ph Ch Vh |Vh |

T
T
Dt = −ρf Pt Ct VtT |VtT |

(4.50)

N
where Ph , Pt are the perimeters of the head and tail in the yz plane. The coeﬃcients Ch ,
T
T
N
Ct , Ch and Ct are empirical drag coeﬃcients in the normal and tangential directions;
N
1: While terms similar to ch Vh or ct VtN do not appear in [42], later workers recognized
the need for a linear correction. Some of these works can be found in [31].
79

and ch and ct are linear damping coeﬃcients.
The above expressions for drag are used to derive the components of drag in the x and
x y
x y
y directions (Dh , D , Dt , Dt ):
h
x
T
D h = Dh

x
T
N
Dt = Dt − Dt y ′

y
N
D = Dh
h

y
N
T
D t = Dt + Dt y ′

(4.51)

Substitution of Eqs.(4.49) and (4.50) into (4.51) yields

x
T ˙
˙
Dh = −ρf Ph Ch R0x R0x
y
N ˙
˙
˙
˙
˙
D = −ρf Ph Ch R0x R0y + θx + ch R0y + θx
h

(4.52)

and

T
x
Dt = − ρf Pt Ct

˙
˙
˙
R0x − θy + R0y + y + θx y ′
˙ ˙

N ˙
˙
˙
+ ρf h Ct R0x − θy + R0y + y + θx y ′
˙ ˙
˙
˙
˙
˙ ˙
+ct R0y + y + θx − (R0x − θy)y ′

˙
˙
˙
R0x − θy + R0y + y + θx y ′
˙ ˙
˙
˙
˙
R0y + y + θx − (R0x − θy)y ′
˙ ˙

y′

y
N ˙
˙
˙
Dt = − ρf h Ct R0x − θy + R0y + y + θx y ′
˙ ˙

˙
˙
˙
R0y + y + θx − (R0x − θy)y ′
˙ ˙

˙
˙
˙
+ct R0y + y + θx − (R0x − θy)y ′
˙ ˙
T
− ρf Pt Ct

˙
˙
˙
R0x − θy + R0y + y + θx y ′
˙ ˙

˙
˙
˙
R0x − θy + R0y + y + θx y ′
˙ ˙

y′

(4.53)

In addition to these distributed viscous forces, the low pressure region at the blunt trailing
end of a slender body is an additional source of drag. Consistent with the literature [31],
80

this communication refers to this as a drag force, rather than a pressure force. These forces,
commonly referred to as base drag, are given by the following expressions:

b
Dh = −

ρf
A Cb V T V T
2 h h h h

,
x=0

b
Dt = −

ρf
A C b V T VtT
2 t t t

(4.54)
x=L

b b
where Ah , At are the cross-sectional areas of the head and tail in the yz plane, and Ch , Ct
are the base drag coeﬃcients of the head and tail. Using Eqs.(4.52), (4.53) and (4.54), the
drag force Fd in Eqs.(4.45), (4.46), and (4.47) is written in component form as follows

Ωh ,Ωt
Ωh ,Ωt

4.5

4.5.1

Fd dx · ˆ =
i
Fd dx · ˆ =
j

Ωh

x
Dh dΩ +

Ωt

x
b
b
Dt dΩ + Dh + Dt

y
y
b
Dt dΩ + Dt y ′(L)
D dΩ +
h
Ωh
Ωt

Dynamics of the Flexible Tail

Equations of Motion

In section 4.4 a generic function y(x, t) was used to describe the motion of the ﬂexible tail.
In this section, we will derive the equation of motion of the tail such that y(x, t) can be
computed. To this end, Fig.4.4 presents the free-body diagrams of diﬀerential elements
within control volumes Ω2 , Ωt and Ω4 . A balance of forces for each element in the x and y
directions is provided below:
81

q S dx

1 ∂p
A2 (p2 + 2 2 dx)
∂x

FI dx
wall of fluid-conveying tube

1 ∂p
A2 (p2 − 2 2 dx)
∂x

1 ∂p
A6 (p4 + 2 4 dx)
∂x

FE dx

q S dx
∂p
A6 (p4 − 1 4 dx)
2 ∂x

dx

(a)

(c)
y
Dt

x
Dt

y

∂Q
Q+ 1
2 ∂x dx
q S dx
1
M + 2 ∂M dx
FE dx
∂x
1
M − 2 ∂M dx
∂x
1
FI dx q S dx
T + 2 ∂T dx
1 ∂Q dx
∂x
Q− 2
∂x
x
1
T − 2 ∂T dx
∂x

(b)
Figure 4.4: Free-body diagrams of (a) a diﬀerential element of Ω2 (tail internal ﬂuid), (b)
a diﬀerential element of Ωt (ﬂexible tail), and (c) a diﬀerential element of Ω4 (tail external
ﬂuid).

82

Tail Internal Fluid:

∂p
∂y ˆ
∂
−A2 2 − qS − FI
i + FI − A2
∂x
∂x
∂x

∂p2 ∂y
∂x ∂x

− qS

∂y ˆ m2 dV
j=
∂x
L dt

(4.55)

where

∂V
¨
˙ j
= R0 + y + 2U y ′ + U 2 y ′′ ˆ + θ x ˆ − y ˆ + 2θ U ˆ − (y + Uy ′ ) ˆ
¨
˙
j ¨ j
i
˙
i
∂t

Flexible Tail:
∂2y
∂y
∂T
∂y
x i+
+ qS + FI
−Q
+ Dt ˆ
− FE
2
∂x
∂x
∂x
∂x
∂Q
∂
−FI +
+
∂x
∂x

∂T ∂y
∂x ∂x

+ qS

∂y
m ∂V
y
j
+ FE + Dt ˆ = t
∂x
L ∂t

(4.56)

where

∂V
¨
˙ ˙i ˙
= R0 + y ˆ + θ x ˆ − y ˆ − 2θyˆ − θ2 x ˆ + y ˆ
¨j ¨ j
i
i
j
∂t

Tail External Fluid:

∂p
∂
∂y ˆ
i − FE + A6
−A6 4 + FE
∂x
∂x
∂x
m4
L

xx
kt

∂V ˆ ˆ
·i i+
∂t
83

∂p4 ∂y
∂x ∂x

yy ∂ V ˆ ˆ
·j j
kt
∂t

ˆ=
j
(4.57)

where

∂V
¨
2
˙ j
˙
˙ i ¨ j
= R0 + y + 2Ue y ′ + Ue y ′′ + Ue y ′ ˆ + Ueˆ + θ x ˆ − y ˆ + 2θ U ˆ − (y + Uy ′ )ˆ
¨
˙
j
i
˙
i
∂t
˙
i
j
− θ2 x ˆ + y ˆ

By combining the x-components of Eqs.(4.55), (4.56) and (4.57), we obtain the following
expression which will be used to simplify the expression that will be obtained by summation
of the y components:

∂p
∂p
∂T
− A2 2 − A6 4 =
∂x
∂x
∂x

xx
mt + m2 + m4 kt
L

xx
m4 kt
x
˙
ax +
Ue − Dt
L

(4.58)

where

¨
˙ ˙
ax = R0 · ˆ − θy − 2θ y +
i ¨

xx
m4 kt
m2
∂y
∂y
˙2
+
U
xx ∂x m + m + m k xx Ue ∂x − θ x
mt + m2 + m4 kt
t
2
4 t
(4.59)

Summation of the y-components of Eqs.(4.55), (4.56) and (4.57) gives:
yy
2
m2 2 m4 kt
∂
∂y
∂4y
y
2 ∂ y
+
(T − A2 p2 − A6 p4 )
=
U +
Ue
Dt − EI
∂x
L
L
∂x4 ∂x
∂x2
yy
yy
yy
m4 kt
m4 kt
m2
∂y mt + m2 + m4 kt
∂2y
˙
+2
Ue
U+
Ue
+
+
L
L
∂x∂t
L
∂x
L

∂2y
+ ay
∂t2
(4.60)

84

where

¨
˙
ay = R0 · ˆ + θx + 2θ
j ¨

yy
m4 kt
m2
˙2
yy U +
yy Ue − θ y
mt + m2 + m4 kt
mt + m2 + m4 kt

(4.61)

The combined tension-pressure gradient term in Eq.(4.60) can be rewritten as:

∂p
∂p
∂
∂y
∂T
(T − A2 p2 − A6 p4 )
=
− A2 2 − A6 4
∂x
∂x
∂x
∂x
∂x

∂y
∂2y
+ (T − A2 p2 − A6 p4 )
∂x
∂x2

Using Eq.(4.58), terms on the right hand side of the above equation can be written as
xx
xx
mt + m2 + m4 kt
m4 kt
x
˙
Ue − Dt
ax +
L
L
L ∂
(T − A2 p2 − A6 p4 )dx + T (L) − A2 p(L) − A6 p4 (L)
(T − A2 p2 − A6 p4 ) = −
x ∂x
b
=Dt

∂p
∂p
∂T
− A2 2 − A6 4
∂x
∂x
∂x

=

Note that the ﬁnal term in the above equation represents the low pressure region at the
blunt trailing edge of the tail, which is commonly termed the base drag. Substitution of the
above expressions into Eq.(4.60), yields the equation of motion for transverse vibration of
the ﬂexible tail:
yy
m2 2 m4 kt
∂4y
2
U +
Ue
EI
+
L
L
∂x4
xx
2
L mt + m2 + m4 k xx
t a + m4 kt U − D x dx − D b ∂ y
˙e
+
x
t
t ∂x2
L
L
x
yy
xx
xx
m4 kt − m4 kt
mt + m2 + m4 kt
x ∂y
˙
+ −
Ue + Dt
ax −
L
L
∂x
yy
yy
m4 kt
mt + m2 + m4 kt
m2
∂2y
∂2y
y
+ ay − Dt = 0
+2
U+
Ue
+
L
L
∂x∂t
L
∂t2
85

(4.62)

The boundary conditions for the ﬂexible tail are the familiar expressions for a cantilever:

∂y(0, t)
=0
∂x

y(0, t) = 0

∂ 2 y(L, t)
=0
∂x2

∂ 3 y(L, t)
=0
∂x3

(4.63)

Note that while the system presented here is similar to that derived in [19], the boundary
condition at x = 0 is not the same; in that work, the authors used an inertial frame at x = 0,
and allowed the rigid head to move in the y-direction. In contrast, since the boundary
conditions for the present system are written in a non-inertial frame which is ﬁxed to the
undeformed neutral axis of the tail, the simple equations for a cantilever are used.

4.5.2

Discretization

A ﬁnite-diﬀerence scheme was used to discretize the PDE in Eq.(4.62), in which the continuous function y(x, t) is considered to have discrete values y1 , y2 , · · · , yN at the locations
x1 , x2 , · · · , xN . The ﬁve-point stencils in Eq.(4.64) were used to approximate the required
spatial derivatives at the i-th point:

∂yi
∂x
2y
∂ i
∂x2
∂ 3 yi
∂x3
∂ 4 yi
∂x4

1
y
− 8yi−1 + 8yi+1 − yi+2 ,
12∆ i−2
1
−yi−2 + 16yi−1 − 30yi + 16yi+1 − yi+2 ,
=
12∆2
=

1
−yi−2 + 2yi−1 − 2yi+1 + yi+2 ,
2∆3
1
yi−2 − 4yi−1 + 6yi − 4yi+1 + yi+2 ,
=
∆4
=

86

(4.64a)
(4.64b)
(4.64c)
(4.64d)

where ∆ is the spacing between the points xi+1 and xi . The central diﬀerence relationships
provided in Eq.(4.64) were used for all points on the ﬂexible tail; “mirror” conditions were
used to handle the edge values and implement the boundary conditions. In this method,
additional points are created at x < 0 and x > L such that evaluating the stencils in
Eq.(4.64) yields the boundary conditions of Eq.(4.63). The boundary condition ∂y/∂x(0) = 0
is implemented by creating the points y0 and y−1 :

y0 = y2

y−1 = y3

The condition y(0) = 0 is implemented by simply setting y1 = 0. The two conditions at
x = L are created by rearranging the expressions for ∂ 2 y/∂x2 = 0 and ∂ 3 y/∂x3 = 0 yields:

−y(N−2) + 16y(N−1) − 30yN = −16y(N+1) + y(N+2)
−y(N−2) + 2y(N−1) = 2y(N+1) − y(N+2)

which can be solved for the values y(N+1) and y(N+2) :

y(N+1) =

15yN − 9y(N−1) + y(N−2)
7

,

y(N+2) =

−30yN + 32y(N−1) + 9y(N−2)
7

The scalar equations, Eqs.(4.45), (4.46) and (4.47) also require integration of y(x, t); to
approximate these integrals, a discrete trapezoidal approximation is used, such that
L
0



y dx ≈ y1 + yN +

1
2

N−1
i=2



yi  ∆

(4.65)

After applying these discrete approximations to the integrodiﬀerential terms, N + 3 ordinary
87

diﬀerential equations corresponding to the three scalar equations remain: Eqs.(4.45), (4.46)
and (4.47) and the N equations that result from discretization of Eq.(4.62). These N + 3
equations are second-order in time, and are transformed into a system of 2N + 6 ﬁrst order
equations. These equations are solved with a fourth-order Runge-Kutta algorithm.

4.6
4.6.1

Numerical Simulations
Simulation Parameters and Procedure

Simulations were conducted using three diﬀerent conﬁgurations for the tail: a ﬂexible tail
with a ﬁnned-tube geometry as shown in Fig.4.5, a rigid tail with identical geometry to that
of the ﬂexible tail, and a rigid tube which does not have the vertical ﬁns which exist on the
other two conﬁgurations. The two conﬁgurations of the submersible with rigid tail structures
were simulated to provide baselines measures of performance which were compared to the
performance of the submersible with the ﬂexible ﬂuttering tail. Experimental studies of the
rigid conﬁgurations were used to obtain good estimates of the viscous drag coeﬃcients used in
the simulations. The values of the parameters used for the ﬂexible tube conﬁguration in these
simulations are provided in Table 4.2. All values are the same for the rigid tail conﬁguration,
except that the Young’s modulus E is irrelevant since the tail does not oscillate.

4.6.2

Estimation of Drag Coeﬃcients

b
b
N
N
T
T
The second-order drag coeﬃcients Ch , Ct , Ch , Ct , Ch , and Ct , and the linear coeﬃcients ch and ct have been estimated using a combination of analysis of experimental results
and extrapolation of published correlations. The tangential and base drag coeﬃcients are
88

z
ﬁn

tube
y

y(x, t)

x

x

h

Figure 4.5: A ﬂuid-conveying beam in the form of a ﬁnned-tube with ﬁn height h.

calculated by analyzing the speeds of the rigid tube and rigid tail conﬁgurations. For these
conﬁgurations, the thrust is well-known, since it comes only from the jet. The rigid state
of the system means that an expression for the balance of forces on the submersible in the
x-direction can be easily expressed as:


A1 ρf Uh (Uh − Ue ) = ρf 

b
Ch Ah
2


bA
C t
T
T
2
+ t
+ Ch Ph Lh + Ct Pt L Ue
2

(4.66)

Equation (4.66) can be written for both the case of the rigid tail and the rigid tube using values of external velocity Ue which are determined experimentally. Two additional equations
89

Table 4.2: Parameters used in simulations
Masses (kg)
mh
4.200

mt

m1

m2

m3

m4

0.045

0.020

0.017

1.800

1.600

h
kθθ

t
kθθ

1.00

1.00

h
kxx

t
kxx

0.10

Added mass coeﬃcients
h
t
kyy
kyy

0.10

1.00

1.00

Center of mass locations (m)
xh
¯
−0.20

xt
¯

x1
¯

x2
¯

x3
¯

x4
¯

0.17

−0.20

0.17

−0.20

0.17

Lengths (m)
Lh
0.390

L

Ph

Pt

h

0.340

0.239

0.163

0.067, 0.072, 0.077

Miscelleneous (SI units)
Uh
4.0

Ah
45.60

U

ρf
1000

4.0

Jzh

0.217
Areas ×10−4 (m2 )

I

E

5.65 × 10−10

1.30 × 106

At

A1

A2

A3

A4

A5

A6

0.495

0.500

0.500

45.60

45.60

46.80

46.80

b
Ct

ch (m/s)

ct (m/s)

0.0790

0.116

0.116

N
Ch

N
Ct

T
Ch

0.0118

0.0039

0.0060

Drag coeﬃcients
T
b
Ct
Ch
0.0020
0.0830

from [21]1 relate the tangential and base drag coeﬃcients:

b
Ch =

0.029
C T Ph Lh /Ah
h

b
Ct =

1: Chapter 3 in that work.
90

0.135
3 (C T S + C T P L)/A
t
t t
h h

(4.67)

In the above equations, the denominator is a measure of the tangential drag on the “forebody” of the body in question, i.e. the portion of the submersible that is ahead of the base.
b
b
The diﬀerent forms of the expressions for Ch and Ct reﬂect the ﬂat aspect ratio of the tail’s
cross-section relative to that of the head; the two-dimensional form of the correlation found
in [21] has been used for the tail. Equations (4.66) and (4.67) may now be solved for the
b
T
T
b
tangential and base coeﬃcients Ch and Ct , and Ch and Ct . The second-order normal
N
T
T
N
coeﬃcients have been set as Ch = 2.0 Ch and Ct = 2.0 Ct Pt /h, per the discussion in
N
[31]. The additional Pt /h scaling on Ct is to account for the diﬀerent length scaling for
the tail in the tangential and normal directions.
The linear coeﬃcients ch and ct can be estimated by using a procedure similar to that
employed by [6]. In that work, Batchelor analytically derived an estimate for the viscous
drag produced by the boundary layer on an oscillating cylinder; this result provides an
estimate for ch . Batchelor’s method [6] can be extended to produce an estimate for the drag
coeﬃcient ct on an oscillating ﬂat plate, that approximates the tail. This procedure yields
estimates of ct = 0.0065 m/s and ch = 0.0108 m/s. It should be mentioned that Batchelor’s
procedure estimates the velocity and velocity gradient at the surface using a method more
appropriate to ﬂow with very low Reynolds number, and the coeﬃcients cited in Table 4.2
have been found to produce oscillations more similar in magnitude to those of the physical
device.

4.6.3

Simulation Results

Using the methods and parameters given above, the equations of motion have been solved
for the ﬂexible tail, rigid tail and rigid tube cases. The forward speed of the submersible
91

using each of these conﬁgurations is shown in Fig.4.6, along with the deﬂection of the tip
of the tail (x = L) and orientation of the rigid head for the ﬂexible tail conﬁguration. The
maximum speed of the ﬂexible tail conﬁguration exceeds that of the rigid tail conﬁguration,
indicating that the ﬂuttering motion of the tail generates a net thrust. A large increase in
the amplitude of oscillation of the tail tip is observed near t = 10 sec - this coincides with
an increase in the speed of the submersible with the ﬂexible tail conﬁguration.

0.023
−0.023

rigid tail

21.20

rigid tube
ﬂexible tail

0.75

23.03

0.05
y(L)

Ue (m/s)

0.85

0.00
−0.05

0.65

0.05
θ

0.55

0.00
−0.05

0

10

20

30

time (sec)

0

10

20

30

time (sec)

Figure 4.6: The subplot on the left shows the forward speed of the submersible with ﬂexible
tail, rigid tail, and rigid tube conﬁgurations. On the right, the displacement of the tail tip
y(L) and orientation θ of the rigid head of the submersible are given for the ﬂexible tail
conﬁguration.

The steady-state speed of the ﬂexible tail conﬁguration is somewhat faster than the
dimensionally identical rigid tail and somewhat slower than the ﬁnless rigid tube; this trend
is consistent with experimental trials conducted on these conﬁgurations presented in section
5.4. The lower speed of the ﬂexible and rigid tail conﬁgurations is the result of the drag
created by the ﬁns aﬃxed to the tube. The additional thrust created by the ﬂuttering action
92

of the ﬂexible tail is not suﬃcient to overcome this drag for the geometry simulated in Fig.4.6.
A ﬂuttering tail conﬁguration can however exceed the speed of a rigid tube - this is shown
later in this section.
As an aside, note that the oscillations of the tail tip shown in Fig.4.6 become bounded as
time increases. This would not be expected in an unstable linear system, in which perturbations grow unbounded as t → ∞. The system simulated is indeed non-linear, unlike the
systems presented in Chapters 2 and 3. The non-linear terms responsible for the bounded
behavior of the tail are the quadratic drag terms presented in section 4.4.7, e.g., VtT VtT .
Since these terms are quadratic with respect to the speed of the tail, increasing the tail’s
amplitude causes a large increase in drag, eﬀectively bounding the tail’s motion.
The traveling wave nature of the ﬂexible tail’s motion is made clear in Fig.4.7, which
provides a set of images of the submersible over one cycle of oscillation. The dynamic
coupling between the tail’s motion and the rigid head is evident from these images where it
is observed that the angle of the rigid head is not constant throughout the cycle. The period
of oscillation for the tail of 0.3 sec is very similar to the experimentally observed value of
≈ 0.33 sec.
Additional simulations have been performed using diﬀerent tail geometries. To conform
to assumption A1, all simulated tails have rectangular planforms. Figure 4.8 presents the
forward speeds of three ﬂexible tails with h = 0.077 and lengths L = 0.30, 0.34 and 0.38.
The speed achieved with the L = 0.30 tail is comparable to that achieved with the rigid tail
shown in Fig.4.6. This is because the L = 0.30 tail is on the border of ﬂutter instability this can be veriﬁed from the amplitude of tail oscillation shown in Fig.4.8. The L = 0.34
tail exhibits the highest forward speed of these conﬁgurations, although the behavior of the
93

y

∆x = 0.265

x

Figure 4.7: Sequence of images of the submersible with the ﬂexible tail. The scaling for the
y axis has been adjusted to make the tail oscillations more visible.
L = 0.38 tail also indicates thrust production by ﬂuttering action. The oscillations in the
forward speed of the L = 0.38 tail are the result of that tail temporarily stabilizing. An
analytical treatment of a simpliﬁed version of the submersible showed similar stabilization
behavior (Hellum et al., 2010).
In Fig.4.9, the forward speeds of three L = 0.34 tails are presented; the heights of the
tails are diﬀerent and equal to h = 0.067, 0.072, and 0.077. A tail height of h = 0.082 was
also simulated. The data for this conﬁguration is not shown in Fig.4.9 since it did not exhibit
ﬂutter instability; its forward speed is presented in Table 4.3. The h = 0.072 tail exhibits
the highest forward speed of these conﬁgurations, although the behavior of both the 0.067
and 0.077 tails also indicate thrust production by ﬂuttering action. Note that the steady
state speed of the h = 0.072 tail is Ue = 0.933, meaning that this geometry is faster than
the rigid tube conﬁguration.
A summary of the data in this section is provided in Table 4.3. It is interesting that
the two ﬂexible tail conﬁgurations with minimal oscillations, (h, L) = (0.077, 0.30) and
94

Table 4.3: Diﬀerent tail geometries and their forward speeds and
amplitudes.
Illustration Conﬁguration
L
h
y(L)
Flexible Tail
0.34 0.077 0.023
Fig.4.6
Rigid Tail
0.34 0.077
Rigid Tube
0.34
Flexible Tail
0.30 0.077 0.001
Fig.4.8
Flexible Tail
0.34 0.077 0.023
Flexible Tail
0.38 0.077 0.020
Flexible Tail
0.34 0.067 0.025
Fig.4.9
Flexible Tail
0.34 0.072 0.044
Flexible Tail
0.34 0.077 0.023
Flexible Tail
0.34 0.082
0

steady-state oscillation
Ue
0.889
0.851
0.894
0.857
0.889
0.876
0.874
0.933
0.889
0.849

(h, L) = (0.082, 0.34) have forward speeds which are close to, but not equal to, that of the
rigid tail conﬁguration with (h, L) = (0.077, 0.34). This is because of the change in surface
area and associated drag of these two conﬁgurations relative to the rigid tail.

y(L)

0.95

L = 0.30
L = 0.34
L = 0.38

0.75

y(L)

Ue (m/s)

0.85

0.005
−0.005

L = 0.30

0.05
−0.05

L = 0.34

y(L)

0.65
0.55
0

10

20

30

time (sec)

0.05
−0.05

L = 0.38

0

10

20

30

time (sec)

Figure 4.8: The subplot on the left shows the forward speed of the submersible with L = 0.30,
0.34 and 0.38 where h = 0.077. On the right, the displacement of the tail tip y(L) is shown
for these tail conﬁgurations.
95

h = 0.072

h = 0.077
h = 0.067

0.85

0.005
−0.005

h = 0.067

0.05

0.75

y(L)

Ue (m/s)

y(L)

0.95

−0.05
h = 0.072

y(L)

0.65
0.55
0

10

20

30

time (sec)

0.05
−0.05

h = 0.077

0

10

20

30

time (sec)

Figure 4.9: The subplot on the left shows the forward speed of the submersible with h =
0.067, 0.072 and 0.077 where L = 0.34. On the right, the deﬂection of the tail tip y(L) is
shown for these tail conﬁgurations.

4.7

Validation

The most robust form of simulation validation is comparison to experimental data, which
will be presented in the following section. While the complexity of the equations presented
in this chapter precludes comparison to an analytical solution, it is possible to compare
a numerical solution of a simpliﬁed version of the system to a analytical solution of that
simpliﬁed system. This approach cannot provide validation of system’s full dynamics, but
can at least provide assurance that the numerical method and base terms are correct. Here, a
comparison between the numerical and analytical solutions of a cantilever will be undertaken.
The analytical method used will be the same as that presented in Chapter 2. The tail
simulated will be the same as the “baseline” case from the conﬁguration study in the prior
section (h = 0.772 m, L = 0.34 m). Consistent with the modeling from Chapter 2, the
96

4
3

y(L)

[mm]

2
1
0
−1
−2
−3
−4
0

2

4
time

6
[s]

8

10

Figure 4.10: Displacement of the tail tip y(L) in mm for a cantilever, U = Ucr = 2.372
m/s. Neutral stability is evident, since the oscillations of the tail are neither growing nor
shrinking over time.
eﬀects of the external ﬂow will be neglected. For this system, the non-dimensional mass
ratio, critical velocity, and critical frequency are

m2
= 0.274
mt + m2
m2 1/2
LUcr = 6.619
ucr =
LEI
mt + m2 1/2
ωcr =
Ωcr L2 = 13.801
LEI
β=

where Ucr is the internal velocity at the onset of instability, and Ωcr is the radian frequency
of oscillation at the onset of instability. For this system, the equations above predict Ucr =
2.372 m/s. The full equations of motion were constrained such that the rigid head was held
stationary, and solved with U = Ucr .
97

Figure 4.10 shows the displacement of the tail tip y(L) of the simulated tail. Notice that
the amplitude is constant over the length of the trial; since the system being simulated has
no energy dissipation (drag), a constant amplitude of oscillation means that the simulation
predicts a neutrally-stable system for u = 6.619. The period of oscillation is 0.827 sec., which
yields a critical frequency value of ωcr = 13.79, which agrees with analytical predictions to
within 0.1 percent. This agreement, coupled with the experimental comparisons presented
in the next chapter, provides some conﬁdence in the simulation method presented in this
chapter.

98

Chapter 5
Experimental Studies and Control

5.1

Introduction

Devices which moved like ﬁsh have been produced by the academic community since the
1930’s, with a series of results by James Gray. Gray, using mechanical [13] and electric/galvanic means [14], was able to induce swimming motions in the bodies of dead eels.
A great deal of more recent work has been done in the Draper Laboratories at MIT, which
produced the well-known RoboTuna [5, 43]. RoboTuna and its sister RoboPike use articulated tails covered by ﬂexible sheaths. The tails of these vessels are individually actuated
such that the resulting waveforms closely resemble those of living ﬁsh, a technique known as
“biomimicry”.
Articulated systems comprise the bulk of the experimental platforms; a non-biomimetic
articulated system was built by McMasters and coworkers [27] which used a novel geared
mechanism to produce a ﬁsh-like motion using only a single servomotor. Saimek and Li
[35] designed a two-link tail with a motor at the base of the ﬁrst link to achieve heave and
99

pitch motion similar to that described in Anderson, et al. [2]. Yu, et al. [45] constructed
an articulated system and designed tail motion for three types of turning behaviors. The
work on articulated systems has not been conﬁned to “traditional” servomotor transmissions;
Guo, et al. [16] and Chen, et al. [10] used electrically active polymers for actuation.

An alternate approach was taken by Alvarado and Youcef-Toumi [1], who described a
system which used a ﬂexible tail excited with a single actuator at the base. Paidoussis’ work
[28] on a system similar to the one proposed in this communication also falls into this class
of “underactuated” tails. An intermediate form between the inﬁnitely continuous beam-like
tails and discrete systems was described by Harper et al.[18], who constructed an articulated
tail which was actuated in series with lateral and torsional springs. A properly designed
underactuated tail has natural dynamics which are similar to the desired tail waveform,
which allows the use of a lower bandwidth apparatus. Underactuated systems are also likely
to have lower transmission losses, since there are unactuated degrees of freedom. Helpful
natural dynamics and unactuated degrees of freedom do lead to certain other diﬃculties,
however. A ﬂexible tail which naturally vibrates with a good waveform for straight-line
swimming will need to resist those dynamics in order to maneuver. Control of unactuated
degrees of freedom is likewise more diﬃcult, since the system’s natural dynamics must be
taken into consideration.

This chapter is organized as follows. Section 5.2 describes an experimental platform
designed at Michigan State University, and provides some preliminary results from that
device. Section 5.3 extends the analysis provided in the previous chapter to include timevariable internal ﬂow. A method which can be used to turn the submersible is described.
100

5.2

Experiments

A series of tests has been performed in the diving well of Michigan State University’s intramural pool with an experimental realization of the submersible described in the previous
chapter. This prototype is shown in Fig.5.1. The submersible’s speed was determined for
the ﬂexible tail, rigid tail and rigid tube conﬁgurations which are described in the previous
section. The presence of the tube, as opposed to simply removing all attachments from the
tail barb, was required in order to stabilize the vehicle; it was found that removing the tube
led to a pitch instability of the hull. Note that the simulations of the rigid tube conﬁguration
described in the previous section account for the drag on this tube, as opposed to simply
simulating the “no tail” case, as might seem natural.

antenna for wireless
communication
intake nozzle

fluid-conveying tube

tail fin

recharging port

0.82 m
Figure 5.1: Exterior of the current prototype SPI, with marked features of interest. For
interpretation of the references to color in this and all other ﬁgures, the reader is referred to
the electronic version of this dissertation.

The prototype’s ellipsoidal Delrin hull is sealed via an integral silicone gasket; the hull
is of monocoque construction, in that there is no separate “frame” upon which internal
101

components are mounted. The intake nozzle is also made of Delrin, and has an elliptical
proﬁle to reduce losses. The ﬂuid is expelled from the hull through a hose barb upon which
the tails are ﬁtted. The prime mover is a Rule inline bilge pump rated at 500 gallons per
hour. We have found that the major and minor losses of the ﬂow system limit the ﬂow
rate to approximately 40% of this rated value. This prime mover and the control electronics
are powered by 20 NiMH AA cells. The control electronics consist of an ATmega1284
microcontroller which provides a PWM signal to the prime mover’s driving circuitry and
receives wireless communications from the surface; this communication method requires
that testing occur at or near the surface. The hull was designed to be excessively buoyant
with these components; to achieve nearly neutral buoyancy, sticks of lead are added to the
interior. This allows us to move the center of mass of the hull away from the centerline to
counteract roll moments. The center of mass can also be moved fore and aft, in order to
change the pitch of the submersible.

The submersible’s ﬂexible tail (L = 0.34 m) is made of two 1/20 inch thick latex sheets
bonded to a latex tube (5/16 inch ID, 3/8 inch OD) using waterproof silicone RTV. The
speciﬁc gravity of the latex is ≈ 0.96, and is approximately neutrally buoyant. The tail’s
geometry is that of a ﬁnned tube, as shown in Fig.4.5, and the height of the tail is h = 0.077
m. The rigid tail, which has identical span and length, is made of an identical latex tube
as that of the ﬂexible case, bonded to two sheets of polypropylene with waterproof silicone
RTV. The polypropylene has a speciﬁc gravity of ≈ 0.95 and is nearly neutrally buoyant.
The rigid tube is made of a polycarbonate tube (3/8 inch ID, 7/16 inch OD) which has been
found to ﬁt over the barbed ﬁtting at the rear of the submersible without leaks. A small
aluminum nozzle (5/16 inch ID) has been ﬁtted into the end of the rigid tube to ensure that
102

the exit velocity of this conﬁguration is the same as the ﬂexible and rigid tail cases.
The SPI prototype described here was designed primarily by Mr. Paul Streﬂing, and his
M.S thesis [40] contains additional details regarding its construction.

5.2.1

SPI Performance

Figure 5.2: High-speed images of the submersible with the ﬂexible tail acquired over a period
of 1 sec, soon after it started from rest in the MSU swimming pool. The average speed of the
submersible in these images was ≈ 0.4 BL/s, less than the maximum speed of ≈ 1.15 BL/s.
Note the large deﬂections of the tail which occur during acceleration.

Figure 5.2 is sequence of images taken over 1 sec of operation; these images were acquired
while the submersible was accelerating from rest, during which time pictures could be taken
more normal to the free surface to reduce glare. The reduced forward speed during acceleration also leads to a longer period of oscillation of the tail; Fig.5.2 indicates a period of 0.8
seconds, compared to a period of ≈ 0.33 sec. observed when the submersible is at top speed.
The prototype currently lacks control surfaces, which required the trials to be conducted
using a guide string. This string runs through a tube aﬃxed to the dorsal ﬁn; there is
considerable play in the interface between the tube and string, and small rotations of the
rigid hull are visible during operation in the pitch, roll and yaw axes. The average speed
103

of the submersible is measured by analyzing 30 frame/sec video taken during operation;
the guide string has been marked at two locations 9.14 m apart, and the times at which the
submersible crosses these marks can be determined, and the speed calculated. The tests were
conducted relatively near the surface of the pool, to enable wireless control of the throttle
at the start of the run. In order to limit the loss of energy due to free surface eﬀects, the
mass and center of mass of the rigid head were adjusted such that the measured portion of
the run was entirely below the surface.
The average speed of each conﬁguration is given in Table 5.1. The reader will note that
only the speed for each conﬁguration is cited, rather than individually listing the thrust and
drag of each conﬁguration. This is consistent with the approach taken in the prior chapter;
for a self-propelled swimming body the sources of drag and thrust cannot be conveniently
separated, unlike a conventional propeller-driven craft.
Table 5.1: Mean speeds of the three conﬁgurations in experiments
Conﬁguration Trials Speed [m/s] Speed [BL/s]
Rigid Tube
5
0.918
1.176
Rigid Tail
10
0.878
1.126
Flexible Tail
4
0.889
1.139

The relatively low number of trials is the result of selecting trials for which the prototype
is traveling in a relatively straight line, and not interacting with the guide string; these runs
can be determined after the fact, by using the video data. Figure 5.3 provides essentially
the same data as Table 5.1 with included error bars. It is clear that the general trend found
in simulation (Fig.4.6) is veriﬁed by experiment; a ﬂuttering ﬂexible tail has superior performance to a rigid tail of identical dimension, but does not necessarily exceed the performance
of a rigid tube. The speed of each conﬁguration in simulation is also similar to the value
104

Speed [m/s]

0.95

0.9

0.85

0.8
Rigid

Flexible

Tube

Figure 5.3: Mean speeds of the conﬁgurations listed in Table 5.1. The error bars express the
95% conﬁdence intervals.
found by experiment, though this is to be expected, since experimental values were used
to calibrate the drag coeﬃcients used in simulations. The ≈ 0.33 sec. period of oscillation
varies slightly from the 0.306 sec. period found in simulation. However, the shape of the tail
is similar, in that the tail appears to conform to a second mode waveform. Note also that
the experimental period is diﬃcult to determine precisely from the relatively slow video used
in these experiments.

5.3
5.3.1

Eﬀect of Time-Varying Internal Fluid Velocity
Additional Terms

The equations of motion can be modiﬁed to account for a time-variable conveyed ﬂuid
velocity. The following changes will be made in the derivation presented in Section 4.4:
105

Head Internal Fluid: Ω1

¨
Replace: rxy = 0

¨
˙ i
With: rxy = Uhˆ

Tail Internal Fluid: Ω2
¨
Replace: rxy = (¨ + 2U y ′ + U 2 y ′′ )ˆ
y
˙
j

¨
˙i
˙ j
With: rxy = Uˆ + (¨ + 2U y ′ + U 2 y ′′ + U y ′ )ˆ
y
˙

Expanded per Eqs.(4.18) and (4.20), the above replacements require that the following terms
be added to the scalar equations of motion, Eqs.(4.45 ), (4.46) and (4.47):

Eq.(4.45) :

˙
m1 U˙h + m2 U

Eq.(4.46) :

˙ y(L)
m2 U
L

Eq.(4.47) :

m ˙ L
m1 rcy U˙h + 2 U
(x − rcx )y ′ − (y − rcy ) dx
L
0

The equation of motion for the ﬂexible beam must also be expanded to account for the
time-variable internal ﬂuid velocity. The rate of change of momentum used in Eq.(4.55) is
therefore:

∂V
¨
˙ j
˙i
˙
= R0 + Uˆ + y + 2U y ′ + U 2 y ′′ + Uy ′ ˆ + θ x ˆ − y ˆ + 2θ U ˆ − (y + Uy ′ ) ˆ
¨
˙
j ¨ j
i
˙
i
∂t

With the addition of these terms, Eq.(4.58) is replaced with the following expression:

∂p
∂p
∂T
− A2 2 − A6 4 =
∂x
∂x
∂x

xx
mt + m2 + m4 kt
L
106

xx
m4 kt
m2 ˙
x
˙
ax +
U+
Ue − Dt
L
L

Equation (4.60) is replaced by
yy
2
m2 2 m4 kt
∂
∂4y
∂y
y
2 ∂ y
+
(5.1)
Dt − EI
(T − A2 p2 − A6 p4 )
=
U +
Ue
∂x
L
L
∂x4 ∂x
∂x2
yy
yy
yy
˙
˙
m2 U + m4 kt Ue ∂ 2 y
m2 U + kt m4 Ue ∂y mt + m2 + m4 kt
∂2y
+2
+ ay
+
+
L
∂x∂t
L
∂x
L
∂t2
where the tension-pressure gradient term in Eq.(5.1) is replaced by:

∂p
∂p
∂T
∂y
=
− A2 2 − A6 4
∂x
∂x
∂x ∂x
xx
xx
mt + m2 + m4 kt
m4 kt
m ˙
x ∂y
˙
Ue − Dt
ax + 2 U +
L
L
L
∂x
xx
2
L mt + m2 + m4 k xx
t a + m2 U + m4 kt U − D x dx + D b ∂ y
˙
˙e
−
x
t
t ∂x2
L
L
L
x

∂
∂x

(5.2)

Substitution of Eq.(4.58) into Eq.(5.1) yields the modiﬁed equation of motion:

yy
m2 2 m4 kt
∂4y
2
+
U +
Ue
EI
L
L
∂x4
xx
2
L mt + m2 + m4 k xx
t a + m2 U + m4 kt U − D x dx − D b ∂ y
˙
˙e
+
x
t
t ∂x2
L
L
L
x
yy
xx
xx
m4 (kt − kt )
mt + m2 + m4 kt
x ∂y
˙
Ue + Dt
ax −
+ −
L
L
∂x
yy
yy
m4 kt
mt + m2 + m4 kt
m2
∂2y
∂2y
y
+ ay − Dt = 0
+2
U+
Ue
+
L
L
∂x∂t
L
∂t2

(5.3)

˙ ∂y term appears for U . This is consistent
˙
It is interesting to note that no analog of the Ue
∂x
with a previous derivation of the equation of a ﬂuid-conveying with time variable internal
ﬂow [32]. For internal ﬂow, there is no diﬀerence between the eﬀective mass acting in the x
107

yy
xx
and y directions, meaning that the (kt − kt ) term is eﬀectively zero for internal ﬂow.

5.3.2

Numerical Simulations

˙
˙
Addition of the terms associated with U and Uh requires that two additional states be added
to the discretized equations of motion, such that

˙
˙ ˙ ˙ ˙
U = fu (X0 , Y0 , θ, yi , X0 , Y0 , θ, yi , t)
˙
˙ ˙ ˙ ˙
Uh = fh (X0 , Y0 , θ, yi , X0 , Y0 , θ, yi , t)

˙
˙
Note that if U and Uh are explicit functions of only time, no additional states need to be
added to the discretized equations, since U and Uh can be computed a priori. We will also
make another assumption:
A10. At any time, the mass ﬂux through the inlet Γ1 is equal to that through the control
surface at the tip of the tail, Γ2 . This precludes ﬂuid storage within the rigid head of
˙
˙
the submersible, and implies that Uh = (A2 /A1 )U.
˙
The functional description of U is as follows:

˙
U = αy(kL)
˙

U(0) = 4.0

(5.4)

where α is a scaling parameter, and y(kL) is the velocity of a given point on the ﬂexible
˙
tail. Functions of this type have several favorable properties, not least of which is the ease
of implementation within the simulation. While we have yet to do so, it is also easy to
experimentally implement the sensing required for this type of function, since it requires a
single point measurement at a point on the beam. This is in contrast to a function requiring
108

knowledge of one or more spatial derivatives, or an approximation of the integral of the beam
position.
It is easy to understand why a function in this form can cause the vehicle to turn; since
y(kL) is signed, the thrust from the jet may be unbalanced with respect to the centerline of
˙
the vehicle. The degree of imbalance is related to the point k at which the measurement is
taken, because of the traveling nature of the beam’s waveform. Essentially, by moving the
point k nearer to the head of the vehicle, we can mimic a phase shift between the observed
˙
value of y(kL) and U.
˙
In the ﬁgures which follow, the tail is the “baseline” tail described in the previous section:
L = 0.34, h = 0.077. In addition, the value for U has been subjected to upper and lower
saturation such that 3.8 ≤ U ≤ 4.2; removal of these limits aﬀects neither system stability
nor the qualitative results, but it is diﬃcult to keep the average velocity between cases the
same without imposing them. Furthermore, similar limits are likely to be imposed by any
physical pump.
Figure 5.4 gives a comparison of the constant velocity case where U = 4.0 and the nonconstant velocity case where α = 40, k = 0.75. The total rotation of the vehicle θ is given in
the left axes, and the right axes depict the forward speed of the vehicle Ue and the conveyed
ﬂuid velocity U(t). It is clear from the left axes of Fig.5.4 that the U = U(t) case is subject
to a net moment over the course of its travel, making a clockwise turn of ≈ 1 rad. relative to
the U = const case. The middle right axes show the forward speed of the vehicle for both
cases, and Ue for both cases can be seen to be very similar; the top right axes show that
there are larger ﬂuctuations in the forward speed of the U = U(t) case. The conveyed ﬂuid
velocity is shown in the bottom right axes, and the increase in U(t) amplitude and eventual
109

0.895
Ue [m/s]

0.4

0.2
U=const

−0.2

−0.4

Ue [m/s]

[rad]

0.88
0.875
27

0

θ

0.89
0.885

5

28

10

28.5

15

20

29

25

30

U=U(t)
[m/s]

−0.6

U

−0.8

−1
0

0.9
0.8
0.7
0.6
0.5
0

27.5

5

10

15
20
time [s]

25

30

4.2
4
3.8
0

2

4
time

6
[s]

8

10

Figure 5.4: Comparison of the constant-velocity U = 4.0 case to non-constant velocity case
where α = 40, k = 0.75

saturation is clearly visible. Note that the time scale of the U(t) plot is 1/3 that of the Ue (t)
plot.
Figure 5.5 compares the behavior of three values of k: 0.25, 0.75 and 1.0. The left axes
show θ for each case. The right set of axes shows the orientation of each at t = 20 sec.,
marked with a dashed line on the θ axes.
The diﬀerent behavior of each value of k is evident from each case’s θ curve. Measuring
the beam velocity at diﬀerent points places U(t) in diﬀerent phase with respect to the
orientation of the tip of the tail, due to the traveling wave nature of the tail’s waveform. A
measurement point of k = 0.75 corresponds to a large negative moment, k = 1.0 a smaller
negative moment, and the direction of turn can be reversed by measuring at k = 0.25. A
physical realization of this type of control is likely to simply reverse the sign of Eq.5.4 in
order to reverse the direction of the moment, and it is likely that such a realization will use
110

0.4
k=0.25
0.2
k=0.25

[rad]

k=1.0

−0.2

θ

0

−0.4

k=1.0

−0.6

k=0.75

k=0.75

−0.8

−1
0

5

10

15
20
time [s]

25

30

Figure 5.5: Comparison of three values of k: 0.25, 0.75 and 1.0. The scaling parameter
α = 40 for all cases.

a measurement point near x = L, to take advantage of the larger deﬂections there.
The sequence of images shown in Fig.5.6 have been acquired over one cycle of the tail
where U = U(t) with λ = 0.75 - note that the frame of reference pictured is the non-inertial
(xy) frame. Compared to Fig.4.7, there is a distinct asymmetry to the tail’s motion; this
asymmetry is most clear at the tip of the tail, which spends a larger fraction of the cycle
on the right half of the submersible’s centerline. This asymmetry is the likely reason for the
turning behavior - as the angle of surfaces Γ2 and Γ6 change throughout the cycle, ﬂuid
which enters through the inlet surfaces Γ1 and Γ3 is turned with respect to the centerline.

111

Figure 5.6: An image sequence of the submersible during a turn. Note that the dashed line
is collinear with the non-inertial x axis.

112

Chapter 6
Concluding Remarks and
Future Work

6.1

Dynamics of a Pipe Conveying Fluid with a NonUniform Velocity Proﬁle

Chapter 2 assessed the dynamics of cantilever pipes conveying ﬂuid with a fully developed
non-uniform velocity proﬁle. The equation of motion derived in that section is tractable,
requiring only the use of a single empirical parameter µ to account for the dependence of
ﬂuid momentum on the square of the ﬂuid velocity. Previous analyses made the assumption
of a uniform velocity proﬁle, implicitly assuming µ to be unity. While this assumption is
reasonable at high Reynolds number, µ = 1 is the minimum value possible, approached only
in the limit of inﬁnite Reynolds number. That is, the momentum ﬂux for real ﬂuid ﬂows
will always be greater than the uniform case, and the uniform proﬁle assumption becomes
monotonically less accurate as the Reynolds number decreases.
113

It was shown that the dependence of u, the non-dimensional velocity relevant to the
stability of the pipe, is not uniquely dependent on the Reynolds number. Since the Reynolds
number determines the velocity proﬁle, a uniform velocity proﬁle may not be generally
assumed near the onset of ﬂutter instability, despite the requirement of high u. This is
particularly relevant for laminar ﬂow, where the value of µ reaches its maximum value; as
the scale of “engineering applications” shrinks, the dynamics of pipes conveying laminar ﬂow
with high nondimensional velocity u may become relevant.
The stability characteristics of a sample pipe were assessed with our updated model.
In this pipe, the ﬂuid became turbulent at relatively low u with the eﬀect that µ ≈ 1.05
when u approached values necessary to achieve ﬂutter instability. This proximity to the
uniform model causes similarity between the predictions over much of the parameter space,
a similarity which is to be expected since the predictions of the uniform model have been
qualitatively veriﬁed by experiment. There are signiﬁcant diﬀerences in certain sensitive
regions of the parameter space, such as near β = 0.3 and β = 0.7. A pipe with dimensions
such that the conveyed ﬂuid is laminar at the onset of ﬂutter instability exhibited markedly
diﬀerent stability characteristics from the uniform case. Both the critical velocity and critical
frequency were signiﬁcantly diﬀerent over large regions of the parameter space.

6.2

Simpliﬁed Dynamics of a Fluid-Conveying, FluidImmersed Pipe Aﬃxed to a Rigid Body

Chapter 3 derived the equations of motion for an immersed ﬂuid-conveying tail aﬃxed to a
rigid body. It was shown that both the cantilever and free pipe can be expressed as special
114

cases of the rigid body boundary condition. The neutral stability curve in ui -ue space was
computed over a wide range of values of the rigid body mass fraction. The serpentine nature
of these curves renders it diﬃcult to make statements about the entire parameter space;
nonetheless, some local trends can be observed. At low values of ue , a tail aﬃxed to a rigid
body of larger mass become unstable at lower values of ui , but this trend is reversed for higher
values of ue . The system’s stability was also found to be largely insensitive to additional
rigid body mass above µ = 25, i.e., the system’s behavior is largely indistinguishable from
that of a cantilever.
Estimates of the sign of the thrust produced by the ﬂuttering tail and the eﬃciency of
that thrust have also been computed via methods ﬁrst proposed by Lighthill [23]. Although
this analysis is less robust than that proposed in Chapter 4, some general statements about
the suitability of this type of propulsor can be made. A relatively high value of ui is required
to drive a thrust-producing ﬂutter instability. Although a system of this type may permit
ﬂutter instability for ui = 0, instability driven entirely by the external ﬂow does not provide
a thrust-producing motion. This ui = 0 limit case is intuitively obvious, since a ﬂapping
ﬂag does not produce thrust. The regions of the neutral stability curves where thrust is
produced therefore have boundaries primarily deﬁned by ui , since the internal ﬂuid ﬂow is
the only energy input to the system.
It is also clear that a less massive rigid head leads to thrust-producing instability at
lower values of ui , though these motions occur at lower forward velocity. Note also that
the eﬃciency of the produced thrust is relatively insensitive to ue over a wide range of ue
values. Similar observations have been made for live ﬁsh [34], which move with a waveform
reminiscent of the traveling waveform generated by a ﬂuid-conveying pipe. It is therefore
115

heartening that one of the great advantages of ﬁsh-like propulsion is preserved.
The analysis proposed in this chapter was conﬁned to the mass of the rigid body, but it
is possible to investigate the role of many other parameters for a vehicle of this type. Besides
the parameters identiﬁed in this work, a terminal nozzle might be used at the free end, or
a ﬁsh-like planform used for the shape of the tail. Analysis of the latter system must be
performed numerically, and will likely be an adjunct to later work along the lines of that
presented in Chapter 4.

6.3

General Dynamics of a Fluid-Conveying,
Fluid-Immersed Pipe Aﬃxed to a Rigid Body

Chapter 4 presented the planar equations of motion for a ﬂuid-conveying submersible with
a ﬂuttering tail without many of the limiting assumptions used in Chapter 3. In particular,
the only assumption made regarding the accelerations or rotations of the non-inertial frame
attached to the submersible are that the accelerations and rotations occur within a single
plane. This derivation therefore comprises the scalar equations for the motion of this noninertial frame as well as the equation of motion for the ﬂuid-conveying tail within the frame.
These equations appear here for the ﬁrst time in the academic literature likely because most
applications for ﬂuid-conveying pipes, such as pipeline vibrations or the simple “garden hose”
problem, do not require analysis in a non-inertial frame.
These equations are solved for a variety of tail geometries, and one geometry is presented
which produces thrust in excess of a rigid tube of equal length; this provides evidence of what
our research group has termed “synergistic propulsion”, in which there is net thrust provided
116

by the ﬂuttering motion of the tail, in addition to the thrust of the ﬂuid jet. This particular
ﬁnding has yet to be veriﬁed by the experimental work described in Chapter 5, though the
speed of a device with a ﬂuttering tail has been shown to exceed that of the device with a
dimensionally identical rigid tail. This ﬁnding is the same as that of Paidoussis [28], since
that work compared the ﬂuttering and non-ﬂuttering tails by placing an additional stiﬀener
on the ﬂexible tail.
The methods laid out in this chapter are likely to be the basis for much of the future
work done on this topic; indeed one such “future” topic, control of a ﬂuttering ﬂuid-conveying
submersible, is touched upon in Chapter 5 of this communication. This type of empiricallybased numerical analysis is perhaps a bridge between the analytical methods typiﬁed by
the work in Chapter 3 and a full-blown CFD simulation. In addition to further work on
control, extension of these methods to non-uniform tail geometries is planned for the future.
A ﬁsh-like planform, with a large reduction in area near the tip of the tail, has been shown
[25] to be highly eﬃcient for producing thrust on the basis of slender-body theory. A beam
with a ﬁsh-like planform will be relatively weak at a location near its its end, and then much
stiﬀer in that portion which mimics the caudal ﬁn. This should lead to a large, out-of-phase
angle between this “ﬁn” and the forward portion of the tail, a condition which is seen in
swimming ﬁsh.
The large angles and large amplitudes associated with non-uniform tail geometries require
that the eﬀects of large displacements be accounted for. The equations of motion of a ﬂuidconveying Euler elastica have been presented before [4] [39], though neither of the cited
works considered either external ﬂow or a non-uniform pipe/beam geometry. Extension to
these conditions is planned for future work.
117

6.4

Experimental Studies and Control

Chapter 5 describes a prototype version of a ﬂuttering ﬂuid-conveying submersible that
was produced at Michigan State University. Much of the credit for the prototype’s design
and construction is due to Mr. Paul Streﬂing, and his M.S. thesis [40] may be used as a
more complete source of information about the platform. Preliminary ﬁndings qualitatively
verify the model proposed in Chapter 4. A great deal of work remains to be done on the
experimental hardware. Communication with the device as currently constructed is not
possible through water, which means that it is diﬃcult test in a way which removes freesurface eﬀects. Underwater communication is likely to be added soon, as is buoyancy control,
to permit running at greater depth.
The equations of motion derived in Chapter 4 have also been extended to permit the
use of a time-variable internal ﬂuid velocity, and a simple function for this velocity has been
shown to turn the submersible. This function requires the velocity measurement of a single
point on the tail; changing this measurement point leads to a phase shift in internal ﬂuid
velocity which can be exploited to turn the submersible in both directions. The potential
ability to maneuver by changing the internal ﬂow rate is a useful feature, since it would allow
control of the submersible without additional control surfaces, such as ﬁns.
The function proposed for the internal ﬂuid velocity produced reasonable turning rates,
but it is likely that there is signiﬁcant potential for improvement by using information about
additional states of the system. In particular, Eq.(5.4) contains no information about the
˙
global coordinates X, Y, θ in its deﬁnition for U (t); knowledge of the state of the non-inertial
frame will certainly be required for trajectory tracking. Note also that the turn shown in
Fig.5.4 takes place over many oscillations of the tail, and the image sequence shown in Fig.
118

5.6 shows a tail which is oscillating asymmetrically, but not wildly so. It is therefore likely
that the device is nowhere near its limits of its turning performance; a functional description
˙
for U (t) which leads to a more asymmetric oscillation of the tail is likely to increase the
turning rate of the vehicle.

119

BIBLIOGRAPHY

120

BIBLIOGRAPHY

[1] PV Alvarado and K Youcef-Toumi. Design of machines with compliant bodies for
biomimetic locomotion in liquid environments. ASME J. Dynamic Systems, Measurement, and Control, 128:3–13, 2006.
[2] JM Anderson, K Streitlien, DS Barrett, and MS Triantafyllou. Oscillating foils of high
propulsive eﬃciency. Journal of Fluid Mechanics, 360:41–72, 1998.
[3] H. Ashley and G. Haviland. Bending vibrations of a pipeline containing ﬂuid. Journal
of Applied Mechanics, 17:229–232, 1950.
[4] AK Bajaj, PR Sethna, and TS Lundgren. Hopf bifurcation phenomena in tubes carrying
ﬂuid. SIAM Journal of Applied Mathematics, 39:213–230, 1980.
[5] DS Barrett, MS Triantafyllou, DKP Yue, and M Wolfgang. Drag reduction in ﬁsh-like
motion. J. Fluid Mechanics, 392:183–212, 1999.
[6] GK Batchelor. An Introduction to Fluid Mechanics. Cambridge University Press (London), 1967.
[7] R Benedict. Fundamentals of Pipe Flow. John Wiley and Sons, 1980.
[8] RB. Bhat and H. Wagner. Natural frequencies of a uniform cantilever with slender tip
mass in the axial direction. Journal of Sound and VIbration, 46(2):304–307, 1976.
[9] CE Brennan. A review of added mass and ﬂuid inertial forces. Technical report, Naval
Civil Engineering Laboratory, 1982.
121

[10] Z Chen, S Shatara, and X Tan. Modeling of biomimetic robotic ﬁsh propelled by an
ionic polymer metal composite caudal ﬁn. IEEE/ASME Transactions on Mechatronics,
15(3):448–459, 2010.
[11] HK Cheng and L Murillo. Lunate tail swimming propulsion as a problem of curved
lifting line in unsteady ﬂow. Journal of Fluid Mechanics, 143:327–350, 1984.
[12] FE Fish. The myth and reality of gray’s paradox: implication of dolphin drag reduction
for technology. Bioinspiration & Biomimetics, 1:R17–R25, 2006.
[13] J Gray. Studies in animal locomotion. i. the movement of ﬁsh with special reference to
the eel. Journal of Experimental Biology, 10:88–104, 1933.
[14] J. Gray. Studies in animal locomotion vi. the propulsive powers of the dolphin. Journal
of Experimental Biology, 13:192–199, 1936.
[15] R.W. Gregory and M.P. Paidoussis. Unstable oscillation of tubular cantilevers conveying
ﬂuid. i. theory. Proceedings of the Royal Society (London), 293:512–527, 1966.
[16] S. Guo, T. Fukuda, and K. Asaka. A new type of ﬁsh-like underwater microrobot.
IEEE/ASME Transactions on Mechatronics, 8(1):136–141, 2003.
[17] MJ Hannoyer and MP Paidoussis. Instabilities of tubular beams simultaneously subjected to internal and external axial ﬂows. ASME Journal of Mechanical Design,
100:328–336, 1978.
[18] KA Harper, MD Berkemeier, and S Grace. Modeling the dynamics of spring-driven
oscillating-foil propulsion. IEEE J. Oceanic Engineering, 23:285–296, 1998.
[19] AM Hellum, R Mukherjee, and AJ Hull. Dynamics of pipes conveying ﬂuid with nonuniform turbulent and laminar velocity proﬁles. Journal of Fluids and Structures, 26
(5):804–813, 2010.
[20] AM Hellum, R Mukherjee, and AJ Hull. Flutter instability of a ﬂuid-conveying ﬂuidimmersed pipe aﬃxed to a rigid body. Journal of Fluids and Structures, In Press, 2011.
[21] SF Hoerner. Fluid-Dynamic Drag. Self Published, 1958.
[22] GV Lauder, EJ Anderson, J Tangorra, and Madden PGA. Fish biorobotics: kinematics
and hydrodynamics of self-propulsion. Journal of Experimental Biology, 210:2767–2780,
2007.
122

[23] M.J. Lighthill. Note on the swimming of slender ﬁsh. Journal of Fluid Mechanics,
9:305–317, 1960.
[24] M.J. Lighthill. Hydrodynamics of aquatic animal propulsion. Annual Review of Fluid
Mechanics, 1:413–446, 1969.
[25] M.J. Lighthill. Aquatic animal propulsion of high mechanical eﬃciency. Journal of
Fluid Mechanics, 44:265–301, 1970.
[26] MJ Lighthill. Large-amplitude elongated-body theory of ﬁsh locomotion. Proceedings
of the Royal Society (London), B179:125–138, 1971.
[27] RL Mcmasters, CP Grey, JM Sollock, R Mukherjee, A Benard, and AR Diaz. Comparing the mathematical models of lighthill to the performance of a biomimetic ﬁsh.
Bioinspiration & Biomimetics, 3:1–8, 2008.
[28] MP Paidoussis. Hydroelastic icthyoid propulsion. AIAA Journal of Hydronautics, 10:30–
32, 1976.
[29] MP Paidoussis. Marine propulsion apparatus. US Patent 4,129,089, December 1978.
[30] M.P. Paidoussis. Fluid-Structure Interactions: Slender Structures and Axial Flow, Volume 1. Academic Press, 1998.
[31] MP Paidoussis. Fluid-Structure Interactions: Slender Structures and Axial Flow, Volume 2. Academic Press, 2004.
[32] M.P. Paidoussis and N.T Issid. Dynamic stability of pipes conveying ﬂuid. Journal of
Sound and Vibration, 33:267–294, 1974.
[33] MC Potter and JF Foss. Fluid Mechanics. Great Lakes Press, 1982.
[34] J. Rohr and F. Fish. Strouhal number and optimization of swimming by odontocete
cetaceans. SSC San Diego Biennial Review, pages 175–178, 2006.
[35] S Saimek and PY Li. Motion planning and control of a swimming machine. Int. J. of
Robotics Research, 23:27–53, 2004.
[36] J Saroudis. Propulsion by oscillating plates. B. Eng. Honours Thesis, Department of
Mechanical Engineering, McGill University, Montreal, Canada, 1974.
123

[37] WW Schultz and PW Webb. Power requirements of swimming: do new methods resolve
old questions? Integrative and Comparative Biology, 42:1018–1025, 2002.
[38] AP Seiranyan. Collision of eigenvalues in linear oscillatory systems. Journal of Applied
Mathematics and Mechanics, 58 (5):805–813, 1994.
[39] M Stangl, J Gerstmayer, and H Irschik. An alternative approach for the analysis of nonlinear vibrations of pipes conveying ﬂuid. Journal of Sound and Vibrations, 310:493511,
2008.
[40] P.C. Streﬂing. Experimental platform for studying the behavior and performance of
a submersible propelled by a ﬂuid-conveying ﬂuttering tail. Master’s thesis, Michigan
State University, 2011.
[41] L. Tang, M.P. Paidoussis, and J. Jiang. Cantilevered ﬂexible plates in axial ﬂow: Energy
transfer and the concept of ﬂutter-mill. Journal of Fluids and Structures, 326:263–276,
2009.
[42] GI Taylor. Analysis of the swimming of of long and narrow animals. Proc. of the Royal
Society (London), A214:158–183, 1952.
[43] MS Triantafyllou, GS Triantafyllou, and DKP Yue. Hydrodynamics of ﬁshlike swimming. Annual Review of Fluid Mechanics, 32:33–53, 2000.
[44] T. Yao-Tsu Wu. Hydromechanics of swimming propulsion. part 3. swimming and optimum movements of slender ﬁsh with side ﬁns. Journal of Fluid Mechanics, 46:545–568,
1971.
[45] J Yu, S Wang, and M Tan. A simpliﬁed propulsive model of bio-mimetic robot ﬁsh and
its realization. Robotica, 23:101–107, 2005.

124