Cozens PD 1987 PHD Thesis

NUMERICAL MODELLING OF THE
ROLL DAMPING OF SHIPS DUE TO
VORTEX SHEDDING
by
Paul D e n n i s C o z e n s , M.A. (Cantab.)
D e p a r t m e n t of A e r o n a u t i c s
Imperial C o l l e g e
P r in ce C o n s o r t Road
L o nd on SW7 2BY
A th es is s u b m i t t e d for the d e g r e e of
Do ct or of P h i l o s o p h y of the U n i v e r s i t y of London
and for the
D i pl om a of M e m b e r s h i p of the Imperial Co llege
September, 1987
2
Flow visualization of the vortex shedding from the bilges of a rolling
rectangular cross-section barge with bilge keels

3
SUMMARY
The vortex shedding from the bilges of rolling ships has been simulated
using three computational models.
The first model involves the application of the Discrete Vortex Method to
the simulation of oscillatory flow around infinite wedges, in particular
the square isolated edge with bilge keel and/or edge rounding. Computed drag
forces were included, via a matching technique, in a ship motion response
computer program, so as to yield a method for predicting the roll damping of
transportation barges due to vortex shedding. Comparisons with experiment
for rectangular cross-section barges with various bilge geometries revealed
that whilst the method was successful at computing roll damping for barges
with sharp bilges, it overestimated damping increasingly at higher bilge
radii.
The second model is an extension of the basic Discrete Vortex Method to

predict the vortex shedding from the bilges of rolling vessels of quite
general cross-section subject to certain restrictions. Computations were
performed for a ship representative of current research interest. A tendency
to overestimate roll damping was explained in terms of weaknesses in the
Discrete Vortex Method itself and in its application to shedding from
three-dimensional bodies.
Finally, a more advanced computer model of the full two-dimensional unsteady

Navier-Stokes equations of motion has been developed which has enabled the
prediction of the roll damping of barges with rounded bilges. This model
represents, through its novel solution approach, a considerable advance in
the simulation of vortex shedding especially for configurations where the
separation point is required to be calculated as part of the computational
solution. Comparisons with experiment for a variety of cases clearly
established the validity of the method.
Flow visualisation experiments were undertaken to study the vortex shedding

processes from barges with bilge keels and/or bilge rounding. It is shown
that both the Discrete Vortex Method and the full Navier-Stokes solver
predict flow patterns for oscillatory flow which correspond closely to
experiment.
The three methods described here are capable of quite general application to
many types of vortex shedding problems. Suggestions for possible future
areas of work are described.
4
ACKNOWLEDGEMENTS
It is with a deep sense of gratitude that I write these
acknowledgements. The very considerable rewards of the research work
itself have been added to immeasurably by the kindness of all those at
Imperial College who have made it possible. In the first instance, I
count myself most fortunate that I should have been supervised by
Professor P. W. Bearman and Dr. J. M. R. Graham who have both been an
unquenchable source of ideas, inspiration and guidance.
I am most grateful to BMT Ltd. and the Marine Technology Directorate of
the Science and Engineering Council for funding this project. I would
also like to thank Dr. R. G. Standing and Mr. A. K. Brook of BMT Ltd.
for their advice and support.
I am indebted to Mr. J. O'Leary for his patient assistance in setting up
the experiments, Mrs. R. Fairhurst for typing many parts of my thesis,
and Mr. J. D. Powell for proof reading.
My thanks also to Dr. M. J. Downie, Dr. J. P. Felix and Mr. P. Kalkanis
for their useful advice and the many helpful discussions, and to all my
friends and colleagues who have patiently helped me through these last
few months.
Last, but by no means least, my heartfelt thanks to my wife, Elizabeth,
for the support without which this thesis could never have even been
completed.
5
CONTENTS
Page
SUMMARY 3
ACKNOWLEDGEMENTS 4
CONTENTS 5
NOTATION 11
1. INTRODUCTION 14
PREFACE 14
1.1 The Prediction of the Motions of a Vessel in a Seaway 15
1.1.1 The prediction of vessel response in waves using
strip theory 15
1.1.2 The roll response of a vessel at zero forward speed 16
1.2 The Discrete Vortex Method 19
1.2.1 Early development of the Discrete Vortex Method 19
1.2.2 The application of the Cloud-in-Cell technique to the
Discrete Vortex Method 24
1.3 The Solution of the Full Two-Dimensional Unsteady
Navier-Stokes Equations 26
1.3.1 Eulerian solutions of the Navier-Stokes equations 27
1.3.2 Lagrangian solutions of the Navier-Stokes equations 29
1.3.3 Comparison of Eulerian and Lagrangian schemes 31
1.3.4 Hybrid Lagrangian/Eulerian solutions of the
Navier-Stokes equations 31
1.4 The Application of Vortex Methods to Vortex Shedding
from Vessel Hulls 33

6
2. DESCRIPTION OF THE PROCESSES OF VORTEX SHEDDING
IN OSCILLATORY FLOW 37
2.1 Experimental Method 38
2.1.1 Experimental apparatus 38
2.1.2 Test procedure 38
2.1.3 Flow visualization techniques 40
2.2 Photographs of the Vortex Shedding Processes 41
3. THE DISCRETE VORTEX METHOD 46
3.1 Transformation 46
3.1.1 Plain isolated edge transformation 47
3.1.2 Plain isolated edge with bilge keel transformation 48
3.1.3 Full numerical solution of the Schwarz-Christoffel
transformation 49
3.2 Basic Concepts of the Discrete Vortex Method 54
3.3 Treatment of the Nascent Vortex 56
3.3.1 Satisfying the Kutta Condition 56
3.3.2 Determination of the position of the nascent vortex 57
3.4 A Scheme for Amalgamating Vortices 59
3.5 Starting the Flow Pattern 61
3.6 Decaying the Vortex Pairs 61
3.7 Calculation of Pressure Coefficient 63
3.8 Calculation of Vortex Force 64
3.9 Extension of the Discrete Vortex Method to Model the
Vortex Shedding from Rolling Ships 64
3.9.1 Vessel in pure sway 64
3.9.2 The representation of roll 65

7
3.9.3 Method of solution for the velocity distribution
due to a source distribution 67
3.9.4 Application of the theory to a rolling vessel 67
3.9.5 Calculation of the vortex roll damping moment 70
3.9.6 Range of validity of the method 71
4. RESULTS FROM THE DISCRETE VORTEX METHOD AND COMPARISON WITH
EXPERIMENT FOR ISOLATED EDGE FLOWS 73
4.1 Assessment of the Sensitivity of the Method to Changes
in the Run Parameters 73
4.1.1 Effect of vortex decay rate 73
4.1.2 Effect of the position of the nascent vortex 74
4.1.3 Effect of time step size and amalgamation rate 75
4.2 Results for Square Edge with Bilge Keel and
Optional Rounding 77
4.2.1 Vortex positions and pressure distributions 77
4.2.2 Comparison between theoretical and experimental
vortex size 78
4.2.3 Time dependent force coefficient 79
4.2.4 Drag and inertia coefficients 79
4.3 Results for Square Edge withRounding 81
4.3.1 Vortex positions and pressure distributions 82

8
5. APPLICATION OF THE DISCRETE VORTEX METHOD TO ROLLING
SHIPS AND BARGES 85
5.1 Barge Motions and the BMTIMP Barge Rolling Computer
Program 85
5.1.1 Treatment of the roll damping contributions
by BMTIMP 87
5.1.2 Matching the isolated edge flow to flow round the
rolling barge 88
5.2 Results for a Rectangular Cross-Section Barge with
Bilge Keels 91
Rounded Bilges 93
Rounded Bilges and BilgeKeels 95
5.5 Results for a Typical Ship inRoll Motion 96
5.5.1 The roll damping of a circular cylinder with
bilge keels 97
5.5.2 The Fisheries Protection Vessel "Sulisker" 98
5.5.3 Vortex positions and velocity vectors 99
5.5.4 Roll damping results 103
6. DEVELOPMENT OF A HYBRID MESH / POINT VORTEX NAYIER-STOKES
SOLVER FOR THE PREDICTION OF THE FLOW AROUND ISOLATED EDGES 109
6.1 Fundamental Equations 110
6.1.1 Primitive equations 110
6.1.2 Stream function and vorticity transportequations 111
6.1.3 Boundary conditions 112

9
6.2 Finite Difference Forms 112
6.3 Outline of the Basic Inviscid Cloud-In-Cell Method 114
6.4 Outline of the Hybrid Moving Vortex Diffusive Method 116
6.5 Transformation 120
6.5.1 Utilisation of the general numerical Schwarz-
Christoffel transformation 120
6.5.2 Setting up the mesh 120
6.6 Treatment of the Diffusion Term in the Solution of the
Navier-Stokes Vorticity Transport Equation 123
6.6.1 Solution of the diffusion equation 123
6.6.2 Vorticity boundary condition 125
6.7 Solution of the Poisson Equation 128
6.7.1 Stream function boundary condition 130
6.8 Determination of Mesh Velocities 131
6.9 Convection of Vortices 132
6.10 Interpolation Between the Point Vortices and the Mesh 133
6.11 Transfer of Mesh Circulation onto the Point Vortex System 133
6.11.1 Amalgamation of vortices 134
6.12 Calculation of Vortex Force 135
7. RESULTS FROM THE HYBRID MOVING VORTEX DIFFUSIVE METHOD 138
7.1 Test Case of an Oscillatory Boundary Layer over a
Flat Plate 138
7.1.1 Analytical solution for an oscillatory boundary layer 139
7.1.2 Comparison of the analytical solution with the
computed results 142
7.2 Vortex Shedding from an Isolated Square Edge in
Uniform Flow 144

10
7.2.1 Definition of Reynolds number in uniform flow 144
7.2.2 Vortex positions and streamlines 145
7.3 Results for Oscillatory Flow Past a Square Edge
with Rounding 146
7.3.1 Definition of Stokes number in oscillatory flow 146
7.3.2 Choice of computational mesh 147
7.3.3 Choice of run parameters 148
7.3.4 Flow results 152
7.3.5 Time dependent vortex force coefficient 158
7.3.7 Effect of variation of Stokes number 162
7.3.8 Sensitivity to mesh density and timestep size 162
7.4 Application of the Hybrid Moving Vortex Diffusive Method
Results to the Roll Damping of a Rectangular Cross-Section
Barge with Rounded Bilges 163
8. CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK 165
8.1 Conclusions Concerning the Discrete Vortex Method 165
8.1.1 Isolated edge results 165
8.1.2 General hull cross-section results 166
8.1.3 General conclusions 169
8.2 Conclusions Concerning the Hybrid Moving Vortex
Diffusive Method 170
8.3 Conclusions Concerning the BMTIMP Computer Program 172
REFERENCES 175
FIGURES 190
11
NOTATION
A, B coefficients used in the Discrete Vortex Method calculations
A],, B2 half-coefficients used in the Discrete Vortex Method

2 2
calculations
radiation roll damping coefficient

B1
critical roll damping coefficient
Bc
element of matrix of damping coefficients
Bjk
vortex roll damping coefficient
Bv
drag coefficient
CD
vortex force coefficient
CFv
inertia coefficient
Cm
pressure coefficient
CP
roll damping moment coefficient (see Equation 5.5.1)
^rdml
total vortex roll damping moment coefficient, obtained by
^rdm2
summing c ^ g over the ên9th of the ship
complex number used in the determination of L

Cz z
complex number used in the determination of L
c? e
F aerodynamic force
inertia force
Fm
vortex force
Fv
Keulegan-Carpenter number ( = U^T / d)
Kc
flow length scale in the physical plane
Lz
flow length scale in the transform plane
Lc
M transformation constant
vortex roll damping moment

Mv
Re Reynolds number
St Stokes number ( = Re / K )
c
12
T period of oscillation of flow
U,V freestream flow velocities in the transform plane
UQ freestream velocity amplitude in the physical plane (used
in oscillatory flow)
freestream velocity in the physical plane
V-j- transpiration velocity
Vz flow velocity scale in the physical plane
W complex potential
a bilge keel span

i. L.
a^ itn complex ordinate in the transform plane

*hh
a.j* itn complex ordinate in the physical plane
bv sectional vortex roll dampingcoefficient (see Equation
3.9.10)
crdm2 sectional vortex roll damping moment coefficient (see
Equation 5.5.2)
d cylinder diameter
i /-I
1 moment arm
a
ls length of ship section
lv measure of vortex size
p pressure
r bilge radius
t time
t‘ time since the beginning of a particular time cycle
u,v local velocities in the real and imaginary directions
x,y real and imaginary components of complex ordinate

13
z complex physical plane coordinate
zR height of roll centre above ship base
r L.
a. itn transformation angle in the physical plane
a first moment of the vorticity field

th
strength of m vortex
r. . circulation at mesh point (i,j)

T9 J
6 interior edge angle (degrees)
C complex transform plane coordinate

J L
nk kin mode of barge motion ( is roll amplitude)
X representation of edge angle [=(360-6)/180]
y coefficient of viscosity
v kinematic viscosity
p density
a source strength
x shear stress
$ real part of complex potential
\fj. . stream function at mesh point (i,j)

19 J
Q. . vorticity at mesh point (i,j)
19 J
O) roll angular frequency
14
CHAPTER 1
INTRODUCTION
PREFACE
Accurate methods of assessing and predicting the roll damping of ships
are becoming increasingly important for their effective design and
efficient operation. Determination of roll damping has largely so far
been conducted by experimentation. A useful by-product of such
experiments has been the development of empirically or semi-empirically
based models of roll damping. Now, however, with the advent of fast
digital computers it is possible to develop wholly theoretical models
for many types of flow processes. In particular, recent advances in
computational 'vortex shedding methods have begun to enable the
theoretical prediction of ship eddy-making roll damping. Solutions
produced by such computational models would provide a cheap and
effective alternative to model testing. The work presented in this
thesis concerns firstly the development and use of one such vortex
shedding computer model, the Discrete Vortex Method, for the prediction
of the roll damping due to vortex shedding of flat-bottomed barges with
various bilge geometries. Secondly, the model has been extended to
predict the roll damping of vessels of quite general cross-section
subject to certain restrictions. Finally, a more advanced computer model
of the full two-dimensional unsteady Navier-Stokes equations of motion
has been developed which has enabled the prediction of the roll damping
of flat-bottomed barges with rounded bilges to be carried out. In
addition to roll damping data, the methods have been used to provide
15
valuable information on other aspects of the vortex shedding processes.
Much of the work presented here concerning the Discrete Vortex Method is
also presented in Cozens, Bearman and Graham (1986).
1.1 The Prediction of the Motions of a Vessel in a Seaway
1.1.1 The prediction of vessel response in waves using strip theory
The advent of large tankers and high-speed dry-cargo ships has brought a
greater awareness of the importance of the correct prediction of the
various vessel responses to wave loading. Since the pioneering work of
Korvin-Kroukovsky and Jacobs (1957) there have been major advances in
prediction methods. The basis of this early work is the utilisation of
strip theory, where the three-dimensional problem is reduced to a series
of two-dimensional problems by the division of the ship into a number of
transverse strips. The response of each section to wave loading is
calculated, taking into account the damping due to wave radiation.
Notable improvements have been devised by Smith and Salvesen (1970) who
used close-fit methods of geometry representation to predict head-seas
heave and pitch motions quite accurately even for high-speed hulls with
large bulbous bows. Salvesen, Tuck and Faltinsen (1970) employed Frank's
close-fit method (Frank, 1967) in the development of a theory for the
prediction of heave, pitch, sway, roll and yaw motions as well as
wave-induced vertical and horizontal shear forces, bending moments and
torsional moments for a ship advancing at constant speed with arbitrary
heading in regular waves. Comparisons with experimental data, where
available, showed satisfactory agreement.

16
Despite the major advances outlined above, the more lightly damped
motion responses, and particularly roll near resonance, are often poorly
described by conventional strip theory, mainly due to the neglect of the
effects of viscous damping. Salvesen, Tuck and Faltinsen, for instance,
reported that roll response could only be computed with reasonable
accuracy in near-resonance conditions if a viscous roll damping term
were included in addition to the radiation damping. The same also holds
true for the prediction of roll by the strip theory of Ogilvie and Tuck
(1969) which was based on slender body theory. Since the roll motion may
typically have a natural frequency within the range of the wave
spectrum, it is therefore sensitive to the magnitude of the roll
damping.
Many versions of strip theory are capable of predicting vessel response
at forward speed. Although there is great practical interest in the
response of vessels at forward speed in waves, the work on roll damping
in this thesis is restricted to zero forward speed. However, it would
certainly be possible to make empirical corrections to the results
presented here to take account of the effects of forward speed. This
might require an appropriate correction to roll damping at every strip
constituting the vessel.
1.1.2 The roll response of a vessel at zero forward speed
The roll response of vessels at zero forward speed may be affected by a
number of factors. Potential strip theory predicts a damping due to wave

17
radiation which has a linear dependency upon roll amplitude. Such a
theory, however, fails to take account of the non-linearities in the
vessel response. The possible sources of non-linearity at zero speed are
skin friction and eddy-making damping, which give rise to non-linear
damping moments, and non-linearities in the roll restoring moment.
As far as the roll restoring moment is concerned, Denise (1983) suggests
that for shallow-bottomed vessels such as transportation barges linear
potential theory may be invalid because of "non-linearities in the
external forces, other than damping, due to water-structure interaction
in the splash zone". Robinson and Stoddart (1986) have confirmed that
for transportation barges at moderate to high roll angles there are
indeed non-linearities in the restoring moment which may have an
important influence on the vessel response. For deeper draught vessels,
however, such non-linearities are not believed to be generally
significant.
Kato (1958) concludes that the effect of skin friction stress is
negligible for most practical cases. One should expect it to be
relatively more important at model scale than full scale. Patel and
Brown (1985) have estimated that the skin friction damping for a model
barge of 0.8m beam affects the roll response at resonance by only 1%, so
that in the case of a barge even at model scale skin friction stress may
be considered negligible.
Usually it is argued that the most significant source of non-linearity

18
in the damping is the force associated with the shedding of vortices
from the vessel hull and/or its appendages. The importance of vortex
shedding in the roll damping of vessels was first noted in connection
with bilge keels, which are a strong source of vortex shedding. Froude
(1874) hypothesised a non-linear damping term, based on experiments with
a flat plate oscillating in water. White (1895) carried out experiments
to determine the additional damping produced by fitting bilge keels.
Gawn (1940) and Dalzell (1978) performed and analysed tests on
battleships and found that there was a quadratic or cubic variation of
damping with roll amplitude, indicating substantial non-linearities in
the damping term.
Ensuing from earlier experimental work, Kato (1966) developed a
semi-empirical formula for the determination of the roll damping of
vessels due to the effect of bilge keels. Ikeda et al. (1978) have
developed a semi-empirical prediction method for the roll damping terms
due to friction, wave, eddy and lift components for ordinary ship forms.
Agreement with experiment for the simple shapes tested was fair except
for the wide beam/draft variation. The eddy-making prediction method was
based on experiments with twelve two-dimensional sections representative
of container and cargo ships. A near-quadratic damping was found, that
is, damping was roughly proportional to the squares of both roll
amplitude and frequency. However, unappended sections with sharply
defined flow separation points were not tested. Since the method was
empirically based, its validity must presumably be restricted to vessels
of largely similar cross-section to those used in the development of the
method.
19
1.2 The Discrete Vortex Method
In response to a need for more accurate and more generally applicable
methods for the prediction of the eddy-making component of ship roll
damping, researchers have turned to develop wholly theoretical models
for its prediction. Research has been directed, until the later stages
of the work reported here, towards the Discrete Vortex Method.
1.2.1 Early development of the Discrete Vortex Method
Many detailed reviews of the historical development of the Discrete
Vortex Method (DVM) have been published - see, for example, Clements and
Maul! (1975), Leonard (1980) or Jaroch (1986).
Rosenhead (1931) was the first researcher to introduce the concept of
the discrete vortex representation of vortex sheets. He modelled
successfully the Kelvin-Helmholtz instability of a two-dimensional
infinite perturbed vortex sheet by replacing one wavelength of the
vortex sheet with a distribution of evenly spaced discrete potential
vortex elements along the sheet, thereby ignoring the effects of
diffusion at finite Reynolds number. However, subsequent work by Hama
and Burke (1960) repeating Rosenhead's work with a larger number of
vortices and smaller time steps showed that whilst the vortices did
indeed form into clusters, as in Rosenhead's original calculations,
their paths were extremely irregular and occasionally crossed, an
occurrence which is physically impossible for vortex sheets. Hama and

20
Burke discovered an error in Rosenhead's calculations, which, when
rectified, improved the stability of their calculations. Rosenhead's
work was extended by Abernathy and Kronauer (1962) who examined the
stability of two initially parallel shear layers using the discrete
vortex model. It was found that the two layers interacted to form the
familiar von Karman vortex street pattern.
The discrete vortex approximation was first applied to the sheet shed
from an elliptically loaded wing by Westwater (1935). The rolling up
process followed closely that predicted analytically by Kaden (1931). A
re-examination of Westwater's work by Moore (1974) revealed difficulties
which were similar to those discovered with Rosenhead's model especially
once many point vortices and much smaller time steps were employed. Most
reporters on the work of Westwater have not mentioned the somewhat
irregular appearance of the vortex sheet even in his calculations. This
irregularity points to the fact that even with large time steps and few
vortices instabilities in the vortex sheet still can occur. Moore
concluded that it was the method of discretising the vortex sheet into
point vortices which caused the instabilities. At certain times during a
typical calculation the discrete point vortices approached each other
closely; a spuriously large interaction then sometimes occurred which
disrupted the calculation. Moore achieved suppression of the
instabilities by successively amalgamating the innermost vortices of the
spiral. Chorin and Bernard (1973) chose an alternative approach for
suppressing the instabilities. They gave each point vortex a cut-off
radius, so that the velocity at the centre of the vortex was zero rather
than infinite as was the case for the pure potential vortex. The cut-off
21
thus prevented spuriously large interaction between pairs of vortices.
Clements (1973) applied the Discrete Vortex Method to sharp-edged bluff
bodies. His model of flow behind a square based block modelled
separation by the introduction of potential vortices into the flow at
the sharp shedding edges, their initial velocity and strength depending
only on the flow velocity at the "outside edge of the boundary layer".
Clements' model predicted the form of the shedding pattern and the
Strouhal number to within reasonable agreement. The model was
successfully extended in Clements and Maull (1975), where vortices were
placed near the separation points, their strength now being determined
by application of the Kutta Condition, which was used as the "inviscid"
analogue of boundary layer separation.
Kuwahara (1973) and Sarpkaya (1975) performed calculations for an
inclined flat plate, transforming the plate into a circle using a
Joukowsky transformation. Sarpkaya estimated the strength of the
nascent, newly introduced vortices from the equation:
dr/dt = i u % h 1.2.1
where the velocity Us^ was taken as the mean velocity of the latest four
vortices introduced into the shear layer. But Kiya and Arie (1977) have
argued convincingly that Sarpkaya should not have based his formula on
the velocity "U but rather "Umax", the velocity at the outer edge of
the shear layer. The position of the vortices was estimated from
satisfaction of the Kutta Condition. The generalised unsteady Blasius
theorem (Milne-Thomson, 1968) was employed for the derivation of the
normal force, which was predicted 20 - 25% higher than in experiment.

22
Strouhal number was predicted satisfactorily.
Further computations of the separating flow past flat plates at
different angles of attack and different freestream conditions have been
carried out by Kiya and Arie (1977), with separation positions for the
nascent discrete vortices fixed at a distance "as" away from the edge of
the flat plate. They carried out a careful analysis of the optimum value
of ag . Kiya, Arie and Harigane (1979) calculated the flow past a nearly
normal flat plate using vortex decay to obtain results in good accord
with experiment. Kamemoto and Bearman (1978) carried out an extensive
investigation into the parameter a$ and suggested a suitable range of
values within which predicted quantities did not vary substantially.
Even bluff body flows where no geometrically fixed separation point
exists have been modelled using the Discrete Vortex Method. In his
pioneering work on the application of the DVM to bluff body flows,
Gerrard (1967) modelled the shedding of vortices behind a circular
cylinder in uniform flow. Rather than placing the nascent vortices near
the flow separation points, Gerrard placed them roughly one cylinder
radius out from the centre line of the flow and one radius downstream of
the centre of the cylinder. Calculated lift and drag curves showed time
dependent oscillations, which could in many cases be interpreted to
yield a representative Strouhal number. Further discrete vortex
representations of the flow past circular cylinders have been carried
out by Sarpkaya (1968), who simulated the early roll up of the shear
layers. Stansby (1977) modelled with some success the flow past a
circular cylinder both by fixing the separation points at an average

23
position estimated from experiment and also by predicting the separation
points by "using the condition that the separating streamlines are
tangential to the surface". In reality the actual shedding positions are
unsteady and Reynolds number dependent. Thus fixing the separation
points at an average position estimated from experimental evidence leads
to certain inaccuracies; no convincing argument can be found for
determining them by the tangential streamline method. On the other hand,
these simplifying assumptions did yield useful results. Stansby also
calculated oscillatory flow past the cylinder by fixing the separation
points.
In bluff body flows where no geometrically fixed separation point
exists, discrete vortex models are most often used either, as described
above, by fixing the separation points in advance from, for example, the
results of experimental analysis, or else by interaction with some form
of boundary layer analysis. Deffenbaugh and Marshall (1976) matched an
outer discrete vortex model with an inner viscous flow region,
calculating the separation points by a solution of the unsteady laminar
boundary layer equations. Pleasing agreement with other work was
obtained in the early stages of the flow, but at later times the
separation points moved unrealistically far forward, adversely affecting
lift. More successful for long times were the calculations of Sarpkaya
and Schoaff (1979) who employed a heuristic model of vortex decay to
reduce the levels of circulation in the vortex sheet. It was their
observation, as it is of many workers with the Discrete Vortex Method
(see, for example, Clements and Maull, 1975), that the concentrated
vortices of the wake contained typically 30% more vorticity than would
24
be found in experiment. Vortex decay was employed as a means of reducing
the circulation strength in the wake nearer to experimental values.
Sarpkaya and Schoaff also utilised a method of rediscretisation of the
vortex sheets which had been pioneered by Fink and Soh (1974). Fink and
Soh demonstrated that a source of instability in the roll-up of vortex
sheets which had been noticed, for example, by Hama and Burke (1960),
could be removed by rediscretisation of the vorticity. Rediscretisation
involves the placing of each vortex at each time step at the midpoint of
its segment. One might view this as another method of smoothing out the
instabilities referred to above caused by the discretisation inherent in
the potential vortex model.
1.2.2 The application of the Cloud-in-Cell technique to the Discrete
Vortex Method
One of the main limitations of the Discrete Vortex Method as described
so far is that it depends on the calculation of the influence on each
point vortex of the rest of the point vortices as determined by the
Biot-Savart law. For a computation with N vortices this procedure
requires 0(N2 ) operations to compute the required velocities. This
effectively limits the maximum number of vortices in a typical
computation on a mainframe computer to hundreds rather than thousands.
The Cloud-in-Cell (or Vortex-in-Cell) technique which is described in
Chapter 6 retains the Lagrangian treatment of the vortex field with its
attendant advantages in terms of its ability to represent without
smearing (or without introducing the levels of artificial viscosity
associated with conventional grid methods) small features of the flow

25
moving with the fluid through a large region, but solves the Poisson
equation for the velocity field on a fixed Eulerian Mesh. The technique
evolved out of methods such as the Particle-in-Cel1 method developed at
Los Alamos in 1955. A comprehensive review of that method may be found
in Harlow (1964). A development of the PIC method is the FIuid-in-Cel 1
method of Gentry, Martin and Daly ( 1966). The advantage of the CIC
hybrid Lagrangian/Eulerian approach is the greatly reduced computational
cost for small loss in accuracy. For a mesh with M intersections the
number of operations is 0(N ) + 0(Mlog2M). Reasonably sized calculations
might be speeded up by two orders of magnitude by use of the CIC method.
Christiansen (1973) was the first to report the use of the CIC technique
in studying the motion of a two-dimensional incompressible, inviscid and
homogeneous fluid. He computed successfully such fundamental cases as
the Kelvin-Helmholtz instability and the formation of the Karman vortex
street from two parallel shear layers, echoing the work of Rosenhead
(1931) and Abernathy and Kronauer (1962) respectively. Baker (1979)
applied the method to the evolution of the vortex sheet from an
elliptically loaded wing and a wing with a flap deployed. Naylor (1982)
has employed the method for uniform flow past a normal flat plate.
Basuki (1983) used it to calculate unsteady flow over aerofoils and
cascades.
In the present work the Cloud-in-Cell method has not been used in
conjunction with the Discrete Vortex Method. This is because alternative
procedures for reducing the computational expense were found. The
importance of the CIC method here lies in its use in the solution of the
full Navier-Stokes equations. This will be discussed in the next
section.
26
1.3 The Solution of the Full Two-Dimensional Unsteady Navier-Stokes
Equations
There are two main limiting factors in the application of the Discrete
Vortex Method to the computation of vortex shedding from ship hulls.
1) The Discrete Vortex Method, with its assumption of an infinitely thin
vortex sheet, is incapable of representing the effects of either
diffusion or turbulence in the vortex sheet.
2) The Discrete Vortex Method cannot predict the flow separation point.
This is not a problem in cases where the separation point is fixed by
geometry. Shedding off a flat plate or bilge keel are typical examples
where it can be used with a high degree of success. On the other hand,
rounding of the shedding surface, such as on a circular cylinder or
barge with rounded bilges, poses many problems as to the location of the
shedding point. One solution which is feasible for small amounts of
rounding or more generally where the separation point is known not to
move markedly is to fix the point either from the results of experiments
or by a heuristic approach. An alternative would be to carry out a
boundary layer calculation to determine the point of separation. An
obvious weakness of this is that boundary layer methods fail at and very
near the separation position so that it cannot be predicted with much
certainty.
In many cases, then, it is clear that some model of the vortex diffusion
process, not only in the wake but also near the separation point, is
27
required. Such models are described in the succeeding sections, and
involve solution of the unsteady two-dimensional Navier-Stokes
equations.
1.3.1 Eulerian solutions of the Navier-Stokes equations
A useful survey of solutions to the Navier-Stokes equations may be found
in Roache (1972). When an Eulerian mesh method is employed the
Navier-Stokes equations are often solved in their vorticity transport
form (see Chapter 6), although solutions of the "primitive"
Navier-Stokes equations (the so-called "u,v,p" solutions) are also
found. The vorticity transport equation is parabolic in form. Since the
computational expense of solving parabolic equations is often large much
effort has been expended on producing fast and efficient solution
schemes. As early as 1957, Richtmyer (1957) had presented more than ten
numerical schemes for the solution of parabolic equations. Solution
schemes for these equations can be divided into two categories, explicit
and implicit. Explicit schemes advance to a new time step using only
variable values determined at the old time step; implicit methods
advance using variable values from the new time step as well. In general
explicit methods are simpler to program but computationally more
expensive since they are limited by stability arguments to short time
steps; implicit methods are more complicated but generally
computationally cheaper since time steps can often - at least in
principle - be arbitrarily large. One important feature of all solutions
of the Navier-Stokes equations of this type is the requirement of small

28
mesh size near the body side in order to enable adequate representation
of the boundary layer. This has important implications on computational
expense since the time step size will usually be governed for stability
or accuracy reasons by the smallest mesh cell size.
A frequently referenced example of computations carried out with an
explicit scheme is the work of Thoman and Szewczyk (1969) who used a
first order explicit scheme with "upwind" differencing to improve the
stability of the calculations. Their computations, for uniform flow past
a circular cylinder, were remarkable in that a Reynolds number of
nearly one million was achieved. Accuracy of the solution was poor,
however, because the upwind differencing scheme introduced "artificial"
viscosity into the computations, which had the effect of damping the
flow and in many cases preventing vortices from detaching from the
cylinder.
For many applications the Alternating Direction Implicit (ADI)
second-order finite difference solution scheme for parabolic equations
due to Peaceman and Rachford (1955) is considered to be the most
suitable. Peaceman and Rachford themselves tested the ADI method by
solving the heat flow equation. By comparison, a typical explicit
procedure required 25 times as much computer time. Borthwick (1986)
employed the ADI method in a comparison with Thoman and Szewczyk's work
on the circular cylinder. The ADI method produced much more realistic
results. Nevertheless, computing restrictions have, until recently at
least, generally limited the ADI method to computations involving
Reynolds numbers of hundreds rather than thousands.

29
1.3.2 Lagrangian solutions of the Navier-Stokes equations
In contrast with the Eulerian solutions to the vorticity transport
equation a number of Lagrangian schemes have been developed which
possess quite different properties. These have generally evolved from
the inviscid Discrete Vortex Method which solves the convection portion
of the vorticity transport equation, but they also have the effects of
diffusion taken suitably into account. Three alternatives for the
modelling of the diffusion term have been investigated. Generally, since
a large number of vortices is required to represent adequately the
effects of diffusion, the Cloud-in-Cell technique has been employed to
reduce computation time.
1) Representation of the shear layer with discrete vortex "blobs" (see
Leonard, 1980). Diffusion can be taken into account by allowing the
blobs to increase in size with time. This approach has been little used.
2) Diffusion onto neighbouring point vortices (see Raviart, 1986). This
requires a large number of particles to ensure that diffusion is
adequately represented. Problems can be envisaged with the imposition of
the wall vorticity boundary condition. No numerical tests or
implementations of the algorithm are known.
3) Addition of a random walk to the discrete vortices at each time step
to represent statistically the effects of viscosity. Chorin (1973)
introduced this method for flow past a circular cylinder. Agreement with
30
experiment was good at Re = 1000, but worse at both higher and lower
Reynolds numbers. At the end of the computations there were 300 vortices
in the flow. Milinazzo and Saffman (1977) found that to achieve
acceptably accurate results for the viscous decay of a uniform vortex
using this scheme the number of vortex elements had to be large, and in
direct proportion to Reynolds number. It is difficult to know how many
vortices are required for a given degree of accuracy, but it seems
likely that Chorin's computations did not contain enough vortices for
adequate representation of either the initial separation of the vortices
or the wake.
Stansby and Dixon (1983) employed the random walk to compute shedding
from a circular cylinder in both uniform and oscillatory flow. In
applying the random walk method, with certain modifications, to various
bluff body flows, Lewis and Porthouse (1983) have produced results which
seem qualitatively plausible, although comparisons with experiment were
limited, and not particularly good, except for the Blasius boundary
layer. Several problems have been encountered in the use of the random
walk method. First, there is an upper Reynolds number limit in terms of
computational cost imposed by the requirement that, for a given degree
of accuracy, the number of vortices must increase as Reynolds number
increases. Secondly, there is a lower Reynolds number limit which is due
to the inadequacy of the random walk model at low Reynolds number.
Thirdly, the random walk model is incapable with present computing power
of realizing the transition to turbulence, which would require the
representation of much smaller length and time scales. In fact, it
appears for the above reasons that the random walk method is applicable
31
in quite a narrow Reynolds number range, although inadequate work has
been carried out on the method to enable this statement to be
quantified.
1.3.3 Comparison of Eulerian and Lagrangian schemes
The main advantage of a Lagrangian scheme is that it is able to resolve
small details of the flow moving with the fluid through a large region.
It is economical for flows with a strong convection of vorticity because
the vortex particles are the computational representation of vorticity.
Areas away from the wake are not represented simply because there is no
such requirement. Large time step size can be used because there are
apparently no stability restrictions either in the flow or with the
boundary conditions. One of the most serious drawbacks of a Lagrangian
scheme is its inability to cope with large distortions or large amounts
of diffusion in the fluid, whereas an Eulerian scheme proceeds without
difficulty in these situations. The obvious disadvantage in particular
with the random walk method is the limitation on Reynolds number and
achievable accuracy imposed by computing restrictions on the number of
vortices in the flow.
1.3.4 Hybrid Lagrangian/Eulerian solutions of the Navier-Stokes
equations
It would appear from the foregoing discussion that Lagrangian vortex
schemes are suited for solving the convection problem, but solution of
the diffusion problem is not as accurate or efficient as with a typical

32
Eulerian scheme. A combination of Lagrangian convection and Eulerian
diffusion schemes could therefore eliminate many of the problems of both
methods. Such a combination was first described by Garder, Peaceman and
Pozzi (1964) for the solution of high Peclet number miscible
displacement equations such as might be encountered in oil recovery
problems. The Peclet number is analogous in these equations to the
Reynolds number in the vorticity transport equations. Price, Cavendish
and Varga (1968) showed that such a Moving Point Method (MPM) avoided
the problems of numerical viscosity and spurious oscillations of the
solution so frequently encountered in pure finite difference solutions,
and was more accurate than those solutions for sufficiently large Peclet
number. They indicated, however, that the MPM was not convergent under
mesh refinement. Farmer (1985a) has demonstrated the convergence of an
improved MPM method and showed for a typical test case that it is much
more accurate than either 1st or 2nd-order Eulerian upwind differencing
schemes. In an excellent review of modern moving point techniques,
Farmer (1985b) describes two other techniques, the Modified Method of
Characteristics and the Hybrid Moving Point Method.
The methods described by Farmer for miscible displacement problems have
obvious application to the solution of the Navier-Stokes vorticity
transport equation. The Moving Point Method especially could be employed
to extend the Discrete Vortex Method in its Cloud-in-Cell formulation to
cope with the effects of diffusion. However, solution of the vorticity
transport equation entails some different problems to the miscible
displacement equations, especially in the treatment of the boundary
conditions. The MPM as adapted to the solution of the transport of

33
vorticity (we shall call this the Hybrid Moving Vortex Diffusive Method,
or HMVDM) has been adopted in this present study. No previous research
in this area can be found. The HMVDM represents a new approach to the
solution of the Navier-Stokes equations.
1.4 The Application of Vortex Methods to Vortex Shedding from Vessel
Hulls
This thesis concerns the application of the Discrete Vortex and Hybrid
Moving Vortex Diffusive Methods described above to the problem of vortex
shedding from rolling vessel hulls. The research work took as its
starting point the earlier work of Downie, Bearman and Graham (1984 and
1985) on the roll damping of barges including the effects of vortex
shedding. Typically, such barges are large ocean-going flat-bottomed
barges which are utilised in the offshore oil industry. These vessels
are usually towed to their destination bearing bulky cargos such as
sections of oil rig structures. During transportation the cargo
structures are secured to the vessel deck with welded steel fastenings.
The fastenings and the cargo itself have to be designed to withstand
inertia forces due particularly to roll in a seaway. Both cargo modules
and fastening designs are strongly dependent on accurate predictions of
barge roll motions.
Downie et al. considered only barges of rectangular cross-section for
the critical design case of a barge with zero forward speed in beam
waves. Such a situation might arise due to engine failure or loss of tow
lines causing the barge to turn beam to waves. The work was an
34
application of the earlier research using the Discrete Vortex Method on
the drag due to vortex shedding for sinusoidally oscillating flow around
a semi-infinite square edge (see Graham, 1980). A matching technique was
derived to match the flow calculated by the Discrete Vortex Method to
the flow round the barge cross-section. Use of such a matching technique
is a powerful tool since it means that discrete vortex calculations for
a particular edge geometry need be performed only once, for all time,
regardless of the overall barge geometry. The success of this technique
led to the development by the present author of a Discrete Vortex Method
for more general barge bilge shapes, particularly the isolated square
edge with rounding and/or bilge keel (see Cozens, Bearman and Graham,
1986). The work is also reported in this thesis. Subsequently the method
was extended to cope with complete ship cross-sections, without the
matching technique. Finally, for a barge with rounded bilges, where the
Discrete Vortex Method was considered to be generally inappropriate, the
Hybrid Moving Vortex Diffusive Method was developed.
Other work on the roll damping of rectangular cross-section barges using
the Discrete Vortex Method has been carried out by Patel and Brown
(1986) following the earlier research of Brown, Eatock Taylor and Patel
(1983) on barge motion responses excluding the effects of vortex
shedding. Patel and Brown were prevented by computational considerations
from carrying out a discrete vortex computation with enough vortices or
for a sufficient duration to enable an accurate answer to be obtained.
Nevertheless, damping predictions were in many cases in fair agreement
with experiment. They also presented damping results for a barge with
rounded bilges, but, though this is not stated explicitly, it is

35
believed that these results are empirically based.
Robinson and Stoddart (1986) have presented results for a rounded bilge
barge which are based on the rectangular cross-section calculations of
Downie et al. (1984) coupled with an empirical correction for the
effects of edge rounding based on the experiments of Ikeda, Himeno and
Tanaka (1977).
Faltinsen and Sortland (1987) have presented results for a rectangular
cross-section barge in pure sway. The effect of bilge radius and bilge
keels was also tested. Pleasing comparisons with experiment were made,
and calculations performed for a selection of the above parameters for
barges of differing aspect ratios. This work is an interesting
application of the Discrete Vortex Method at reasonable computational
cost to shedding computations around complicated shapes.
Few workers have employed the Discrete Vortex Method for the modelling
of vortex shedding from ships with more general cross-sections. Soh and
Fink (1971) employed it for the impulsive start of the flow around a
heaving ship with bilge keels. No damping results were presented. Ikeda
and Himeno (1981) have modelled the sway of various Lewis-form cylinders
representative of ship shapes. Shedding patterns appeared plausible but
again no damping results were presented. They did not attempt to model
roll motions.
Odabasi, Incecik and Edwards (1985) have developed a method which
predicts the roll damping of ships due to vortex shedding via a process
of local matching to the results of Discrete Vortex Method computations

36
for simpler shapes. In order to extend the method to edges with
rounding, calibrations against experiment were carried out to provide
empirical constants which were used in the damping predictions.
Brook (1986) has carried out comparisons on selected vessels between
experimental dampings, Ikeda's roll damping prediction method (Ikeda et
al., 1978) based on experimental observations and the method due to
Odabasi, Incecik and Edwards (1985). The latter method predicted
dampings for naked hull configurations which were in reasonably good
agreement with experiment and which were in some cases significantly
better than the predictions using the Ikeda method. However, the method
tended to overpredict the effect of bilge keels relative both to
experiment and to Ikeda's method for comparisons carried out for the
Fisheries Protection Vessel, the Sulisker.
No research involving the solution of the full Navier-Stokes equations
for the rolling of vessels can be found.
The present work on the modelling of vortex shedding from vessel hulls
represents two major innovative contributions. First, the Discrete
Vortex Method has been developed to predict the vortex shedding from a
variety of barge bilge shapes via a specially modified matching process,
and also from complete general ship cross-sections in roll (or sway).
Secondly, a new computational technique for solving the Navier-Stokes
equations has been developed. Computations have been carried out for
certain cases where bilge rounding would normally prevent Discrete
Vortex Method calculations from being performed successfully.

37
CHAPTER 2
DESCRIPTION OF THE PROCESSES OF VORTEX SHEDDING IN OSCILLATORY FLOW
It is convenient to describe the nature of vortex shedding in
oscillatory flow in this Chapter so that reference can be made later in
the thesis to the various flow patterns which are observed. The
discussion here is related to observations of the vortex shedding from
the bilges of a rolling model barge for which flow visualisation
experiments have been carried out as part of the current research.
A programme of experiments was carried out about four years ago at
Imperial College to visualize the flow round the bilges of a rolling
rectangular cross-section barge. The most important results of these
experiments are presented in Downie, Bearman and Graham (1984). A vortex
shedding pattern was visualized, and this was in reasonable agreement
with theoretical predictions as far as vortex paths and vortex sizes
were concerned. In the current set of experiments the same basic barge
model was roll tested with bilge keels of different spans and bilges of
different radii in order that visualisations could be made of the
fundamental shedding patterns involved and also of the effects of
different combinations of edge rounding and bilge keels.

38
2.1 Experimental Method
2.1.1 Experimental apparatus
The basic barge model is that used in the previous experiments, with the
addition of bilge keels and/or rounded bilges. Details of the model are
shown in Figure 2.1. The model was tested in the wave tank at Imperial
College. This tank has a working section width of 0.61m and depth of lm,
although for the present tests the depth of the water was only 0.6m. For
these tests, the wave-maker was not used; instead, the barge was force
rolled. A schematic representation of the experimental apparatus is
shown in Figure 2.2.
2.1.2 Test procedure
Forced-roll tests were carried out at a variety of combinations of
angular frequency, roll amplitude, bilge keel span and bilge radius. For
the square edge with bilge keel most of the testing was performed at a
roll angular frequency of 4 rad/sec and roll amplitude of 5.7°. Other
tests were carried out at a higher frequency (6 rad/sec) and a higher
amplitude (9.9°). For the rounded edge, the two different bilge radii
were tested with no keels and keels of varying span at an angular
frequency of 1.5 rad/sec and an amplitude of 11.4°. Details of the
figures discussed in this chapter are shown in the table below.

39
Barge beam = 0.36 m, length = 0.61 m and draught = 0.07 m
Fig. no. Bilge rad. (mm) Span (mm) Ang. freq. (rad/s) Amp. (deg
2.3 13 13 1.5 11.4
2.4 0 0 4.0 9.9
2.5 0 25 6.0 5.7
2.6 0 76 4.0 5.7
2.7 25 0 1.5 11.4
2.8 13 0 1.5 11.4
2.9 13 6 1.5 11.4
2.10 13 13 1.5 11.4
2.11 0 76 4.0 5.7

40
In Chapter 7 a formulation for a characteristic value of the Stokes
number is given. For the tests reported here Stokes number on that basis
is of the order of 5000. A typical equivalent full scale Stokes number
would be of the order of 5 million. It is not completely clear what
differences in roll damping there would be due to Stokes number between
model scale and full scale. If the main component of damping is skin
friction then the differences would be large; otherwise the major
unknown factor is the effect of Stokes number on the separation point.
2.1.3 Flow visualization techniques
The flow was visualized using two separate techniques - polystyrene
beads and dye.
1) Polystyrene beads
Very small polystyrene beads with specific gravity near that of water
were injected in solution into the water at low speed. A thin sheet of
intense light in a plane parallel to that of the barge end-plate was
shone towards the source of the beads. This had the effect of
illuminating beads in the shedding plane. This particular method of
visualization is best in situations where the vortex shedding is strong
and where significant detail is not required. It is capable of yielding
a good picture of overall flow patterns.

41
2) Dye
A concentrated solution of dye was injected into the flow through fine
hypodermic tubing. Illumination was from the further side of the barge.
This method of visualization is best for lower local flow speeds. In
this case the flow processes near the shedding edge can be observed in
great detail. Care was taken in the experimentation to avoid possible
interference between the tubing admitting the dye into the flow and the
flow patterns themselves. Dye was ejected into the flow from points
near, but not at, the shedding edge. Concentrated dye was used which
permitted very low ejection speeds.
2.2 Photographs of the Vortex Shedding Processes
Photographs were taken of the vortex shedding as part of the flow
visualization process. The photographs were used for two main purposes -
to gain a better qualitative understanding of the flow pattern, and to
measure vortex sizes to enable an assessment of the quality of the
theoretical vortex predictions. In this section a representative sample
of the photographs is shown and their qualitative importance briefly
analysed.
Figure 2.3 shows a close-up photo of a typical single vortex shed from a
keel. Apart from the main vortex "spiral" and its central "core", it can
be seen that there are a number of smaller, secondary vortices forming
on the main vortex sheet.

42
Figure 2.4 shows a series of photographs representing one complete time
cycle of flow round a square edge. Flow visualization was performed
using polystyrene beads. At the beginning of the time cycle, with the
barge at peak amplitude, a large vortex may be seen to the right of the
barge. Its pair is just beginning to form at the edge on the underside
of the barge. By photo 2 this nascent vortex has grown considerably,
whilst the first vortex of the pair has been swept right round to a
position below and to the right; by photo 3 it is below and to the left,
indicating its speed of travel. At the same time in photo 3 (peak
amplitude at the opposite extreme of roll) the second vortex of the pair
has now ceased growing, and there is a new nascent vortex at the right
of the shedding edge. It appears, in general, that vortex pairs are of
roughly equal, but opposite, strength. Since for the square edge
equalisation of vorticity does not necessarily occur at the end of the
growth of the second vortex, the pairing process splits the second
vortex sheet leaving a smaller amount of residual vorticity which may be
seen below and to the right in photo 4 and to the right in photo 5 of
the nascent vortex. This is soon engulfed by the nascent vortex. Photo 4
shows the nascent vortex having grown rapidly; the structure of the
previous pair is now indistinct, although its effect is clear enough.
Photos 5 and 6 show the further growth of the nascent vortex.
Figure 2.5 shows a set of photographs representing one time cycle. The
figure refers to a bilge keel of 25mm span at the higher roll rate; the
visualization was performed using dye. Photo 1 shows the situation at
peak amplitude. One vortex of each pair has already been cast off, and
43
the second vortex can be seen just starting to grow. It grows through
photo 2, until by photo 3 the elements of the pairs are almost evenly
matched in strength. Notice how the structure of the right hand pair
appears to be much better preserved than that of the left hand, although
this effect could possibly be due to some asymmetry in the ejection of
dye into the flow. In photo 4 (peak amplitude) both pairs have left the
edge completely and there are signs of a new nascent vortex at each
edge. These grow considerably through photos 5 and 6, until by photo 1
again peak amplitude has been reached (and the end of the time cycle)
and the second vortex of each pair begins to be shed. Figure 2.6 is a
bead visualization of flow round the 76mm keel for a complete time
cycle. All the events described above for Figure 2.5 can be seen
happening here in considerable detail, with the barge moving from one
peak of amplitude (photo 1), through to the other (photo 6) and then
back again. Shedding off each edge is 180 degrees out of phase with the
other. Vortex pairs convect nearly straight downwards, at an angle of
about 45 degrees to the plane of the bilge keels. Figures 2.5 and 2.6
show that the presence of a bilge keel makes the vortex structure more
stable in comparison with the square edge (Figure 2.4). In the two cases
shown, there is little or no residual vorticity at the end of the
pairing process. The pairs convect away from the vessel surface, whereas
for the square edge convection occurred along the barge surface. There
is a bistability in the shedding process. Vortex pairs may either
convect upwards, towards the free surface, or downwards. In general
(see, for example, Figure 2.6) in the cases where the vortex pairs would
have been strong enough to have been affected by the free surface, they
have convected downwards away from it.

44
Figure 2.7 shows the flow pattern for an edge with rounding but no keel.
A "double" shedding pattern is observed, in that two strong vortex pairs
are shed per cycle. The two pairs can be seen about to convect away from
the edge in photos 1 and 5. In the presence of rounding, the shedding
appears weaker than for the square edge; the vortex pairs remain very
near the barge surface as they convect away. The double shedding pattern
is a result of the residual vorticity noted in Figure 2.4 being large
enough to force a second strong pairing of vortices. In fact the
residual vorticity has become as large as the "primary", thus setting up
two shedding patterns apparently equal in strength. Figure 2.8 shows an
intermediate stage where residual vorticity is greater than for the
square edge, but not as great as for the larger bilge radius. A double
shedding pattern per cycle is observed, with pairs forming at photos 1/2
and 5/6, but the first pairing, to the left, is stronger than the
second.
In Figure 2.9, where the keel span is half the bilge radius, the
"single" shedding pattern observed in earlier figures is regained,
despite the relatively small size of the keel. However, small secondary
disturbances can still be seen. Figure 2.10 shows the flow situation
with a keel span twice as great. The general appearance of the flow is
now much like that which would be expected from an edge without any
rounding. The figure shows particularly clearly the various stages in
the shedding process. In both Figures 2.9 and 2.10 a small amount of
residual vorticity is apparent at the end of the main pairing process.
This pairs up with a portion of the new nascent vortex to form a weaker
secondary pair. This is clearly seen in photos 2-5 of Figure 2.10.

45
Figure 2.11 shows the same flow situation as that in Figure 2.6. Two
separate dye sources have been used here, so that the degree of
two-dimensionality of the shedding can be assessed. Slight differences
can be observed with the far vortex appearing to have convected further
from the edge than the nearer one. This example typifies the
observations made for all the different bilge configurations. Only very
occasionally could there be noted any significant difference in, say,
the convection paths.

46
CHAPTER 3
THE DISCRETE VORTEX METHOD
The application of the DiscreteVortex Method to the flow around
isolated edges has been described in detail by Graham (1980), who drew
on the earlier work of Pullin (1978). This chapter describes its
development in the first place for more complicated isolated edge
geometries, particularly the edge with rounding and/or bilge keel.
Vortex shedding from such shapes can be directly related to vortex
shedding from the bilges of barges in roll or sway by a matching
process. Section 3.9 deals with the extension of the method to cope with
the vortex shedding in roll ofgeneral symmetric ship hull
cross-sections with an arbitrary separation position and arbitrary roll
centre.
3.1 Transformation
Two alternative methods of representing the influence of the body were
considered. Firstly, it is possible to represent the body by a
distribution either of sources and sinks or vortices. A good example of
the use of the Discrete Vortex Method with such a singularity
distribution may be found in Sakata, Adachi, Saito and Inamuro (1983).
The second alternative is the use of a transformation technique. An
important advantage of the transformation technique is usually higher
accuracy especially for geometries with sharp edges such as will often
be encountered in vortex shedding problems; one drawback is the fact

47
that transformations are more difficult to implement for complicated
shapes.
Because of the possibility of greater accuracy of representation, for
the current research a numerical Schwarz-Christoffel transformation was
employed which enabled isolated edges of arbitrary geometry to be
transformed to a straight line segment as shown in Figure 3.1. The
general transformation can be written as:
n a-j/ir
dz
= M n U - an.) 3.1.1
He i=l
where the symbols are as defined in the figure. M is called the
transformation constant and defines the scaling relationship between the
physical and transform planes. The normal sign convention for angles is
adopted, that is, angles are taken to be positive for an anticlockwise
rotation when the body is circled in an anticlockwise sense. The full
numerical integration of Equation 3.1.1 which is necessary when a large
number of vertices is involved has been discussed in detail by Davis
(1983), and will be treated in Section 3.1.3. It is noted first,
however, that for simple edge geometries versions of Equation 3.1.1 can
be derived which yield analytical solutions of the integral. Two such
geometrical configurations are discussed below.
3.1.1 Plain isolated edge transformation
Such a geometry with its transformation is shown in Figure 3.2. For this
problem, with only one vertex, this being at the origin, Equation 3.1.1
48
reduces to:
3.1.2
where:
x ~ 2 “ T M 1 + “l ^ 3.1.3
Integration of Equation 3.1.2 and setting the constant of integration to
zero leads to a familiar simple form of the Schwarz-Christoffel
transformation:
The second derivative of the transformation, which is required for the
determination of Routh's correction (Clements, 1973), is:
3.1.5
dC2
3.1.2 Plain isolated edge with bilge keel transformation
Such a geometry with its transformation is shown in Figure 3.3. For this
problem, with three vertices in total, Equation 3.1.1 reduces to:
^ -1
j}§ = Mc(C2 + b)7 3.1.6
Integration of Equation 3.1.6 and setting the constant of integration to
zero leads to the following form of the Schwarz-Christoffel
transformation:
z 3.1.7
49
The second derivative of the transformation is:
a _2
= M{c2 + b)? [(C2+b) + (A-2)?2] 3.1.8
d?2
3.1.3 Full numerical solution of the Schwarz-Christoffel transformation
In cases where many vertices are present, such as in the representation
of rounding for the isolated edge case, analytical solutions are not
available and instead a computational solution must be sought. The
problem is centred on the determination of the a^'s which are initially
unknown.
Solution for the a^'s proceeds as follows. Without loss of generality
one of the vertices in the z-plane can be set at the origin, as can one
of the vertices in the transform plane. Further, the first and last
vertices can be set to +i and -i in the transform plane. In order to
determine the remaining values of a^, Equation 3.1.1 must be integrated
to give z in terms of C. In this work an extension of the midpoint
integration rule has been used to determine z ,, from z , where
Az = z - z is small. This requires the division of each facet into

m m+1 m n
sub-elements, the number of the sub-elements required depending on the
facet size and the nature of the transformation in the region.
In order to obtain the unknown values of a. and M, the following
iterative procedure was adopted.

50
1) Guess initial values for the unknown constants. M is complex and can
be conveniently set to (1 + 0 i); the a.'s are chosen such that
ai-l < a i < a i+l anc* t*ie ipterva^ between each a., is equal.
2) Integrate Equation 3.1.1 using the midpoint integration rule between
appropriate limits as follows:
n (?-ai)(“ i/lr)+1 ’m+1

Az
m
= MAC
m
n (ot^/TTj+l
3.1.9
i=l ACm
In the physical plane the vertices will then not be in their correct
locations. Therefore adjust M so that the first and last vertices lie in
their correct position.
3) Use the adjusted value of M and the calculated positions of the
physical plane vertices to provide better estimates of the a^'s.
4) Return to Stage (2) and perform a new iteration assuming the new
values of M and a^. Repeat the iteration procedure until convergence has
been achieved.
More detail on the above process may be found in Davis (1983).
Once converged values of a.. and M have been obtained, then the various
properties of the transformation can be evaluated.

51
First derivative of the transformation. For geometries where symmetry
is involved Equation 3.1.1 can be rewritten in a more efficient form:
j§ = M U - a r)a r/lT V (c2-a?)ai/7T 31K 1 0
where r = (n+l)/2, n odd.
Second derivative of the transformation. Little computational
simplification can be achieved by using symmetry arguments in the case
of the second derivative of the transformation, which is accordingly
derived as:
"n a ./tT n (“ i/ir)

d 2z 3.1.11
= M n (?-a.)
d£2 i=l 1
Determination of z from C. Equation 3.1.9 should be employed to
determine z for given C, using an appropriate number of integration
sub-elements.
Determination of c from z. An iterative scheme was developed for the
determination of c for given z. The scheme involved the calculation of
successive approximations to c using a "shooting" scheme. In the first
part of the scheme, C q s the first estimate of the required value of C,
was set to unity. In the second stage, the corresponding value of the
i 61 d
physical plane ordinate, z ^ = r0 ]d e » was computed by the method
of the previous sub-section. In the third stage, an improved estimate of

52
C, £pew» was provided from the following equation:
3.1.12
rold
where z = re is the actual physical plane ordinate corresponding to c
and R is a relaxation factor. Optimum performance in most conditions was
achieved by under-relaxing the solution by a factor of one half, i.e.
putting R=|. Iteration between stages two and three was then carried
out, with C , , successively being equated to C , until the difference
between z and z ^ was less than the required tolerance.
Flow length scale. Flows around semi-infinite plain (i.e. single vertex)
isolated edges possess no physical length from which a length scaling
can be determined. Further, because of the semi-infinite nature of the
body, there exists no freestream velocity except in certain special
cases. Pullin (1978) derived a formulation of length scaling for plain
isolated edges in uniform flow, this scaling being based on properties
of the transformation, the "freestream" flow velocity in the transform
plane, and the time scale of the flow (i.e. the time from the start of
the flow). Self-similar solutions were thus provided for plain isolated
edges of specified edge angle. Graham (1980) extended the concept to
oscillatory flow, where the time scaling was now provided by the period,
T, of the oscillations. In the present research a generalised expression
for the length scale (l_z ) of the flow in the physical plane around
isolated edges of arbitrary edge geometry has been derived from the
velocity scale (Vz ) in the physical plane as follows. A representative
edge velocity in the physical plane will scale on dW/dz. We choose the
53
constant of proportionality in the scaling to be X, where X here
represents the overall edge angle, so that:
dW
vz = x <37
3 . 1.13
where W is the complex potential in the absence of vortex shedding, and
dW/dz is taken as evaluated at Cz « Cz is a complex ordinate of magnitude
Lz and orientation such that it lies on the edge bisector. We can
further define a length scale in the transform plane, L^, where
L = |C |, and is the point in the transform plane corresponding to
Cz in the physical plane. Thus:
dW de
cZ = xQZf . T = X-
dz
3.1.14
SC
Then we can write l_z as:
XV T
a - / it 3.1.15
Lz = IC;
M n (C -a.) 1
i=l 5 1
In the above equation V is treated as the transform plane velocity
amplitude for oscillatory flow. This equation was solved by iteration
between Cz and C^. Since its solution is computationally expensive, for
isolated edge cases where, for instance, keel span, a, or bilge radius,
r, are small compared to the flow length scale an approximation to this
equation can be solved instead as follows. Far away from the edge, i.e.
where |c| » a^ for all a., Equation 3.1.1 can be rewritten as:
. I (Ot./TT)
dz m i-i 1 m ^X-1
= M c 1-1 = Me
3.1.16
54
Therefore, using Equations 3.1.14 and 3.1.15:
fN0-l/(2A--|) X / (2A-1)
L
z rl (VT) 3.1.17
Equation 3.1.17 is actually the expression used for the flow length
scale for oscillatory flow round a plain isolated edge (see Graham,
1980).
3.2 Basic Concepts of the Discrete Vortex Method
The basis of the DVM, as applied here to vortex shedding off bluff
bodies, consists in the representation of the shear layer from a bluff
body by a series of point potential vortices (see Figure 3.4). For such
a distribution of vortices and their images there exists a complex
potential (W) which may be written for the transform plane as:
W(?) = iV? + J - I rm l09^-Sm> logU+s*) 3.2.1

m
which gives:
3W(C) _ iV + i
2tt I m c-Cm 2 tt £I rr £ + C
L *m 3.2.2
35 m m m
The concept of the complex potential enables the determination of the
influence of each of the components of the flow field, namely the flow
itself and the vortices and their images, at a given location. The
equation of motion of the vortex sheet, which is force free and
therefore moves with the fluid, is:

55
n _ 3W 3.2.3
at ' azn
at each point zn of the sheet, where (dW/dz) excludes the singularity
due to the vortex element at z . From Equation 3.2.3 we derive that:
9C
3.2.4
3z
In the case of a core vortex, the vortex which represents the centre of
a continuous spiral or cluster and whose strength grows with time
according to the rate of amalgamation of vorticity into that core, a
modified zero force equation is employed (Graham, 1980). Equation 3.2.4
is rewritten as:
3W] K 2 1 3rn rac

bn1 3.2.5
3t rn at I n Zsj
3zn L3znJ
where zn - zs is the cut joining the core vortex to the rest of the
sheet. Equations 3.2.4 and 3.2.5 are used in conjunction with Equation
3.2.2 for the determination of the convection velocities of the vortices
in the transform plane. With the inclusion of a core vortex and Routh's
correction (Clements, 1973), the transform plane convection equation is:
iv+ I rm
ir
n 1 1 n f 9 2 Z nn ,/ n
+ ~ r ' 1 1
at 2tt
m/n
+C* 2ir n c +c* 4 tt
mj ^ ti ^nJ
3z.
_L
1
rn
—L
9Fn
L r z*-z*)
at K n s' 3z. 3.2.6
56
Equation 3.2.6 is employed in the Discrete Vortex Method via a
straightforward one-step first order Euler integration to yield a means
of convecting the discrete vortices through a time step.
3.3 Treatment of the Nascent Vortex
In discrete vortex representations where the vortex sheet is being fed
with new vorticity from a shedding edge, a mechanism is required for
determining the strength and complex coordinate of the new discrete
vortex which represents this additional vorticity over a discrete period
of time. Such a vortex is designated a "nascent" discrete vortex. The
Kutta condition has most often in previous work been invoked to solve
one of these three unknowns (vortex strength and complex coordinate); in
this work it has been used to calculate the strength of the nascent
vortex. There is then a requirement for the determination of the nascent
vortex position; this is described in Section 3.3.2.
3.3.1 Satisfying the Kutta Condition
The Kutta Condition is satisfied when there is a finite flow velocity at
the shedding edge in the physical plane. This implies because of the
singularity of the transformation at the shedding edge that the velocity
at this position in the transform plane is zero (i.e. dW/d£ = 0). The
general equation of motion of the fluid in the transform plane at any
location other than at a vortex itself is Equation 3.2.2. Satisfying the
Kutta Condition leads to an expression for the strength of the nascent

57
vortex, r , placed at C :
n r n
-2ttV - J r r i i yI f 1 1
3.3.1
V Jr m 1 1 1
m?n ^shed Sn ^shed+Srv-
fished ^n ^shed+^n-
where ^sêcj is the position of the shedding edge. For the case where the
shedding edge is at the origin this may be simplified to:
r = 2ttV - l rm 1 n1
n 3.3.2
m^n v^n sn'
3.3.2 Determination of the position of the nascent vortex
Various different methods of solution for the position in the flow of
the nascent discrete vortex were tested, including a sheet vortex
representation of the nascent vortex and a solution of the Brown and
Michael equations for the nascent vortex (Brown and Michael, 1955).
Eventually the simplest method, reported in detail by Kamemoto and
Bearman (1978), was selected. The nascent vortex was placed at a
specified point on the edge bisector, its position remaining invariant
throughout the time cycle. In sinusoidal flow, a suitable position based
on the peak flow velocity was found to be somewhat different from that
for uniform flow, which was analysed by Kamemoto and Bearman. In fact,
the position of the nascent vortex for satisfying the Kutta Condition in
the physical plane (zkut) was taken as:
3.3.3
zkut = c3 Cz
58
where was determined from numerical experimentation to have a value
of about 0.01, although the optimum value varied somewhat according to
geometry. The principal factor in the determination of c^ was the
requirement to maintain an approximately even spacing between the
vortices near the shedding edge. Equation 3.3.3 may then be used to
determine C (i.e. C. ., the transform plane equivalent of z. .) in

n kut K M kut
Equations 3.3.1 or 3.3.2.
It has been noted by Naylor (1982) that instabilities quickly arise in
computations if one nascent vortex is introduced into the calculation
per shedding edge at each time step. Naylor reported that by placing a
new nascent vortex in the flow less frequently, say once every two or
three time steps, these instabilities were reduced. The drawback of this
approach is that computation time is increased without any concomitant
increase in definition of the vortex sheet. The instability problem
appears to be related to the first order time integration scheme used
for convecting the nascent vortex, in that because the nascent vortex is
in a region of very high flow velocities it can be convected an
unrepresentatively large distance. Having more than one convection step
per nascent vortex increases the accuracy of convection. A less
computationally expensive solution to the problem has been adopted in
the present work. Recognising that in the formulation described in
Section 3.3.1 the Kutta Condition is only satisfied at the beginning of
a time step, and not at any other point in the time step while the
nascent vortex is being convected away from the edge, it is clear that
the Discrete Vortex Method must be somewhat insensitive to the precise
location of the position of introduction of the nascent vortex. It is

59
advantageous to place the vortex at a position slightly further away
from the shedding edge than would be required to satisfy the Kutta
Condition at the beginning of the time step, because here flow
velocities are lower and the vortex is actually generally convected a
shorter and more representative distance than otherwise. Furthermore,
introduction of the nascent vortex even slightly further away from the
edge in the manner described above dramatically reduces instabilities.
Location of the vortex in this manner is justifiable on the grounds of a
degree of insensitivity to the position of introduction of the vortex;
numerical experimentation confirms that indeed this location is not
critical. This technique has been employed in all DVM calculations
reported here. The position at which the nascent vortex is actually
inserted into the flow in the transform plane ^ nascent^ 1S 9 lven in
terms of the position in the transform plane where the Kutta Condition
would be satisfied at the beginning of a time step by:
^nascent c4 Scut
c^ was held between 1.25 and 1.4 depending upon the stability
requirements.
3.4 A Scheme for Amalgamating Vortices
In the present model of vortex amalgamation each successive discrete
vortex nearest the central core is amalgamated with it so as to create a
core of increasing vortex strength (see Figure 3.4). Amalgamation has
been employed for two reasons. First, in connection with the tendency to
60
randomisation of the vortex sheet which has been noted by a number of
workers (see Chapter 1) amalgamation has been shown significantly to
retard the randomisation process (Moore, 1974). Secondly, amalgamation
saves computing time since the total number of vortices in the flow can
be reduced. In fact, for oscillatory flow the computation will reach a
"steady state" where the number of vortices in the calculation will
remain roughly constant since the rate of production of new vortices
will equal the rate of amalgamation. Typically, a new vortex from an
initial cluster of, say, 30 vortices would be amalgamated once every 3
or 4 time steps. This would have the effect of limiting the total number
of vortices in a flow with 64 time steps per cycle to about 60 or 70.
In the above model two discrete vortices are amalgamated so that their
new combined position, (x + iy) , is given in terms of their old
positions and circulations, designated by subscripts "oldl" and "old2"
as follows:
3.4.1
Combined circulation is given by:
rnew roldl + rold2 3.4.2
Deffenbaugh and Marshall (1976) have laid out the conditions for correct
amalgamation. These are firstly that it yields an equivalent velocity
field and secondly that it maintains the basic structure of the vortex
distribution. The above method fulfils both of these conditions.

61
3.5 Starting the Flow Pattern
An initial superposition of velocity off the edge was used to start the
vortex pairing pattern at the beginning of the computation. This
velocity was decayed quickly with time so that it only significantly
affected the first time cycle of flow and had no long term effect on the
solution. The potential in the transform plane due to this velocity was:
W(C) = Uinitial 3.5.1
where U. -j was chosen so that vortices from the first half-cycle
were convected well clear of the shedding edge. Thus by the end of two
time cycles the velocity had decayed to only 0.03% of its initial value,
initial*
3.6 Decaying the Vortex Pairs
Many workers with the Discrete Vortex Method have employed a mechanism
of removing vortices if they approach the body too closely - see, for
example, Clements and Maull (1975). Otherwise, these vortices can
convect at high speed close to the body surface due to the influence of
their images and possibly interfere seriously with the shedding
mechanism. A slightly different approach has been adopted with the
present method - all vortex pairs are slowly decayed in strength once
they have reached a certain degree of development, whether or not they
near the body side. It was explained in Chapter 2 that in sinusoidal
flow vortex pairs develop which have elements of nearly equal strength
62
but opposite sign. Thus their decay only has a small effect on total
circulation levels in the flow.
The decay mechanism chosen was programmed to come into effect as soon as
the latter member of the vortex pair had ceased growing and had started
to convect away from the edge. The decay was of the form:
-k (t/T)2
r = r * •. • i e a 3.6.1
initial
The above method was entirely successful in ensuring that vortex pairs
did not approach the edge too closely, or otherwise interfere too much
with the shedding process in cases where they were slow to convect away.
Since the circulation decay rates employed were low (k^ always less than
0.3), and decay was only employed for vortex pairs clear of the edge,
there was a negligible effect on the computed vortex force, as was
demonstrated by the results of computations with different values of k^.
Unlike the work of some other researchers, circulation decay was not
used as an expedient for reducing the magnitude of the vortex force. For
instance, Kiya, Arie and Harigane (1979) in their calculations of flow
past a flat plate normal to the flow noted that the DVM without the use
of vortex decay overestimated the amount of circulation in the wake and
consequently the drag force. They employed vortex decay successfully as
a means of reducing the circulation in the wake nearer to representative
values. In the present work, however, it was accepted as part of the
computations that forces would be overpredicted, it being thought wiser
to take account of this at the final stages of analysis of the results.

63
3.7 Calculation of Pressure Coefficient
Bernoulli's equation for unsteady incompressible flow is:
3<J> 2
p + |p(u2+v2) = P„ + ipUo o 3.7.1
3t
whence we obtain the following expression for the pressure coefficient:
r _ I (u2+V2) 2
3.7.2
p 3t
V2 V2
z z
with:
1 1 H
(u 2+v2) = |"v + J - Jr 3.7.3
2ir £ m T-T C+C 3z
m 1 h h ttv
and:
3£ =
- Real ( J- y r 1— !--- — + — Ur ---
3t l^ i 3t 3t
1 3rn 3.7.4
+ 2? I F &rg(?+C*) - ai"9 ( ^ n)l
where the subscript "n" refers here to the nascent discrete vortex.
Since there is in general no definable freestream velocity in the
physical plane, Vz is employed as a reference velocity. The above method
of determining the d^/dt term follows from the work of Jaroch (1986).
64
3.8 Calculation of Vortex Force
Following from the work of Graham (1980), the vortex force, F , may be
obtained from an approximation of Blasius1 equation for unsteady flow so
that for isolated edge flows:
A vortex force coefficient may be derived in terms of the length scales
in the physical and transform planes as follows:
CFv ■ Fv/(§pL* Lc T'2 ) 3,8,2
The values of force coefficient can be translated into drag and inertia
coefficients for the isolated edge by taking the appropriate Fourier
integral over one cycle of the flow to give:
C 3.8.3
D
and:
C = C „ + B, where B = X f' Cc cos d(|) 3.8.4

m mO tt2 jo Fv T ----
3.9 Extension of the Discrete Vortex Method to Model the Vortex Shedding
from Rolling Ships
3.9.1 Vessel in pure sway
Extension of the basic DVM to cope with complete symmetric ship
cross-sections in pure sway involved mainly modifications to the basic

65
numerical Schwarz-Christoffel transformation described in Section 3.1
above. An algorithm was developed to determine the exterior edge angles
(a.'s ) of a general hull cross-section composed of small segments as
shown in Figure 3.5. A suitable total number of segments for adequate
representation of the geometry without excessive computational expense
was determined to be roughly 15, although the number varied somewhat
according to the degree of complication of the geometry. The a^'s were
obtained from raw sectional data using scalar products. Vector products
were also calculated to assure the angles were attributed the correct
sign.
The flow velocity in the transform plane (V) is found to be related to
the freestream velocity in the physical plane ( I M by the following
relationship:
V = MU00 3.9.1
---------
3.9.2 The representation of roll
The roll of a vessel can be represented by a transformation method as
follows. Firstly, the ordinates themselves should be constantly changed
to take account of the changing position of the vessel; secondly the
motion of the vessel should be represented by a distribution of
singularities (in this case sources and sinks) along the vessel surface.
This second part of the representation may be seen as a modification of
the surface boundary conditions to yield a normal surface velocity
distribution which represents the motion due to roll. For small roll
66
amplitude the changing position of the vessel can be neglected,
resulting in significant economies in programming and computation time.
This simplification has been made here.
Flow issuing normal to the body segment can be represented by a source
distribution of strength per unit length a(z) where the transpiration
(outward normal) velocity is given as:
A general complex velocity, u - iv at z, is given by:
g(z')1d z '1
u iv 3.9.3
21tt t 1-2 '
as shown in Figure 3.6. The source distribution may be transformed to
the transform plane to yield a general equation for the velocity in the
transform plane:
Idz'l dz' dc'

u iv = J2— f° a(z')
Jr v •
3.9.4
tt cTc^* C-C'
as shown in Figure 3.7. If segment ab in the physical plane is a body
segment, and the transformation maps the body into a vertical line in
the transform plane, the transpiration velocities in the transform plane
will be purely real.

67
3.9.3 Method of solution for the velocity distribution due to a source
distribution
For any but the simplest transformations or source distributions,
Equation 3.9.4 has to be solved by numerical means. The process of
numerical quadrature involved is complicated by the fact that Equation
3.9.4 has a singularity at C = V • A sophisticated NAG library routine,
D01AJF, was selected. This is an adaptive quadrature routine, using the
Gauss 10-point and Kronrod 21-point rules. It was tested most
successfully on certain cases for which analytical solutions were known.
3.9.4 Application of the theory to a rolling vessel
The various aspects of the implementation of the theory described in the
sections above are presented below.
1) Roll motion was divided into two portions, pure roll about a fixed
roll centre (which was actually fixed at the base of the vessel on the
vessel centre line), and a sway component adjusted to yield the actual
desired roll centre, as illustrated in Figure 3.8. Freestream velocity,
U^, is given by:
= ojzR rj^sin w t 3.9.5
where w is roll angular frequency and n is roll amplitude. zD is the

k K
height of the roll centre above the vessel base. Roll angle was defined
68
as positive clockwise. The advantages of fixing the roll centre in this
way are twofold. First, it is convenient for calculating the roll
velocities to use a roll centre at the origin of the coordinate system.
Secondly, the computed roll velocities (i.e. due to the source
distribution) at and near the shedding edge will often be nearly
minimised by this arrangement. This means in practice that the order of
accuracy of the computed roll velocities is not required to be as high
as would have been the case otherwise.
2) The source strength on each segment is calculated using the value of
surface (transpiration) velocity at the midpoint of the segment with the
source strength assumed constant over the segment. This simplification
thus requires that velocities only change slowly along the length of the
segment or that segment size is small. The surface velocity (Vy) is
taken as the normal component of the total velocity. Thus the slip
velocity (velocity parallel to the surface) is not accounted for. This
is quite correct for a potential flow model such as the DVM, but would
be incorrect if a full solution of the Navier-Stokes equations were
attempted.
3) The velocities due to roll were required at two separate stages
within the computer program. In the first place, the velocity at the
shedding edge due to roll was required to be calculated once and for all
at the beginning of the program. The amplitude of this velocity does not
change with time. It was required to be calculated with fair accuracy,
since vortex strengths were partly dependent on it. For this velocity
the NAG routine integration technique D01AJF was employed as described

69
above. The value of the transformation derivative (dz/dc) was required
at each integration point. Secondly, velocities due to roll were
required at every time step for each vortex. Since this would have been
computationally very expensive using the NAG routine, an approximate
analytical technique was developed. This technique made the assumption
that the transformation derivative was constant over the segment, and
equal to its value at the midpoint of the segment in the transform
plane. It was possible to derive the following equation for a general
velocity, u - iv, in the transform plane at a point c due to a segment
cd:
acd fdz dz U -o
u iv 3.9.6
2tr (dC dz
cd
The above equation was derived by recognising that in Equation 3.9.4
acc|, (dz/dc)cc( etc. are assumed constant for the entire segment and
represent the value of the appropriate variable at the midpoint of the
segment. The velocity due to all the segments is simply the sum of the
velocities due to each segment. However, the symmetry of the vessel can
be used here with advantage. We know that if the source strength on
segment cdr on one side of the vessel is a £ , the respective source

r 1
strength on the other side (a ^ )
is simply - a ^ . Thus a general
p __ 1
velocity at a point £ due to a segment eel and its image segment cd is:
CL .
r '
Q
dz dz ~(gc-c)(ej - 0 "
O
u-iv= log 3.9.7

2tt
cd ‘ dz cd s ^ ) u ' c -0..
70
4) Flow length scale, Lz , for the rolling ship case, is based on the
total flow velocity at the shedding edge in the transform plane. It
contains both the sway and the roll components.
3.9.5 Calculation of the vortex roll damping moment
The vortex roll damping moment (per unit length), M , is calculated

v5
using the method described in Downie, Bearman, and Graham (1984). This
is essentially a Blasius momentum balance similar to that for the
isolated edge flows (Equation 3.8.1). Assuming small amplitudes of
motion:
Re al{ (z -z R) dz}
- -irJL V r ^ m + ^m) 3.9.8
" " lp 9t l Fm
m 71^ ^shed^
where c ^ ^ is the position in the transform plane of the shedding edge,
and the limits "a" and "b" denote the two junctions of the hull
cross-section with the free surface. The moment arm, 1 , of the vortex
a
damping force about the true roll centre is given as:
____ 1
b Real { ( z - z R) dz}"
cr
Q.
/
N
3.9.9
Ja (^shed> / .
a ^ ^shed^ —
71
Sectional vortex damping coefficient, t>v , is defined as:
M *1
v s
b 3.9.10
v
where 1 is the section length. Total vortex roll damping coefficient,
Bv , is obtained by summing the sectional damping coefficients over the
length of the ship. Implicit in this analysis is the assumption of strip
theory that for flow calculations a 3-dimensional body can be divided
into a series of 2-dimensional transverse strips. In the case of vortex
shedding, this is a fair assumption provided that changes in the
sectional geometry occur gradually over the length of the ship and that
end effects (i.e. at the bow and stern) can be ignored. Actually, in
many practical cases, these conditions are not completely fulfilled. It
is clear, for instance, that for a ship where the majority of the vortex
damping is generated from a few sections near the bows (and possibly
also the stern) end-effects could be important. In this case, strip
theory might be expected to provide an over-estimate of roll damping.
3.9.6 Range of validity of the method
Various simplifying assumptions are incorporated into the above method.
These are:
1) The vessel is symmetric about its centre line.

72
2) The vessel undergoes a roll motion of small angular amplitude and the
length scale of the vortex shedding is small in comparison with the
length scale of the vessel. These two closely related assumptions are
important not only in the implementation of roll itself but also in the
use of the particular formulation of Blasius' theorem for the
calculation of force which is employed here, which requires that the
vortex shedding is localised to the shedding edge.
3) The segment lengths are small in comparison with the perimeter of the
vessel.
4) The shedding edge is sufficiently sharp that the assumption of
separation point fixed on the edge bisector is valid or approximately
valid.
5) In the case of vortex shedding from more than one shedding edge, it
is assumed that the interference between the vortex structures from each
edge is negligible. This enables the method to compute separately the
damping contributions from each source; it is unlikely to be a major
restriction for such practical cases as shedding from bilge keels.
6) It is assumed that it is approximately valid to use a strip theory
approach for the calculation of total vortex damping moment.

73
CHAPTER 4
RESULTS R O M THE DISCRETE VORTEX METHOD AND COMPARISON WITH
EXPERIMENT FOR ISOLATED EDGE FLOWS
4.1 Assessment of the Sensitivity of the Method to Changes in the Run
Parameters
Variations in the run parameters inherent in the method can in certain
circumstances produce significant changes in force levels. An analysis
of the effect of some of these parameters has therefore been carried
out.
4.1.1 Effect of vortex decay rate
Excessive circulation decay rates applied to vortex pairs of slightly
differing strength can affect the strength of the nascent vortex
sufficiently to induce an instability into the calculation. They also
have the effect, if applied to vortices too near the shedding edge, of
reducing the magnitude of the vortex forces. Results have shown that
drag values are practically unaffected for either of these two causes by
the particular implementation of circulation decay in the current method
(Section 3.6). In one particular case (square edge with bilge keel), the
decay rate was reduced tenfold from its normal value with no change in
drag in the third significant figure. Some calculations were performed

74
for the plain square edge with no circulation decay at all. Until the
point when a vortex pair returned near to the edge, drag values were
virtually identical to those obtained in calculations using decay.
4.1.2 Effect of the position of the nascent vortex
The Discrete Vortex Method is to some extent self-adjusting to the
position (Z|<ut) of the nascent discrete vortex for satisfying the Kutta
condition, through a feedback mechanism via the vortex strength in the
shear layer. So long as the nascent vortex is introduced near the
midpoint between the shedding edge and the previously shed vortex,
changes in its position do not affect the drag excessively. For
instance, for shedding off an edge of zero internal angle a doubling of
zkut from usual position made a 1% difference to the results.
However, for a square edge, or edges of even greater edge angle, drag is
only independent of within narrower limits. It appears that an
excessively large value of z ^ has the same sort of action as a small
bilge keel in increasing the drag and modifying the shedding behaviour.
An excessively low value of the quantity leads to vorticity being
convected along the body side under the influence of its image system
rather than being convected into the flow. In a typical comparison for a
square edge, a 40% reduction in z^ ^ from a reasonable value to an
excessively low value created a 15% reduction in drag and instabilities
in the solution. Care was therefore taken to ensure for plain square
edge calculations that the shedding position was as close to the edge as
possible, but not too near as to create the above mentioned problems.
75
4.1.3 Effect of time step size and amalgamation rate
These two parameters are closely related - time step size affects the
overall representation and behaviour of the vortex sheets, whereas
amalgamation rate affects more specifically the behaviour of the vortex
sheets once they begin to convect away from the edge. Amalgamation rate
is defined as the number of vortices amalgamated per cluster per time
step. A low amalgamation rate retains greater definition in the vortex
sheets as they convect out into the field.
Two sets of numerical calculations have been performed to test the
effect of refinement of both of the above parameters. All the
calculations were performed with a value of a/l_z of unity with force
averaged over 8 time cycles. In the first set the number of time steps
per cycle was kept constant at 64, the number used in most of the rest
of the computations; amalgamation rate was decreased from 1 amalgamation
every 2 time steps for a vortex cluster to 1 amalgamation every 4 time
steps. The results of these computations are given in the table below:
Amalgamation Rate Coefficient A
0.500 6.453
0.333 6.552
0.250 6.443
76
In the second set amalgamation was held constant at 1 amalgamation every
3 time steps (the value used in most of the rest of the computations)
whilst the number of time steps per cycle was varied from 40 to 88. The
results of these computations are given in the following table:
Number of time steps per cycle Coefficient A
40 6.428
64 6.552
88 6.537
The results show that for shedding off a keel attached to a 90° edge
changes in the amalgamation rate and number of time steps per cycle over
the range tested (which spans the range used in all of the computations)
have no significant effect on the value of the drag coefficient.
Similarly, no significant effect was observed on the value of inertia
coefficient.
For the plain square edge also, changes in the refinement of the vortex
structure within the range described above did not produce a noticeable
effect.
77
4.2 Results for Square Edge with Bilge Keel and Optional Rounding
4.2.1 Vortex positions and pressure distributions
A series of four pressure distributions, together with the appropriate
plots of vortex position, is shown in Figures 4.1 - 4.4. These present
distributions for a square edge with keel of span a/l_z = 2.2 but no
rounding. The first figure shows the situation at the beginning of a
time cycle (t'/T = 0); the second figure the situation 1/4 of the way
through a time cycle (t'/T = 0.25), and so on. The qualitative agreement
between the theoretical plots of vortex position and experiment may be
assessed by a comparison with Figure 2.6. The shedding from the left
hand keel is in phase with the computation. No experimental pressure
distribution results are available, but the general form of the
distributions seems correct. A suction peak is observed to grow on the
surface adjacent to the nascent vortex cluster as the nascent vortex
(cluster) grows. This is especially obvious in Figure 4.2. The pressures
on the two surfaces should coincide at the shedding edge. In each of the
distributions as the shedding edge is approached the two halves of the
distribution start to come together before diverging again rapidly very
near the edge. This is interpreted as a feature of the Discrete Vortex
Method in that the discretisation of the vortex sheet means that the
expected discontinuity in the velocity distribution at the shedding edge
is "smeared" over a small distance in the theory, with consequent
effects on the pressure distribution.

78
Figures 4.5 show four plots of velocity vectors for the case
r/Lz = a/l_z = 1.0. The plots are for one complete time cycle, each plot
being separated by a quarter of a time cycle. A velocity vector is drawn
from the centre of each discrete vortex. Figure 4.5a shows a vortex
cluster growing at the edge. Away from the edge a vortex pair is
convecting towards the bottom right of the picture. In subsequent
pictures the vortex growing at the edge pairs up with a new nascent
vortex and this pairing also convects away from the edge, following the
previous one.
Figure 4.6 shows two plots of vortex positions half a time cycle apart
for the case of r/L = a/Lz = 0.039. Together with these figures are
photos of the corresponding flows on the barge. In comparison to the
keel span of the previous calculation this is a comparatively short
keel. Nevertheless, even a short keel such as this has a profound effect
on the vortex shedding pattern compared with the shedding off a plain
square edge.
4.2.2 Comparison between theoretical and experimental vortex size
In Figure 4.7 experimentaland theoretical vortex sizes are compared
plotted against non-dimensional keel span for a square edge with keel.
The apparent divergence between theory and experiment at the highest
keel span is most likely due not to any deficiency in the Discrete
Vortex Method but due to substantial uncertainties in the measurement of
vortex size, lv , in experiment. A method of matching the theoretical

79
isolated edge results to the experimental barge rolling results was
employed.
4.2.3 Time dependent force coefficient
Plots of the variation of force coefficient, Cpv , with time are
presented in Figures 4.8. Figure 4.8a is for a square edge with bilge
keel of span a/l_z = 0.2; Figure 4.8b is for a plain square edge. Both
figures show traces taken over four time cycles, beginning six time
cycles into the flow. The more structured vortex roll-up for the edge
with the bilge keel, and more random shedding pattern for the plain
square edge are reflected in the difference between the two force
traces. Peak force in the first case is consistent over the four cycles,
but varies quite dramatically in the second. A sudden trough in the
third cycle of Figure 4.8b indicates the presence of a vortex very near
the body side being convected at high velocity.
4.2.4 Drag and inertia coefficients
In the present formulation, drag coefficient (Cp) is equal to the
coefficient "A" whilst inertia coefficient (Cm ) is related simply to,
but not strongly dependent on, "B" (Equations 3.8.3 and 3.8.4). Values
of the coefficients A and B for varying non-dimensional bilge keel span
may be found plotted in Figures 4.9 and 4.10 respectively. Results are
presented both for a square edge with keel and for a rounded square edge
80
with "bounded" keel (i.e. where a/r = /2-1 - see insets on these figures
for a definition of a bounded keel). Also shown are the equivalent flat
plate values (i.e. square edge with infinitely long bilge keel). Both A
and B should be divided by a factor of 1.62 to provide comparison with
the actual values of flat plate force coefficient (i.e. infinitely long
flat plate) because of the matching technique involved. Keel span in
these two figures is non-dimensionalised with respect to the equivalent
length scale for the plain edge (Equation 3.1.17) for the sake of
simplicity in its use with the matching process. Over the range of keel
spans tested here, length scale was virtually invariant, and equal to
L / The form of the curves for both the coefficients is

z(p lain edge)
almost identical - a sharp initial rise at low keel span, followed by a
slow asymptote to equivalent flat plate values. In the case of the
rounded edge with bounded keel the drag results are, as would be
expected, considerably lower than those for the keel with square edge.
The edge rounding in the former case obviously reduces the flow
velocities near the shedding edge.
Values of the coefficients A and B are 1.995 and 0.051 respectively for
the square edge, and 5.38 and 0.27 for the flat plate. The relevant
values calculated by Graham (1980) are 1.57 and -0.04, and 5.9 and
0.125. Unlike the coefficient A, coefficient B is sensitive to the
numerical parameters in the discrete vortex analysis. Also B is
relatively unimportant in comparison with A. The most meaningful
comparison with Graham is therefore for A. Experimental results for A
for the square edge and flat plate are approximately 1.4 and 4.0 though
there is a degree of uncertainty attached to them. The experiments - see

81
Bearman, Downie, Graham and Obasaju (1985) - were carried out for
two-edged bodies (i.e. the diamond and the finite flat plate); the
results have been matched to the isolated edge via a process described
by Graham (1980). The matching process requires the use of experimental
results from the very low Keulegan-Carpenter number range,
Keulegan-Carpenter number (Kc ) being defined as:
U T
Kc = S - 1 ^ 1
where d is a representative body length scale (the diameter in the case
of a circular cylinder). Such a requirement implies that the scale of
the shed vortices is much less than the body scale.
4.3 Results for Square Edge with Rounding
This case is more difficult because of the uncertainty over where the
separation point(s) are. On the one hand in the limit of vanishing
radius we can fix the separation point with some certainty. On the other
hand, as radius is increased at fixed flow velocity, experiment has
shown that the vortex shedding pattern observed for the plain edge
(single pair per time cycle) gives way to a different pattern (double
pair per cycle - see Figures 2.7 and 2.8), with a constantly varying
separation point. It was not possible to determine from the experiments
that were carried out in what radius range this change took place.
Eventually vortex shedding ceases altogether. The correct simulation of
the double shedding pattern is beyond the scope of the Discrete Vortex
Method as it stands. However, when the radius is small it seems
reasonable that a fair representation of the vortex shedding may be

82
obtained by fixing separation on the edge bisector, although this is not
in general the true separation point. The results would be more
representative the smaller the radius is. The approach of fixing
separation thus has been taken here, the method being used for
increasing radius until the computations became unstable. Alternative
more sophisticated methods of determining the separation point could
have been developed. However, at this stage of the research it was
decided to proceed to full Navier-Stokes solutions of the flows around
these shapes (Chapters 6 and 7) and so it was decided not to devote more
development time to this method.
4.3.1 Vortex positions and pressure distributions
Results for three different edge radii are presented. Figures 4.11 to
4.14 refer to a non-dimensional edge radius (r/Lz ) of 0.039. The four
figures are separated in time by a quarter of a time cycle. The vortex
shedding is very similar to that off a plain square edge (see Figure 2.4
for the experimental results for a square edge) except perhaps that the
nascent vortex cluster keeps closer to the edge. In this case the
present method referred to above of fixing the separation point clearly
works quite well inasmuch as shedding patterns are well-ordered and
appear plausible. Figures 4.15 to 4.17 however show a situation where
this approach does not function so well. Non-dimensional edge radius is
0.109; each figure is separated by a time interval of half a time cycle.
The computation appears to be attempting to simulate the double pairing
pattern referred to in the introduction to this section. The large

83
vortex cluster at the edge in Figure 4.16 has been split in two half a
time cycle later in Figure 4.17, the smaller (residual) portion of the
vorticity pairing up with the new vortex cluster near the edge. Despite
some short periods (such as here) where the computation runs without
major problems, as a whole it is highly unstable and seems to be on the
point of complete breakdown.
Figures 4.18 to 4.21 are for a non-dimensional radius of 0.086. Here the
radius is somewhat smaller and the calculation is more successful. Each
figure is separated by a time interval of a quarter of a time cycle. The
relevant pressure distribution is shown next to each plot of vortex
positions. The general form of the pressure distributions is similar to
those of the plain edge except that the magnitudes of the pressure
values tend to be somewhat smaller. For instance a suction peak is still
discernible in the case of a strong single nascent vortex cluster
(Figure 4.21 - upper surface). There is also a suction peak at the
junction between rounding and straight edge on the surface opposite that
on which the nascent vortex is growing. This is followed by a strong
recompression. This suction peak can be seen in all four cases, but is
strongest when the nascent vortex is strongest (Figures 4.19 and 4.21).
The feature is interpreted as a local acceleration of the flow due to
the curvature discontinuity between the straight edge and the rounding.
84
Because of the incorrect position of the separation point it is thought
that there will be an error in the drag magnitude which will increase
with radius. Nevertheless drag values are presented here because it is
felt that they give some qualitative guide as to its behaviour and also
provide an indication of the limitations of the Discrete Vortex Method.
The drag results presented in Chapter 7 and calculated using the Hybrid
Moving Vortex Diffusive Method are considered much more reliable. Figure
4.22 shows a plot of the coefficient A against non-dimensional edge
radius. A value of zero radius corresponds to a plain square edge. As
radius increases, drag decreases as we should expect. No results are
presented beyond a value of r/l_z of 0.064 since drag results for
calculations at higher r/l_z are considered unreliable. In these cases,
spuriously high values of drag were obtained when the vortex pairing
switched direction, as it was prone to do frequently for the larger
radii. This led to a sharp levelling off of the decreasing drag trend
which was considered spurious. No plot of the coefficient B is presented
since no discernible trend can be found in the distribution. The
computed values are small and, except for the plain square edge,
negative.
85
CHAPTER 5
APPLICATION OF THE DISCRETE VORTEX METHOD
TO ROLLING SHIPS AND BARGES
In the context of the current research, force data derived from the
Discrete Vortex Method were required to be applied to the problem of the
roll damping of vessels due to vortex shedding. Two different
applications were made. Firstly, the isolated edge results for an edge
with keel and/or rounding were applied via a matching process to the
problem of barge roll damping. In this case the results were input into
a specially modified version of the vessel motion response computer
program named BMTIMP in order to obtain roll damping coefficients. This
work will be described in the following sections. In the second place,
the Discrete Vortex Method as developed for complete hull cross-sections
was used directly in the computation of roll damping coefficients for
the Fisheries Protection Vessel, the Sulisker. These results will be
described in Section 5.5. In this section computed vortex shedding
patterns are also presented and discussed.
5.1 Barge Motions and the BMTIMP Barge Rolling Computer Program
A schematic representation of barge motions is shown in Figure 5.1. The
barge undergoes motion with six degrees of freedom. These can be divided
into two coupled sets, one for sway, roll and yaw, the other for surge,
86
pitch and heave. Interest in the current work lies only with roll motion
and the roll damping coefficient, B^. Damping calculations were
performed exclusively for a barge in forced-roll at zero forward speed.
Response to waves was not considered.
The computer program now called "BMTIMP" was developed at Imperial
College under joint support from BMT Ltd. and the SERC Marine Technology
Directorate. It was formulated to enable the calculation of damping
forces experienced by a rolling flat-bottom barge of rectangular
cross-section due to the vortex shedding from its bilges. It was a full
boundary integral method involving the application of linearised free
surface boundary conditions, solved by strip theory. Since such a method
is based on potential flow theory, it cannot itself model flow
separation off the bilges and thus singularities arise there. The
innovation of BMTIMP is in the treatment of these singularities. The
flow is divided into an inner region close to the edge, consisting of
the edge singularity and the consequent vortex shedding, and an outer
inviscid region, consisting of the rest of the flow field. The inner
region was modelled by use of the Discrete Vortex Method applied to an
isolated edge, and was then matched to the outer region which was solved
by the potential theory referred to above. BMTIMP has been well
documented - see Downie, Bearman and Graham (1984). It could originally
cope only with barges of rectangular cross-section, but has been
modified by the present author to extend its scope to rectangular
cross-section barges with bilge keels and/or bilge rounding. The
modifications pertained mostly to the matching process and are described
in Section 5.1.2 below.

87
5.1.1 Treatment of the roll damping contributions by BMTIMP
It is convenient at this point to outline the manner in which the
various contributions to roll damping have been treated in BMTIMP. A
much fuller description will be found in Downie, Bearman and Graham
(1984) who originated BMTIMP. As was explained in Chapter 1, for the
problem of a barge at zero Froude number rolling at small amplitude,
there are only two significant contributions to roll damping, namely
wave (radiation) and vortex damping. Thus the roll damping coefficient,
B ^ , can be written as:
^44 = + 5.1.1
where Bj and By are the wave and vortex damping coefficient components.
Damping coefficient is derived by dividing the damping by the roll
velocity amplitude, defined as coin . The former is analysed traditionally

1+
through the use of strip theory and is independent of roll amplitude.
Generally in previous work the latter, as discussed in Chapter 1, has
either been neglected altogether or been approximated by empirical or
semi-empirical means. Here, we derive By from the discrete vortex
analysis, which provides a completely theoretically based representation
of vortex shedding (see Section 3.9.5); the coefficient possesses linear
or near-linear dependency on roll amplitude and frequency. This gives
rise to the non-linear, near-quadratic vortex damping term.

88
5.1.2 Matching the isolated edge flow to flow round the rolling barge
The flow field about a rolling barge comprises an inner region in the
vicinity of each keel edge, where the effects of flow separation and
vortex shedding dominate, and an outer region associated with the motion
of the barge and the effects of the free surface. The flow in the outer
region may be found from potential flow theory. Provided the length
scale of the flow round the shedding edge is small compared with a
typical barge dimension (e.g. beam or draught) the flow in the inner
region may be determined by a discrete vortex analysis of oscillatory
flow about an infinite wedge with appropriate edge geometry, matched
locally with the exterior flow field. The power of the matching
technique lies in the fact that for one particular edge geometry and
flow length scale only one discrete vortex calculation needs to be
performed regardless of the overall barge geometry involved.
The principles of the matching process have been described by Graham
(1980) and developed for the problem of a rectangular cross-section
barge by Downie, Bearman and Graham (1984). In this section we describe
the extensions necessary to this process to enable results for more
general shedding geometries to be incorporated into the BMTIMP barge
rolling program. In particular we are interested in square edges with
keels and/or rounded bilges.
The matching technique involves matching of the length scale of the flow
round the isolated edge to that of the flow round the rolling barge. For
89
the square isolated edge all flows are self-similar so that the vortex
force for any barge of rectangular cross-section can be obtained by
scaling from a unique value of isolated edge drag coefficient. With the
inclusion of a bilge keel or bilge rounding we introduce a dimension
into the isolated edge calculation. Now there will be a unique value of
drag coefficient for every value of bilge dimension (e.g. keel span or
edge radius) divided by length scale. This point is demonstrated in
Figures 4.9 or 4.22. Therefore additionally in these cases we have to
match the value of bilge dimension divided by length scale for the
isolated edge to that for the rolling barge, and a number of discrete
vortex calculations are required to cover the required range of bilge
dimension to length scale ratios for each geometry.
In applying the matching technique two important simplifying assumptions
were made.
1) It was assumed that, for identical external flow conditions, the
length scale of the flow is independent of edge geometry for the
rectangular barge with bilge keel and/or rounding. This assumption
places an upper bound on the permissible values of (a/!_z ) or (r/l_z ); it
is necessary for an efficient implementation of the matching process
which matches the results for the isolated edge to those for the barge.
The flow length and velocity scales, Lz and V , vary by less than 10% in
total over the range of bilge keel span, a, or bilge radius, r, for
which drag results were presented in Chapter 4. This means that the drag
results presented in this thesis for the edge with bilge keel or bilge
rounding can be incorporated into the barge rolling problem by the

90
matching process as if they were results for the plain edge, provided
the values of the coefficients A and B are selected via the correct
values of a/Lz(barge) or r/l_z(b e) using the graphs in Figures 4.9,
4.10 and 4.22. In the case of the rounded bilge, B is taken as zero
except for r/Lz (barge ) = 0. The above assumption is unlikely to be a
restraint in any physically realistic situation except at very small
roll amplitudes.
2) It is assumed that the bilge keel or rounding are details of the
inner flow and do not affect the external flow around the rolling barge.
This assumption is valid when not only the direct effect of the bilge
appendage on the external flow but also the disturbance of the external
flow by vortices shed from the appendage is negligible. These conditions
are met when keel span or bilge radius and roll amplitude are small.
This is essentially a restatement to include more arbitrary bilge
geometries of the basic matching assumption already inherent in the
method that the flow length scale round the shedding edge is small
compared with a typical barge dimension. The assumption is particularly
useful since it would mean that only that part of the barge rolling
computer program concerned with vortex shedding need be modified.
Modifications to the potential flow solution or the radiation damping,
for instance, need not be made. The assumption is unlikely to be a
serious further restraint on the validity of BMTIMP, since it is only a
restatement of the basic matching assumption of BMTIMP.

91
5.2 Results for a Rectangular Cross-Section Barge with Bilge Keels
Calculations were carried out over a range of roll amplitudes for a
force-rolled rectangular cross-section barge with and without bilge
keels. The basic barge was the Noble Denton Consortium standard barge
Case 3S with a length of 87.8m and height of roll centre above the water
line of 7.17m; the sectional dimensions are shown in the insert in
Figure 5.2. For all calculations the full scale roll period was 10s,
which corresponded approximately to the resonant frequency. Results are
actually presented at model scale, where the linear scaling relationship
is 30:1 full scale to model scale. For more details see Noble Denton and
Associates (1985). Figure 5.2 shows a comparison between experimental
and theoretical roll damping coefficients in the absence of bilge keels.
Although the experimental data has some level of uncertainty attached to
it, it is clear that theory only slightly overpredicts experiment.
Since in the case of the plain 90 degree edge without bilge keel the
vortex damping coefficient is linearly dependent on roll amplitude, it
follows that vortex damping coefficient per unit of roll as a percentage
of critical damping coefficient is independent of roll amplitude. For
the particular case of the barge Case 3S above it takes the value of
56.2% (see Figure 5.3). When a bilge keel is added, the linear
relationship above no longer holds. In fact, vortex damping per unit of
roll is found to be a function of ai/L2 (p *|a -jn ecjge ) s^nce the drag and
inertia coefficients A and B are functions of this quantity. It follows
from this and the expression for length scale for the square edge
92
(Equation 3.1.17 with A = 1.5) that vortex damping coefficient per unit
of roll as a function of critical damping coefficient is therefore a
function of (roll amplitude)/(keel span)4/3. This function is shown
graphically in Figure 5.3. This one graph contains, for this particular
configuration, all the information required to calculate the roll
damping for any combination of roll amplitude and keel span within the
prescribed limits. The information is used in Figure 5.4, where vortex
damping coefficient is shown plotted against roll velocity amplitude for
six different full scale keel spans, 0, 0.1, 0.2, 0.5, 1.0 and 2.0
metres. The roll centre and basic barge geometry are as described for
the barge Case 3S. No experimental data are available for comparison.
The curve for a keel span of zero is, as explained above, a straight
line. Increasing the keel span from zero raises the vortex damping as
should be expected. Smaller keels are relatively more effective. For a
given keel span the increase in damping over the plain edge case is
proportionately greater the lower the roll amplitude. This point is
shown in Figure 5.5, which shows the percentage gain in vortex damping
over the plain edge case against roll velocity amplitude for a set of
five keel spans.
The original computer program BMTIMP is quoted by Downie et al. (1984)
as being valid for maximum roll amplitudes of over 5 degrees in the case
of the barge tested here. The modifications to the program are also
valid over this range of roll amplitudes, provided that the assumptions
contained in Section 5.1.2 remain valid. These provide two restrictions.
First, assumption (1) requires that non-dimensional keel span (a/l_z )
should not be too great (i.e. should not exceed the highest value for
93
which isolated edge drag results have been presented). This would place
a lower bound on the roll amplitude for which the method is applicable,
but is not a practical restriction since even for the largest keel
tested (2m) the minimum admissible roll amplitude was found only to be
roughly 0.2°. The restriction arising from the second assumption is that
the effect of the keel on the external flow should be negligible. This
would provide an upper bound on keel span and roll amplitude. But it is
not thought to be a practical restriction. The assumptions described in
Section 5.1.2 are thought to permit a keel of at least one metre in span
over the entire range of roll velocity amplitudes presented here,
although with a lack of experimental data this figure is likely to be
very imprecise.
5.3 Results for a Rectangular Cross-Section Barge with Rounded Bilges
Calculations were performed for the standard barge Cases 2 and 3 by
matching round edge Discrete Vortex Method results. Both these barges
are essentially the same as Case 3S. Case 3 has rounded bilges of radius
0.48m; in addition, Case 2 has a length of 65.1m as opposed to 87.8m and
a slightly deeper draught. Graphs of predicted and experimental roll
damping coefficients are shown in Figure 5.6 for Case 3 and Figure 5.7
for Case 2. Also shown for the sake of comparison are the BMTIMP
predictions for a similar square-edged barge, and Robinson's empirical
corrections to these predictions to cope with bilge rounding - see Noble
Denton and Associates (1985). As would be expected, the effect of
rounding in the theory is to reduce roll damping. The general trend of

94
the experiment is predicted. However, the reduction does not in general
bring the prediction in line with experiment or with Robinson's
correction. Probably a number of reasons contribute to this, including:
1) There is quite a high level of uncertainty in the experimental
results.
2) The error in the theoretical results due to inadequacies in the
discrete vortex model used here for rounded edges is believed to
increase with non-dimensional bilge radius (see Section 4.3.2). That
section also explains that reliable drag results could not be obtained
beyond a value of (r/L ) of roughly 0.06. The region beyond which
reliable drag results could not be obtained is marked by a dotted line
on each figure.
3) Even for the square edge theory overpredicts experiment; the larger
the non-dimensional bilge radius the more the overprediction is expected
to be.
It is believed that the modifications to the BMTIMP program when used
for rounded bilges do not place any further restriction on its
operation.
Chapter 7 presents damping results for barge Cases 3S, 3 and 2
calculated using the Hybrid Moving Vortex Diffusive Method which was
developed because of the inadequacies of the present discrete vortex

95
model in computing flows round significantly rounded edges. The HMVDM
results are considered to be substantially more reliable for the barges
with rounded bilges (Cases 3 and 2). Discrete Vortex Method results for
these cases are presented in this chapter for the sake of comparison
only.
5.4 Results for a Rectangular Cross-Section Barge with Rounded Bilges
and Bilge Keels
Comparisons have been made with two barge models tested by Ikeda,
Komatsu, Himeno and Tanaka (1977). These are models B and C, both having
a length of 0.8m, beam of 0.28m, draught of 0.112m and roll centre at
the waterline. In each case keel span is half the bilge radius; this
radius is 1cm for model B and 2cm for model C. Graphs of predicted and
experimental vortex damping coefficients (By ) plotted against roll
angular frequency u at fixed roll amplitude are shown in Figures 5.8 and
5.9 for models B and C respectively. The ratio (a/l_z ) is independent of
frequency so that there is a predicted linear relationship between
vortex damping coefficient and frequency. In both cases agreement with
experiment is excellent and better than the comparisons for rounded
bilges without keels. The Discrete Vortex Method predictions are in
general more accurate the smaller the shedding edge interior angle. For
a rounded edge this angle is essentially 180°, but for a bilge keel it
is 0°.
96
Figure 5.10 shows the variation of vortex damping coefficient with roll
amplitude for model C at two fixed angular frequencies, 6.29 and 4.49
rad/sec. The theoretical curves are non-linear since the ratio a/l_z and
hence the computed vortex drag coefficient varies with roll amplitude.
Theory and experiment agree well at both frequencies, but slightly
better at the lower frequency. An examination of the variation of vortex
damping coefficient with frequency for this model (Figure 5.9) might
suggest that this slight apparent inconsistency relates to the
experimental, not the computational, results in that experiment is
locally somewhat higher than theory around a frequency of 6 rad/sec.
5.5 Results for a Typical Ship in Roll Motion
Section 3.9 describes the further developments necessary to enable the
Discrete Vortex Method to compute the vortex damping due to roll for
complete ships. A set of calculations was carried out using this method
for the Fisheries Protection Vessel, the Sulisker. These are described
in Sections 5.5.2 to 5.5.4 below.
As a test of this version of the DVM computations were carried out for
the standard barge Case 3S referred to previously in this chapter (see
Noble Denton and Associates, 1985). Comparison was made between the
direct calculations of damping for the complete vessel geometry and
those obtained via matching to isolated edge calculations. Respectable
agreement was obtained.

97
One further test case was carried out; this is described in the next
section.
5.5.1 The roll damping of a circular cylinder with bilge keels
A roll damping comparison was carried out between experiment and
Discrete Vortex Method predictions for the case of the circular cylinder
with two symmetrically positioned bilge keels rolling in an oscillatory
motion about its central axis.
The experimental results were obtained from model tests carried out by
Mercier (1965). A schematic view of the model cylinder with one of its
two identical keels is shown in Figure 5.11. The cylinder shown has an
outside diameter of 15.2cm. A second set of tests was carried out on a
larger cylinder with a diameter of 30.4cm, and equivalently proportioned
bilge keels. Roll period for the smaller cylinder was 1.46 seconds and
for the larger 2.04 seconds. The cylinders were suspended vertically on
a torsion wire and were given an initial angular displacement. Roll
damping moment was then deduced from the roll extinction (decrement)
curves. Mercier defined a roll damping moment coefficient, Crcjm j, as:
damping moment
^rdml 5.5.1
(P/2) Ak Rc ’ (,■*)*
where A, is the frontal area of both keels, and R is the radius from
k ca
the axis of rotation to the centre of area of the keel.
98
Mercier calculated C ^ j as 16 for the larger model and 15 for the
smaller.
The computation was performed with the version of the DVM as developed
for the prediction of the roll damping of vessels due to vortex
shedding. It was carried out with 80 vortices per time cycle over six
time cycles.
Roll damping moment coefficient results are as follows, based on
Mercier's definition of C , , :
rami
Experiments 16, 15
Discrete Vortex Method (vortex damping) 18.6
Assuming the vortex damping component to be the only significant
contribution to the overall damping - and Mercier shows that this is so
- the Discrete Vortex Method overestimates damping for this case by
between 16 and 24%.
5.5.2 The Fisheries Protection Vessel “Sulisker"
A series of discrete vortex computations was carried out for the
Fisheries Protection Vessel, the Sulisker. These are described in the
following sections; this section outlines the main characteristics of
the Sulisker.
99
The F.P.V. Sulisker has an overall length of 71.33m, beam of 11.6m,
displacement of 1532 tonnes, and nominal draft of 4.6m. A plan of the
ship may be found in Figure 5.12; cross-sectional data are to be found
in Figure 5.13. Detailed information concerning the roll testing of a
l/20th scale model of this vessel may be found in Spouge and Ireland
(1986).
Digitisation of the sectional data was undertaken by BMT Ltd.. The roll
damping experiments on the model of the Sulisker were carried out with a
bar keel of span 0.16m (full-scale dimension) fitted over a significant
portion of the ship (Sections 3 - 18). DVM calculations have been
performed both with and without a bar keel. Comparisons with other roll
damping prediction methods were only available without the keel.
Experimental results have also been obtained for the Sulisker with added
bilge keels. The bilge keels had a span of 0.35m and length of 8.68m.
They were positioned 3.06m below the waterline, 10m apart at midships,
and spanned Sections 12 to 15.
5.5.3 Vortex positions and velocity vectors
The Discrete Vortex Method was employed to compute vortex shedding
patterns for the Sulisker on a sectional basis in the manner described
in Section 3.9.5. The basic computation which was carried out was for a
roll amplitude of 8° with the bar keel and a roll centre 4.65m above the
ship base. Other computations investigated the effect of neglecting the

100
bar keel, the effect of two additional bilge keels, and the effects of
changing the height of the roll centre above the ship base (zR ) and roll
amplitude (n^). All calculations were performed at full scale. Precise
details of these calculations are given in the table below:
Run zR (m) n (°) bar keel bilge keel
(1) 4.65 8 Yes No
(2) 4.48 8 Yes No
(3) 4.65 4 Yes No
(4) 4.48 4 Yes No
(5) 4.65 8 No No
(6) 4.65 8 Yes Yes
(7) 4.48 4 Yes Yes
Because of the limitations of the method, it was not possible in Run 5
(no bar keel) to compute vortex shedding for all the sections since some
possessed no sharp shedding edge or possessed a sharp shedding edge of
internal angle greater than about 110°, and would have yielded
unrepresentative dampings. By far the greatest proportion of damping is
produced by the sharp sections, so that this is no significant
limitation on the validity of the results. This is confirmed by a
comparison with the sectional roll damping results of Brook (1987) which
were obtained using the roll damping method of Odabasi et al. (1985).
This method will be referred to as the "BMT" method and bases its
prediction of vortex roll damping on a process of matching to simpler

101
body shapes for which vortex drag coefficients had previously been
computed by Graham (1980) using the Discrete Vortex Method. In the
current work, computations for different roll centres have been
performed due to an uncertainty over the true position of the roll
centre in experiment. The roll centre was in fact believed to vary
somewhere between the waterline and the centre of gravity (i.e.
between 4.48 and 4.65) though even these limits may be too precise
(Spouge, 1987).
Plots of vortex positions for twelve of the sections for which DVM
calculations for Run 1 were carried out are presented in Figure 5.14.
The plots are all of the flow situation half-way through the 4th time
cycle. The size of the vortex cluster at the edge grows between Sections
1 and 3 due to the increased roll velocity at the edge. Between Sections
1 and 11 edge angle slowly increases. This is accompanied by a relative
weakening of the vortex strengths. The vortex clusters at the edge grow
nearer to the edge; the pairs are less structured as they convect away
from the edge, and tend to convect away nearer the body side. By Section
9 the vortex structure is quite weak; the pairing direction has also
switched, so that, at this time instant, instead of a single large
vortex cluster at the edge, there is a less well-structured vortex pair
convecting away from the edge. The least structured shedding patterns
are observed for Sections 11 and 13, with a few point vortices actually
having convected through the body side for Section 13 due to the
proximity of the vortex structures to the surface combined with a
convection scheme which is only first order accurate. For Sections 15 to
19 the shedding edge becomes sharper, with a consequent strengthening of

102
the vortex structure. The pairing direction preferred by the sharper
edges is re-established at Section 17. Between Sections 19 and 22 the
shedding edge again becomes blunter, whilst the draft of the sections
reduces, with an attendant diminishing of the size of the vortex
clusters.
Figures 5.15 to 5.18 show plots of velocity vectors for a variety of
conditions. The plots in each figure are separated by a quarter of a
time cycle each, starting at the beginning of the 4th time cycle.
Figures 5.15 and 5.16 compare the shedding patterns at Section 17 for
Runs 1 and 3 (roll amplitudes of 8° and 4° respectively). The size of
the vortex clusters is substantially smaller for the lower roll
amplitude as should be expected. Figure 5.17 shows velocity vectors for
the same section and 8° roll amplitude in the absence of a bar keel (Run
5). There is considerable similarity in the shedding patterns between
this run and Run 1, with a strong vortex pairing structure convecting
out into the flow field in both cases, but the absence of a bar keel in
Run 5 implies weaker vortex shedding, with vortex pairs staying closer
to the body side. Figure 5.18 shows velocity vectors for Run 7, Section
14. Vortices are being shed from one of the two bilge keels present. The
Discrete Vortex Method as developed here is only capable of calculating
vortex shedding from one shedding edge. Thus in the case of two or more
shedding edges on the same section, total vortex damping is obtained by
performing calculations for each shedding edge separately and then
summing the resulting individual vortex dampings. This assumes that the
interference from the vortex structures emanating from each shedding
edge is negligible (see Section 3.9.6).

103
5.5.4 Roll damping results
In this section damping results are presented in terms of the sectional
vortex roll damping moment coefficient, c ^ g, defined as:
'rdm2 = b / ( “n ) 5.5.2
Total vortex roll damping moment coefficient, C ^ obtained by
summing the sectional coefficients over the length of the ship.
Figure 5.19 provides a graph of sectional vortex roll damping moment
coefficient, crcjm2 5 ôr ûns * ' ^ anc* 5- ûn * yields dampings which
are consistently higher than those from Run 2, which has a lower roll
centre and would thus be expected to predict lower edge velocities and
therefore lower dampings. Use of the sectional vortex roll damping
moment coefficient permits direct comparison between damping results of
different amplitudes since damping has a near-quadratic dependence on
roll amplitude and c ^ ^ is a function of damping divided by the square
of roll amplitude. A comparison between Runs 1 and 3 on this basis
reveals small but interesting variations of c ^ ^ with roll amplitude.
It appears, in general, that a concave section geometry near the
shedding edge is associated with a lower damping at higher roll
amplitude, and that a convex geometry is associated with a higher
damping. The reason for this may be that a convex geometry effectively
becomes "sharper" the larger the vortex size (and hence the larger the
roll amplitude). This would then imply an increase in damping. The

104
opposite would be true for a concave geometry.
Comparison is also made in Figure 5.19 between the DVM results and those
due to Brook (1987) using the BMT method. Since the BMT method was
employed for a geometry withouta bar keel, the most meaningful
comparison is with Run 5. Both runs were computedwith the same roll
centre and at the same roll amplitude. Significant differences in the
distribution of damping over the ship may be seen, although both methods
agree roughly on the sections where most damping isbeing generated and
also on the magnitudes of thepeak dampings for the fore and aft
portions of the ship.
The sectional vortex roll damping moment coefficient distribution was
summed over the whole hull to yield the total vortex roll damping moment
coefficient The table below shows values of C ^ £ tor the
present DVM calculations (except the bilge keel calculations), for the
method due to Ikeda et al. (1978) - the "Ikeda" semi-empirical method,
and for the BMT method. Experimental results are given which have been
derived from both roll decrement tests and forced-roll tests. A range of
experimental values is shown, this representing the range of scatter in
the experiment. All damping results presented apply to the full-scale
situation.
105
Method Crdm2 ^tonnes 11,2/ rad2 )
Discrete Vortex Method:
Run 1 (bar keel) 22
Run 2 (bar keel) 18
Run 3 (bar keel) 21
Run 4 (bar keel) 17
Run 5 (no bar keel) 14
Ikeda method (no bar keel) 4
BMT method (no bar keel):
4° amplitude 13
8° amplitude 11
Experiment (bar keel):
Roll decrement tests 9 - 1 2
Forced-roll tests 10 - 13
The DVM computations show that both a lowering of the roll centre and a
lowering of roll amplitude from the values used in Run 1 reduce the
total vortex roll damping moment coefficient. A comparison between Runs
1 and 5 shows that the presence of the bar keel increases damping by
60%. Part of this increase can be accounted for by the fact that, in the
absence of a bar keel, computations are restricted, as described in
Section 5.5.3 above, to sections with sufficiently "sharp" shedding
geometries. This will produce the effect of somewhat underestimating

106
total vortex damping.
Interpretation of the relative accuracy of the various results is
difficult because of the uncertainty over the experimental damping and
because neither the Ikeda nor the BMT method include the effects of the
bar keel present in experiment.
Comparison between the Discrete Vortex Method and the other roll damping
prediction methods should be made for the DVM calculations carried out
in the absence of a bar keel (Run 5). Agreement is good with the BMT
method, but extremely poor with the Ikeda method, there being a
difference of a factor of nearly four between that method and the DVM.
One should expect the Discrete Vortex Method to be the most reliable
since it provides the most complete model of vortex shedding. There
appear to be weaknesses in the Ikeda method when employed for geometries
such as the naked hull Sulisker geometries which are substantially
different from those against which it was originally validated.
There are uncertainties in the experimental vortex roll damping results
for two basic reasons. In the first place, it is difficult to break
total roll damping down into its linear (i.e. radiation damping) and its
non-linear (i.e. vortex damping) portions. In the second place, there is
a significant amount of scatter in the experimental total damping
results. There are differences between the decrement and forced-roll
tests (forced-roll tests being believed to be more reliable - Spouge,
1987); also, experimental vortex roll damping moment coefficient is
averaged over a range of roll amplitudes, thus effectively ignoring any

107
variation with roll amplitude. If comparison between the DVM results is
made only with the forced-roll experimental results, then the DVM
overpredicts damping by between 30 and 120%. Of course, agreement might
even be better (or indeed worse) since there could be a greater
variation in roll centre than is allowed for in the DVM calculations
here; also, it would appear that DVM calculations at still lower roll
amplitudes might predict yet lower vortex roll damping moment
coefficients. Thus comparison with overall experimental dampings
produces no conclusive evidence as to precisely how accurate the DVM
predictions are. One should expect the DVM to overpredict
two-dimensional dampings by between 30 and 45% (see Chapter 4). However,
another possible source of error (inherent, in fact, in all three
prediction methods) is the use of a strip approach for calculating total
vortex dampings (see Section 3.9.5). The vortex shedding processes will
be affected by three-dimensionalities particularly near the bows and
stern but the magnitude of this effect is unclear. It is likely that the
use of strip theory will lead to an overprediction of damping.
A table of total vortex roll damping moment coefficient, Cr(j 2 » is
presented below for the Sulisker with the addition of the two bilge
keels.
108
Method S'dirtf ^tonnes 11,2/ rad2 )
Discrete Vortex Method:
Run 6 (bar keel and two bilge keels) 27
Run 7 (bar keel and two bilge keels) 22
Ikeda method (no bar keel, two bilge keels) 7
BMT method (no bar keel, two bilge keels) 18
Experiment (bar keel and two bilge keels):
Roll decrement tests 11
Forced-roll tests 14
Run 6, which is otherwise identical to Run 1 except for the two bilge
keels, shows that the DVM predicts an increment of 5 tonnes m2/ rad2 due
to the bilge keels. This increment does not vary significantly due to
small changes in roll amplitude or roll centre, as can be seen from Run
7, which is otherwise identical to Run 4 except again for the two bilge
keels.
It is not clear what the true experimental vortex damping increment due
to the bilge keels is. One can only conclude that the DVM prediction
appears feasible in the light of the vortex damping increments predicted
by the Ikeda and BMT methods.

109
CHAPTER 6
DEVELOPMENT OF A HYBRID MESH / POINT VORTEX NAVIER-STOKES SOLVER
FOR THE PREDICTION OF THE FLOW AROUND ISOLATED EDGES
The limitations of the Discrete Vortex Method so far as the present work
is concerned have been described in Section 1.3. These are primarily
twofold - the inability of the method to model either turbulence or
diffusion in the wake and the inability to predict the separation point,
which is problematical where that point is not fixed by geometry. In the
present research the problem arose in the calculation of vortex shedding
from the isolated edge with rounding, drag results from which were then
to be utilised via a matching process to predict the roll damping of
rectangular barges with rounded bilges. Further problems can be
envisaged in the prediction of the roll damping of more general ship
cross-sections in cases where there is no sharp edge to fix the
separation point. A method for these cases was clearly required which
solved the full two-dimensional unsteady Navier-Stokes equations. The
Introductory chapter, Chapter 1, explains the choice of the method,
which will be described fully in this chapter, and results from which
may be found in Chapter 7. This method will be designated here the
Hybrid Moving Vortex Diffusive Method, or HMVDM for short.

110
6.1 Fundamental Equations
This section describes the fundamental fluid flow equations, together
with their boundary conditions, which govern two-dimensional unsteady
fluid flow. The HMVDM employs a system of point vortices superimposed
onto a finite difference mesh to produce a numerical solution of the
fluid flow governed by these equations.
6.1.1 Primitive equations
The fundamental equations for unsteady two-dimensional incompressible
flow of a Newtonian fluid with no body forces and constant properties
are the Navier-Stokes momentum equations in the x and y directions and
the continuity equation (Roache, 1972).
6 .1.1
6 .1 . 2
6.1.3
Ill
6.1.2 Stream function and vorticity transport equations
By eliminating pressure from Equations 6.1.1 and 6.1.2 above, and
defining vorticity (ft) as:
ft = 2 1 9u
6.1.4
“ 9x 9y
we obtain the parabolic vorticity transport equation:
9ft ..9ft ..9ft . f92ft 92ft

■7TT = ~ U---- V-r— + V 6.1.5
9t 9x 9y 9x2 9y7
We define the stream function (ip) by:
dip
9 y u; 6 .1.6
so as to rewrite Equation 6.1.4 as an elliptic Poisson equation,
V 2ip= -ft 6.1.7
It is in the form of the vorticity transport equation that the
Navier-Stokes equations are solved in the present work. The vorticity
transport equation consists of the unsteady term, 9ft/9t, the convective
terms, u9ft/9x and v9ft/9y, and the viscous diffusion term,
v(92ft/9x2 + 92ft/9y2 ). It is parabolic in time, and therefore poses a
marching or initial-value problem, whereas the stream-function equation
is elliptic and poses a "jury" or boundary-value problem.

112
6.1.3 Boundary conditions
The above equations are solved subject to certain boundary conditions,
which will vary according to the problem concerned. The boundary
conditions in the solution of the vorticity transport equation are
conditions on vorticity. In the present problem the body side must be
treated with a no-slip condition such that the wall vorticity ensures
zero tangential velocity at the wall. At all the other computational
boundaries slip conditions should be employed such that ft = 0.
The boundary conditions in the solution of the Poisson equation are
conditions on the stream function. The two commonly used boundary
conditions in the literature are either the Dirichlet or the Neumann
(see Roache, 1970), with the Dirichlet typically setting the stream
function constant on the upper (outer) and wall boundaries and to some
inflow function at the upstream and downstream boundaries. Making the
outer boundary a streamline produces the same restriction on the flow as
a wind-tunnel wall. The Neumann boundary condition would be most
commonly used at the outer boundary and involves setting the derivative
of the stream function in the direction normal to the boundary to the
value of freestream velocity there.
6.2 Finite Difference Forms
As a fundamental prerequisite of the computational solution of the above
equations on a finite difference mesh, we require difference forms of
each of the terms.

113
Forms of the first and second order difference formulae are required
which apply to unequal mesh intervals. These can be derived using either
a parabolic curve fit or a Taylor series expansion. The latter method is
followed here since it enables the accuracy of the differencing scheme
to be ascertained. The following analysis takes after the work of
Salvadori and Baron (1961).
We define a function f(x) whose differences we wish to establish at a
point m, given f ,, f and f We define:

r J m-1 m m+1
Ax = x ,, - x ; Ax , = x -x s=Ax /Ax , 6 .2 . 1
m m+1 m* m-1 m m-l* m m-1
(see Figure 6.1). Then we know that:
f = f + f 1 Axm + i f" Ax2 + 1 f"'Ax3 + 6 .2 . 2

m+1 m m m 2 m m 6 m m
and:
: , = f - f' Ax , + i f" Ax2 , - i f"1 Ax3 i 6.2.3

m-1 m m m-1 7 m m-1 6 m m-1
Subtracting Equation 6.2.3 from Equation 6.2.2 we obtain
- f
■f1 = m+1
"iT» m~l ,1 îi 6.2.4
m ” Ax„ + Ax_ n 7 m
'm-1
't - + 0 (Axm)
m
Taking only the first term in the expansion, neglecting higher order
terms, yields an expression for f' which is at least first order
accurate, and second order accurate if s = 1. Eliminating f ‘ between

114
the two Equations 6.2.2 and 6.2.3 we obtain an expression for f"m :
Or- sf
fi, _ 2s m-1,-fm (l+s)+f
v ' m+1- ..
l_flI Ax 6.2.5
- 2 (l-s) + 0 (Ax )
m sTI Ax2 3 m s v ' v nr
m
Like the first derivative, the resultant expression for f" obtained by
taking the first term of the expansion is at least first order accurate.
For a regular mesh spacing, with s = 1, and taking only the first terms
of the expansions, Equation 6.2.4 becomes:
fm' =
fm+1 ' fm-1
.
6 2.6
2Ax
where Ax = Ax„ = Ax_ ,

m m-1
and Equation 6.2.5 becomes:
f , - 2f + f i
m-1_____ m m+1
f" = 6.2.7
m
Ax2
The above two difference forms are found very commonly in the literature
- see, for example, Roache (1972). They are both second order accurate.
6.3 Outline of the Basic Inviscid Cloud-In-Cell Method
The HMVDM computer program referred to above originated as an inviscid
"Cloud-In-Cell" vortex shedding program written by Basuki at Imperial

115
College; it is convenient to summarise the operation of this method at
this point. Although work on this particular program has never been
reported in detail, the method follows closely an earlier Cloud-In-Cell
method also due to Basuki (1983). The later method employs a
rectangular mesh in the transform plane, whereas the earlier method
transforms to a circle. Basuki's earlier computer program is also
related closely to a program reported by Naylor (1982).
The Cloud-In-Cel 1 codes represent modifications to the basic Discrete
Vortex Method described in Chapter 3. These are as follows:
A) At every time step point vorticity is distributed onto the mesh,
which has been set up at the beginning of the computation. The
distribution function used in the above work was a reverse bilinear
interpolation which distributed vorticity from a point vortex onto the
four nearest grid points.
B) At every time step the Poisson's equation (Equation 6.1.7) is solved
on the mesh using the known values of mesh vorticity to obtain mesh
values of streamfunction. Equation 6.1.7 may be solved on a mesh by one
of many fast computational methods; both Naylor (1982) and Basuki (1983)
chose the Fast Fourier Transform method.
C) Mesh velocities are computed from the streamfunction using a finite
difference approximation.
D) Convection velocities are computed at the point vortices from the

116
known mesh velocities using a bilinear interpolation scheme and the
vortices are convected according to their respective convection
velocities.
Thus the above steps replace the previously described influence function
method for the convection of vortices (see Equations 3.2.3 to 3.2.6).
6.4 Outline of the Hybrid Moving Vortex Diffusive Method
The essential features of the HMVDM are outlined here and then treated
in more detail in succeeding sections. The HMVDM represents a joint
development of the basic Cloud-In-Cell vorticity convection model by 0.
M. R. Graham and the present author to take account of the diffusion
term within the Navier-Stokes equations, yielding a method which solves
a finite difference approximation to the Navier-Stokes equations.
The method solves, subject to the elliptic Poisson equation, the
Navier-Stokes vorticity transport equation. Performing the computation
in the transform plane gives rise to a term additional to those found in
the basic vorticity transport equation (Equation 6.1.5), this term being
the square of the transformation modulus (see Lecointe and Piquet,
1986):
3ft fa2^ a2^] -

V
at ax ay az
117
The solution is performed computationally in two distinct parts, each
preceded by a solution of the Poisson Equation:
1) Diffusion
'd2Q n 82Q nl
Qk = n + At v 6.4.2
[lx5- 3y2 J 8z
where the superscript "n" refers to conditions at time t, and "k" to
conditions intermediate between t and t+At, with vorticity having been
advanced one time step only so far as the diffusion term is concerned.
Diffusion is carried out on the mesh using a conventional differencing
scheme.
2) Convection
n+1 BC
Q = ftk - At [ X + v30 1 6.4.3
3y J Bz
where the superscript "n+1" refers to conditions at time t+At.
Convection velocities are calculated on the mesh, but convection is
actually carried out using the point vortices, not the mesh.
This gives rise to the following solution of the vorticity transport
equation:
118
f32ftn a2n"
nn+1 = sin + At v At
3x ly^
r 3 f32fin 3 2fin ] 3 [32ftn 3 2d nl11 3r

- vAt2
_ 3x [ax2 ' a y 2 J+ v 3y [ 3 x 2 3y 2 1Jr 32 6.4.4
The solution is only first order in time, with a diffusion error being
convected at each time step. This error is in addition to the error
caused by forward differencing in time.
The HMVDM has the following stages:
1) Perform transformation.
2) (Stage A of the Cloud-in-Cel1 method described in the section above)
At every time step distribute point vorticity onto the mesh, which has
been set up at the beginning of the computation. The distribution
function used is a reverse bilinear interpolation which distributes
vorticity from a point vortex onto the four nearest grid points.
3) (Stage B of the CIC method described in the section above)
At every time step solve the Poisson's equation (Equation 6.1.7) on the
mesh using the known values of mesh vorticity to obtain mesh values of
streamfunction. The Fast Fourier Transform method (Cooley and Tukey,
1965) is employed in the solution of the Poisson Equation.

119
4) Perform diffusion of the vortices on the mesh subject to the
vorticity boundary conditions.
5) Solve the Poisson's equation (Equation 6.1.7) again as in Stage 3
above, using the values of mesh vorticity computed in Stage 4.
6) Recompute an improved estimate of the wall vorticity and diffuse the
excess wall vorticity out into the field.
7) (Stage C of the CIC method described in the section above)
Compute mesh velocities from streamfunction using a finite difference
approximation.
8) Transfer the extra mesh vorticity due to diffusion onto the point
vortex system.
9) (Stage D of the CIC method described in the section above)
Compute convection velocities at the point vortices from the known mesh
velocities using a bilinear interpolation scheme and convect the
vortices according to their respective convection velocities.
A flow chart describing the operation of the HMVDM is shown in Figure
..
6 2
120
6.5 Transformation
6.5.1 Utilisation of the general numerical Schwarz-Christoffel
transformation
The transformation has been described fully in Section 3.1. For the
present implementation the geometry input has been set up for an
isolated edge with rounding; it would be simple, however, to change the
input to any other reasonable configuration, and this is one of the
possible areas of future work.
Two quantities were required as output from the transformation: dz/dC in
terms of C, and z in terms of C. The derivative values on the body at
the junction between two segments are zero if the interior angle between
the two segments is less than 180°. To avoid an infinite reciprocal of
the derivative its value was calculated at a position 1/10 of the
distance to the next grid point out.
6.5.2 Setting up the mesh
The Schwarz-Christoffel transformation has a particular usefulness in
the setting up of computational meshes. A mesh composed of equal
rectangular elements in the transform plane transforms back to a mesh in
the physical plane which is distorted in such a way as to concentrate
mesh points in the regions where they are most required, i.e. in the
boundary layer, especially near to body slope and curvature

121
discontinuities. Actually, the mesh as described above in fact comprises
a system of streamlines and equipotentials for the uniform flow
situation in the absence of vortex shedding in both the transform and
physical planes, which explains the fact that the mesh is finest near
regions of significant flow acceleration.
For flows at high enough Reynolds number the experience of a number of
workers (see, for example, Peace and Riley, 1983, or Bayliss, Gunzberger
and Turkel, 1982) has shown that a simple mesh composed everywhere of
equal rectangular elements is inadequate. There are two related reasons
for this. Firstly, extra concentration of mesh points is required in and
near the boundary layer. For laminar flow, the higher the Reynolds
number, the thinner the boundary layer, all other aspects of the flow
being identical, and therefore the mesh is required to be more
concentrated near the surface. In fact, it is widely held that at least
three points are required to represent the boundary layer - one at the
surface, one at the edge of the boundary layer, and one actually in the
boundary layer, see Roache (1972). Secondly, the mesh is required to
extend out into the field a substantial distance compared with, say,
some representative flow length scale.
Both of the above requirements could in theory be satisfied by a
transform plane mesh composed of a sufficient number of equal
rectangular elements. The problem would be that for high Reynolds number
the total number of mesh points might be completely unaffordable
computationally speaking. This problem is addressed in the present work
by having a mesh spacing in the transform plane allowed to vary in the

122
direction normal to the surface (y-direction). It would also have been
advantageous, though not so important, to have a variable mesh spacing
in the streamwise (x-) direction. However, the requirements of the Fast
Fourier Transform only permitted regular spacing in this direction.
The y-mesh spacing in the transform plane was chosen to have an
exponential variation:
y, = b [ea *(J_1 * - 1] 6.5.1

J
A typical example of the use of the Schwarz-Christoffel transformation
for the generation of computational meshes is shown in Figure 6.3, with
the y-mesh spacing in the transform plane expanding exponentially
outwards.
An inverse is easily found for Equation 6.5.1:
j = INT [ log(y ./b + 1) / a ] + 1 6.5.2

J
An exponentially varying mesh has been used widely by other workers
(see, for example, Peace and Riley, 1983). For values of a.(j-l)
sufficiently low the mesh distribution is linear. This fact enables
comparison with uniform mesh methods. It would also be possible to use
other mesh functions. If an analytical inverse did not exist a root
finding process would be required which might be time consuming. An
exponential mesh variation in the transform plane will lead to an even
more significant concentration of points near the surface in the
physical plane because of the properties of the transformation mentioned
above.
123
6.6 Treatment of the Diffusion Term in the Solution of the Navier-Stokes
Vorticity Transport Equation
This section deals with the treatment of the diffusion part of the
Navier-Stokes vorticity transport equation, that is Equation 6.4.2
above. This mainly covers Stages 4 and 6 of the HMVDM.
6.6.1 Solution of the diffusion equation
This section refers to Stage 4 of the HMVDM
Equation 6.4.2 may be written in a finite difference form for variable
y-mesh spacing as follows:
n n
ft 2ft? . + ft
+ At v i-1 »j 1 >3 i+1 ,j +
2
AX-j
n
sftn . . ” •(1+s) + ft
2s
s+T *
i ,j-l ___ *>j_____ __ i»j+l
6 . 6.1
where Ay. = y.,, - y.; s = Ay ./Ay. ,

J J+l 3 3 J“1
In fact, the method works mostly in terms of circulation values, where
mesh circulation r. . is here defined as:

•»J
ri,j = • Axi • (Ayj + Ayj-i)/2

.
6 6.2
124
Then Equation 6.6.1 may be rewritten:
ri-1,j-2rL +r^ , j . 2s ^ J - l )
r^ . r1? . + Atv
15 J T >J Ax? ^ Ay?
1
f sr? . i
i ,j-i (1 + S >r i , j , r?,j+ l [ 9C
6.6.3
9z
[<A y j - l + A y j-2> <A y j + A y j- l > (A y j + l + A y j)J
For an exponentially varying mesh,
Ay.i+ Ay. H Ayj * i + Ayi 6.6.4
Ayj-1 + Ayj-2 Ayj + Ayj-1
so that the y-directi on diffusion term can be simplified to:
32fi 2s # s ri,j-l + s ri,j+l

6.6.5
9y:
5+1 Ay2
and this is the form in which the diffusion equation is solved for every
interior mesh point.

125
6.6.2 Vorticity boundary condition
This section refers to Stages 3 to 6 of the HMVDM
The diffusion equation is solved subject to the no-slip boundary
condition at the wall (see Section 6.1.3 above). Boundary conditions on
the other three boundaries are slip conditions with zero vorticity. A
second order finite difference approximation for the wall vorticity
boundary condition due to Woods (1954) is employed. This may be written
in its circulation form, for variable y-mesh spacing, as:
rn+l
pn+l 3( C 1 - ^
ywn+1)
J !w+l
1w AxAy, 6 .6 . 6
w
Ay,w Ax -(Ayw +iyw+i)
where the subscript "w" stands for conditions at the wall, "w+1" for
conditions one mesh point out in the y-direction, and the superscript
"n+1" stands for conditions at time t=n+l, that is one time step forward
from the current conditions. Vorticity is positive anticlockwise. Figure
6.4 provides a sketch which illustrates the above points.
Roache (1972) outlines the problems encountered when applying the
vorticity boundary condition. The most obvious one is that in the above
equation neither streamfunction nor mesh circulation are known initially
for the new time step (n+1). Mesh solutions of the vorticity transport
equation which are second order in time will usually adopt some
126
iterative procedure between this equation and the Poisson Equation so
that improved estimates of streamfunction and circulation may be
successively obtained. Typical of such a method is the Alternating
Direction Implicit Method used, for example, by Peace and Riley (1983).
In the HMVDM this iteration process is not feasible on economic grounds
and would be most difficult to implement computationally, and so an
alternative approach was used.
1) (Stage 3 of the HMVDM) At each new time step Poisson's Equation is
solved first, using values of r1? . , i.e. circulation from the previous
T »J
time step, but with the correct streamfunction boundary conditions for
kl
the new time (n+1). Obtain a first estimate of streamfunction, ip. ..
i jJ
2) (Stage 4 of the HMVDM) Obtain a first estimate of Tw by the method
described below. Call this estimate rw . In the subsequent diffusion
process an amount of circulation is diffused from the wall. It was noted
during numerical experiments with the method that if this exceeds a
critical amount instabilities develop. This is due to the fact that
there is a feedback mechanism in operation with the method which means
that an under- (or over-) specification of the boundary condition at one
time step leads to an over- (or under-) specification at the next.
Whether or not this feedback mechanism is stable or not depends largely
on the size of the error in the wall circulation. In order to ensure
that instabilities do not arise in this first estimate of wall
circulation, an adjustment factor, R, is applied to the wall circulation
formula (Equation 6.6.6) to yield the first estimate of the wall

rkl
circulation, , using the value of r from the old time step:
127
kl .kl,
,kl 3(iAw+1 -ip )
rw 1 n - R) C i
w
AxAy
w
6.6.7
Ay2
4x(iV AW
A suitable value of R has been found by numerical experimentation to be
2. Diffusion is then carried out throughout the mesh as described in
Section 6.6.1, using the first estimate of wall

3
circulation, rwkl , to
yield estimates of the mesh circulation elsewhere, r*kl . .

•>J
3) (Stage 5 of the HMVDM) This is followed by a second solution of
Poisson's Equation, with these revised values of mesh circulation, to
yield ^ ..
i>j
4) (Stage 6 of the HMVDM) A second estimate of wall circulation is
obtained as follows:
,kl
3 ( ^ +r ^ )
r wk w+1
AxAy,
w
6 .6 . 8
Ayw Ax(Ayw+Ayw+i )
that is, the correct boundary condition is applied on this second pass.
5) (Stage 6 of the HMVDM) The new wall circulation is employed to

b
determine via the diffusion process, but further estimates
y
elsewhere in the mesh are not carried out, so that n • is elsewhere
*>J
X J. Tiki
set to T . .*
•»J
Thus the above scheme makes two estimates of both streamfunction and
128
wall circulation per time cycle. Two estimates are also made of the
circulation one mesh point out from the wall, which is crucially
important since this is where new discrete vortices are introduced into
the flow. Note that the above scheme only advances circulation and
streamfunction to the intermediate stage (k), rather than the new time
step (n+1). The new time step is achieved through convection of the
circulation.
Roache (1972) mentions the existence of instabilities encountered in the
use of the second order wall vorticity boundary condition. In fact
numerical experimentation with the current method has established a
clear link between these and time step size. It has also established
that the same instabilities can exist even for the first order method.
The strength of the current scheme lies in its exploitation of the
feedback mechanism at the wall to enable a much more stable scheme than
would normally be formulated and thus to enable for little extra
computational cost larger time steps to be used.
6.7 Solution of the Poisson Equation
This section refers to Stages 3 and 5 of the HMVDM
Solution of the Poisson Equation (Equation 6.1.7) is carried out using a
Fast Fourier Transform technique. Basuki (1983) has described the
application of the technique to the Cloud-in-Cell method with a circular
mesh; here we describe its use for rectangular meshes, with appropriate
extensions for variable y-mesh spacing.

129
The Poisson Equation may be written in a discretised form as:
-2T. .
V zip 6.7.1
AxCAy^.+Ay^^)
where Ay. and A y . , are as defined in Section 6.2. Taking the Fourier
0 J 1
Transform in the x-direction, we obtain:
-2G
k jj
6.7.2
-k2 \ . j + dy2
k,j AxfAyj+Ayj.-,)
where 'P='i/(k,j) is the Fourier Transform of ^(i,j)
i -i kx ,
ip.e dx,
k is the wave number, and G is the Fourier Transform of circulation,
We can write a finite difference form of the left hand side of Equation
6.7.2 above using the finite difference approximation derived in
Equation 6.2.5:
J r r ’1'k,j+l - ( - k 2A y ? - 2s) T k i j + ! £ • Tk
2Gk,.i Ay.i 6.7.3

ix(Ayj+Ay._1)
130
Notice that as always y. is invariant of its x-location. This equation

J
when applied successively to each interior mesh point forms a series of
tridiagonal matrices which can be solved by a Gaussian elimination
technique. In this work the elimination technique described by Varga
(1962) is employed, with:
-2Gk , ^ j
vj ^kj+l + V ^ k , j + ^j^k,j-1 6.7.4
Ax(Ay .+Ay ._.j)
where v. = 2s_ , y. = 2s^_ , X- = -(k2Ay^ + 2s)

J s+1 J s+1 J ^
As indicated in the previous section, solution of the Poisson Equation
is performed twice per time step, once to provide values of stream
function for the wall vorticity boundary condition, and once to provide
values of stream function for the evaluation of mesh velocities.
6.7.1 Stream function boundary condition
The Dirichlet boundary condition with stream function set to zero at the
inner boundary and the product of freestream velocity and y-distance on
the outer boundary (see Figure 6.5), was adopted due to its ease of use
with the Fast Fourier Transform method. Due to the periodic nature of
the FFT, it is not possible to specify upstream (inlet) and downstream
(outlet) boundary conditions. The FFT automatically sets streamfunction

131
values at the inlet equal to those at the outlet (i.e. a periodic
boundary condition). This imposes the restriction on the HMVDM that any
significant transport of vorticity should be performed well clear of the
inlet and outlet boundaries. The use of the Dirichlet boundary condition
on the upper boundary requires that this boundary should be sufficiently
far from regions of significant vorticity that its effect on the vortex
shedding processes is negligible. This can be verified by testing the
sensitivity of results to outward mesh extent.
6.8 Determination of Mesh Velocities
We can determine the mesh convection velocities from the values of
streamfunction computed on the mesh from the Poisson Equation. If the x-
and y- direction convection velocities in the transform plane are u and
v respectively, then we can write from the definition of streamfunction
(Equation 6.1.6):
6. 8.1
and:
.
6 8.2
The 3C/3z term appears because of the transformation. It has been shown
to be of the above form by Basuki (1983). Using Equations 6.2.4 and

132
6.2.6 we can write finite difference forms of the above equations as:
u. . 6.8.3
i »J
Ayj + Ayj-i
and
2
_ -1, j " î+1,j 6.8.4
Vi »J dz
2Ax
6.9 Convection of Vortices
The point vortices are convected at each time step in the transform
plane using the known mesh velocities as calculated above. An
interpolation scheme is required to determine the contribution of the
mesh velocities at each point. Generally, such an interpolation scheme
might pertain either to the nearest 1 , 4 or 9 mesh points. Here we
employ a bilinear interpolation scheme which employs the nearest four
points to give an expression for the u and v components of the vortex
convection velocity in terms of the u and v components of the mesh
convection velocities:
4 A
uvortex ~L?-j UL* X
and
133
V L
vvortex " VL * X
For an explanation of the above equations see Figure 6.6.
6.10 Interpolation Between the Point Vortices and the Mesh
At the beginning of each time step it is required to redistribute the
circulation on the vortices onto the mesh so that the diffusion process
can be carried out on the mesh. An inverse procedure to that of
Equations 6.9.1 and 6.9.2 is employed:
A,
vortex
6 . 10.1
for L = 1 to 4, where A^ are the areas as shown in Figure 6.7, Ay is the
area of the relevant mesh box, and Tvortex is the strength of the point
vortex.
6.11 Transfer of Mesh Circulation onto the Point Vortex System
The diffusion process as described in Section 6.6 above amends the
circulation values by an incremental amount at every mesh point at every
time step. Some method of efficiently transferring this incremental mesh
circulation onto the moving vortex system needs to be provided which
minimises the number of new vortices created at each time step so that
134
computational costs do not become excessive. The most expensive scheme
available would be simply to place into the flow at each time step
vortices whose initial location is at the mesh vertices themselves. A
more efficient method is to amalgamate this new circulation onto
pre-existent point vortices; this method is employed in the current
work. However, when there are no point vortices near enough to the
relevant mesh vertices, no amalgamation takes place, in order that
adequate definition of the vortex structure may be retained. In these
circumstances vortices are released directly into the flow from the mesh
vertices. The above scheme is described in more detail in the following
section.
6.11.1 Amalgamation of vortices
A first order vortex amalgamation scheme has been devised as follows.
Each mesh box is divided equally into four smaller boxes (see Figure
6.8a). A domain of four hyperbolic sections is constructed as shown in
the figure. A sweep is carried out through all the current point
vortices. During this sweep each point vortex is associated with its
nearest grid point, and the mesh vortex there (i.e. the new circulation
due to diffusion) is amalgamated with it if the point vortex lies inside
the domain and if the two circulations share the same sign. The extent
of the inclusion zone in this first sweep is indicated in the figure. In
a second sweep through the vortices each mesh vortex is amalgamated with
its nearest like-signed point vortex provided the point vortex lies
within the region shown in Figure 6.8b. Amalgamation of like-signed
vortices ensured that the total number of vortices in the flow was kept
135
within reasonable bounds. Vortices of opposite sign were not
amalgamated, in order to avoid disturbing the small scale structure of
the vortex sheet.
Vortices are amalgamated in such a way that the strength and position of
the resultant vortex is as given in Section 3.4.
6.12 Calculation of Vortex Force
The force due to vortex shedding is determined by means of a momentum
balance formed from Blasius' theorem in the way described in Chapter 3
for the Discrete Vortex Method.
It is necessary to show that this theorem is applicable to viscous as
well as inviscid flows. Wu (1981) has shown, for two-dimensional viscous
flows, that the aerodynamic force, F, exerted by a fluid on solid bodies
immersed in and moving relative to the fluid is equal to the inertia
force due to the mass of fluid displaced by the solid bodies (F ) plus a
term proportional to the rate of change of total first moment of the
vorticity field in the solid bodies and the fluid (i.e. the "vortex"
force, F , which, in viscous flows, is here defined as including the
contribution of the boundary layer to the force summation). Thus an
expression can be obtained for Fv :
da
Fv = - Pgj 6-12-1
136
where a is the first moment of the vorticity field, defined by:
a r x ft dR 6 .1 2 . 2
where r is a position vector, and R^ is a circular region r < L. L is
sufficiently large so that contains all the solid regions, and that
the total vorticity and total vorticity moment outside R^ are
negligible, ft is the vorticity field in the infinite region R.
As applied to the problem of vortex shedding off isolated edges, for a
transform plane with the freestream acting in the y-direction (parallel
to the body side) the above Equations 6.12.1 and 6.12.2 yield the
following expression for vortex force:
F xftdxdy 6.12.3
v
Converting vorticity in this equation into circulation and writing the
equation in a discretised form, we obtain for a set of n point vortices:
F x b.12.4
v V ' m m
m=l
which is identical to Equation 3.8.1 except for the use here of the real
part (xm ) of the transform plane point vortex complex ordinate. The real
part of the position of its image is simply (-x ). Equation 6.12.4 may
be rewritten in its mesh form by recognising that the mesh ordinates do
not change with time, but instead that the mesh circulation does:
137
NX MY
= ~ iPI I x — r 6.12.5
i SJ 3t i,j
i=l j=l
Equations 6.12.4 and 6.12.5 not only furnish proof that Blasius' theorem
in its present formulation can be used for viscous flows, but also
confirm the validity of the particular implementation used by Graham
(1980) and adopted for this present work. In the next chapter we examine
more closely the use of the formula for one particular test case, namely
the Stokes boundary layer (oscillatory boundary layer over a flat
plate), and then examine the application of the HMVDM to the flow past
isolated edges.
138
CHAPTER 7
RESULTS FROM THE HYBRID MOVING VORTEX DIFFUSIVE METHOD
7.1 Test Case of an Oscillatory Boundary Layer over a Flat Plate
A lengthy validation process was embarked upon to check the operation of
the Hybrid Moving Vortex Diffusive Method. This included testing of the
various subroutines in the computer method, especially those connected
with the Fast Fourier Transform Poisson Solver. In the next stage of the
validation, various practical flow cases of increasing complexity were
computed. These culminated in computations on the oscillatory boundary
layer over a flat plate, otherwise known as the Stokes boundary layer.
The oscillatory boundary layer was thought to be a suitable test case
for the method for two reasons:
1) It has direct relevance to the more general problem of vortex
shedding from isolated edges with edge rounding in oscillatory flow.
Indeed, it may be thought of as representing the case of infinite edge
radius.
2) There is an analytical solution available so that comparison can be
made with a known solution.
However, the Stokes boundary layer is not a good test of the convective
part of the algorithm, although the HMVDM, originating as it does from

139
the DVM, was believed to have a good algorithm for representing
convection.
7.1.1 Analytical solution for an oscillatory boundary layer
Schlichting (1960) states the analytical solution for the flow past an
infinite oscillating flat plate. This solution can be adapted to that
for a fixed wall and freestream of velocity:
= Uq sin(cot) 7.1.1
For the above freestream the analytical solution of the oscillatory
boundary layer is:
u / U q = - [e ^y .sin(oat-ky) - sin(o)t)] 7.1.2
where:
k = / ( oj / 2 v )
Now shear stress,
t = y3u/3y 7.1.3
= yl^k e ^ [sin(a)t-ky)+cos(cot-ky)]
t = /2y Ugk e ^ sin(wt-ky + tt/4) 7.1.4

140
Wall shear stress,
xw = /2 y U^k sin(o)t+ tt/ 4 ) 7.1.5
In other words, wall shear stress is tt/4 out of phase with the
freestream velocity. Force per unit width on the plate is given by:
Fv= tw . (domain length) 7.1.6
Vorticity can be simply derived from the above analysis by recognising
that there will be no changes in the v-velocity in the x-direction (in
fact, v = 0 throughout the flow). Then:
= - 3u/9y 7.1.7
Q = _ /2Uq k e ^ sin(wt-ky + tt/4) 7.1.6
It has been shown in the previous chapter that the Blasius momentum
balance method of calculating drag is applicable to viscous as well as
inviscid flows. In the case of viscous flows the momentum balance
includes the boundary layer contribution to drag as well as the
contribution due to vortex shedding. This can be demonstrated for the
specific case of the oscillatory boundary layer, where the only
component of drag is boundary layer (or skin friction). In terms of wall
shear stress as opposed to force, Equation 6.12.4 gives us for a
freestream parallel to the x-direction that:

141
d v
T,, = - p T-r ) Q _ Ay. ym 7.1.9
w 8t L m ‘'m m
m
Writing this in the continuum form,
✓oo
9v 9u "
Tw " p 9t 9x " 9y, •y dy
'0
3v
But 8x = 0 S0 that:
9u .
V P 9t y 9y dy
Integrating by parts,
u dy 7.1.10
Tw = " p 9t
The full Navier-Stokes equation of motion for the u-velocity is:
9u , 9u ^ 9u 1 9p , 92u 92u
7.1.11
9t u 9x v 9y ■p 9 7 + V 9x2 W
with:
2E = v = 0
9x 9x
So that:
/• o o 0
92u
T = -pV dy
W
142
^ _ au
w y 8y 7.1.12
w
which is the expression for the shear stress to be found in Equation
7.1.3. Thus the method used in this work for calculating vortex force
(Equations 6.12.4 or 6.12.5) is shown to include and correctly predict,
for this particular viscous flow case, the boundary layer (skin
friction) force.
7.1.2 Comparison of the analytical solution with the computed results
Results of two test computations are presented here, together with the
appropriate analytical results. The first was carried out on a 32 by 32
mesh with 80 time steps per cycle at a kinematic viscosity (v) of 0.1,
freestream velocity amplitude (Uq ) of 1.0, and angular frequency (oj) of
0.277. In the computation the y-mesh spacing was chosen to expand
outwards with each mesh spacing being 1.051 times the previous. This
represents quite a moderately expanding mesh. The x-mesh spacing was
selected to be large enough so as to produce a computation where, far
from the edges of the computational domain, changes in flow quantities
in the x-direction were minimal. The mesh is shown in Figure 7.1a. The
above aim was achieved, as Figure 7.1b shows. In Figure 7.1b velocity
vectors have been plotted for every y-wise mesh location and every
fourth x-wise location across the entire computational domain. The
vectors represent the instantaneous situation at the end of two cycles
of flow. The flow at the edges, as would be expected, is different from
that nearer the middle of the domain, but one can see that there is a
143
substantial region where flow changes in the x-direction are minimal.
Fiaure 7.1c shows an enlarged plot of the y-wise variation of velocity
at a station near the middle of the domain at the end of two time
cycles, that is, when freestream velocity is zero. The full line
represents the analytical result; the symbols represent the calculated
result. The analytical result has been calculated with a phase shift of
-0.011 radians with respect to the computed, so as to compensate for a
small error in computed freestream velocity by matching the two
freestream velocities together. This phase shift represents 0.14 of one
computational time step, which is considered low for a computer method
which is only first order accurate in time. The level of agreement
between the two velocity profiles is excellent.
Figure 7.2 shows a comparison between analytical and computed
time-dependent wall shear stress for two cycles of flow beginning at the
end of the first half-cycle. The computation was carried out with a
somewhat more refined mesh and 320 time steps per cycle. The agreement
between the two results is clearly very good. The main criticism of the
computation is that it is up to about 2% of a time cycle out of phase
with the analytical result.
On the evidence presented above it would appear that the HMVDM copes
most effectively with a diffusion dominated flow such as the oscillatory
boundary layer. This is encouraging, since the strength of the HMVDM is
in its treatment of convection dominated problems, and any weakness that
might exist in the method is likely to lie in its treatment of
diffusion.
144
7.2 Vortex Shedding from an Isolated Square Edge in Uniform Flow
The case of uniform flow past an isolated square edge has previously
been computed using the Discrete Vortex Method both by the present
author and also by Pu11in (1978). Because of the possibility of
comparison with DVM results it therefore affords a good test of the
HMVDM.
7.2.1 Definition of Reynolds number in uniform flow
One feature of the plain isolated edge flows being considered in this
thesis is that they possess no physical length on which scaling may be
carried out, though such a scaling is possible in the transform plane.
Equation 3.1.15 gives an expression for the length scale, l_z , of the
flow in the physical plane based on the properties of the
transformation. Such a length scale can be used in the definition of a
flow Reynolds number, but it should be realised that such a Reynolds
number is not based on the size of the vortex structure, but rather on a
length scale which may be significantly larger than, although
proportional to, the vortex size. Consequently, Reynolds numbers may
appear to be quite high. A Reynolds number based on l_z would be:
Re = Lz Vz /v 7.2.1
145
where Vz and l_2 are defined as for oscillatory flow (Equations 3.1.13
and 3.1.15) save that here V is the uniform flow velocity in the
transform plane and T is the time since the beginning of the
computation.
7.2.2 Vortex positions and streamlines
The HMVDM was run for the isolated square edge in uniform flow using a
64 X 32 mesh with a mesh spacing which produced a high degree of
refinement near the shedding edge, as Figure 7.3 shows. Figure 7.4 is a
plot of the streamline pattern and vortex positions produced after 100
time steps, together with results calculated by Pullin (1978) using the
Discrete Vortex Method (an essentially inviscid model). For the HMVDM
run Reynolds number based on length scale was 3100. The boundary layer
ahead of the vortex cluster was permitted to grow from the surface at a
position five mesh points from the left hand (inlet) boundary. Flow
velocity at the inlet and outlet boundaries was set to attached flow
values. It was ensured that vortex extent was small in comparison with
mesh size. In the plot showing vortex positions the cross symbols
represent positive circulation, and the triangles negative circulation.
Symbol size is representative of the value of circulation.
The comparison between the two results is extremely good, the slightly
flatter appearance of the vortex in the Navier-Stokes solution arising,
it is believed, from the action of viscosity in the HMVDM, although its
effects are unlikely to be too pronounced in this case.

146
7.3 Results for Oscillatory Flow Past a Square Edge with Rounding
A series of computer runs were carried out using the HMVDM to predict
the vortex shedding due to sinusoidally oscillating flow past an
isolated square edge with various degrees of edge rounding. Results
presented in this section thus echo the DVM computations of Section 4.3;
the HMVDM is shown here to produce results of superior quality and
reliability. The ultimate aim of the computations was to produce drag
coefficient data which could be utilised via a matching process in the
prediction of the roll damping of rectangular cross-section barges with
rounded bilges. Section 7.4 treats this subject in detail; this section
deals with the computations and the direct application of their results.
7.3.1 Definition of Stokes number in oscillatory flow
In oscillatory flow we employ the viscous parameter, Stokes number (St),
otherwise known as the 3 parameter (Bearman et al., 1985), which is
defined as:
St = Re / Kc 7.3.1
where both Reynolds number and Keulegan-Carpenter number are here
defined in terms of flow length and velocity scales. Thus we obtain
that:
St = l
J i = Lz2 / ( vT) 7.3.2
v
147
One of the most important parameters in the computations was the Stokes
number. Ideally, one would require results at Stokes numbers
representative of both model and full scale flows. A typical model scale
Stokes number based on flow length scale was of the order of five
thousand for the rolling model barge described in Chapter 2, whilst a
typical full scale Stokes number might be of the order of five million.
It is not known in what Stokes number range transition to turbulence
takes place, but it will depend on a combination of Stokes and
Keulegan-Carpenter numbers; at higher Kc transition will occur at lower
St. It is presumed that at full scale the effects of turbulence are
important. The higher Stokes number is therefore unachievable with the
HMVDM in its present form because it lacks a turbulence model. The lower
Stokes number is, however, achievable. In the computations which follow,
all "high" Stokes number runs were carried out at Stokes numbers
somewhat in excess of ten thousand. Two "low" Stokes number runs at a
Stokes number of 4600 are reported, to demonstrate the effect of varying
Stokes number.
7.3.2 Choice of computational mesh
Correct choice of mesh was most important and was complicated by two
conflicting considerations. On the one hand adequate definition of mesh
in and near the boundary layer was required. On the other hand the
complete mesh was required to be as large as possible to maximise the
possible path lengths of the vortex pairs and to ensure that the
boundaries were far enough from the body and the vortex clusters for the
148
boundary conditions to be correctly imposed. Because of the
computational expense it was not possible to employ meshes which were
large enough to enable computations to be continued for more than about
two cycles of flow without the outer boundaries imposing a restraint,
since vortex pairs would tend, over a number of cycles, to convect
towards these boundaries.
Four meshes were constructed; these are designated meshes A to D, and
may be found plotted in Figures 7.5 to 7.8. Meshes A and B are plotted
to nearly the same scale, whereas mesh C is about one quarter of the
size of mesh A and mesh D one seventh of the size.
7.3.3 Choice of run parameters
Duration of computations
One critical parameter is the duration of the computations in terms of
the number of time cycles. Ideally, one should want computations to be
carried out over a large number of cycles, for example, ten as is the
case in the Discrete Vortex Method calculations. Then it would be
possible to average drag over many cycles, excluding at least the first
cycle, when drag will invariably be lower. However, with the HMVDM it is
not possible to average over several cycles, since the cost of computing
more than two or so cycles of flow would be prohibitive. Further, the
extent of the mesh is not sufficient to permit longer runs. A more
extensive mesh would require more mesh points which would entail further
computational expense.
149
Thus, a strategy alternative to obtaining drag coefficients averaged
over several time cycles had to be adopted in the present work. It was
possible to perform integrals of drag over half-cycle, instead of
whole-cycle periods, by replacing Equation 3.8.3 by the following
equation for the A "half-coefficient":
ri2
3tt r • 27ft ft]
CFv b,n T d
t \0 W
and similarly for the B half-coefficient.
This approach enabled a fuller assessment of the development of drag
coefficient with time to be made (see Section 7.3.6). It was
particularly important to exclude those coefficients from the average
which appeared to have been affected by the starting conditions.
Experience with the method showed that it was sufficient to exclude the
coefficients calculated during the first cycle. It was therefore
possible to limit the important computational runs to between 2 and 2\
cycles, with drag coefficient usually being averaged over the second
time cycle only.
Non-dimensional edge radius
Edge radius (r) itself is a redundant parameter in the set-up of the
computations and is always set to unity. The important parameter in this
context is non-dimensional radius (r/l_z ), which determines the effective
edge radius in terms of flow length scale. Computer runs were carried
out at non-dimensional radii of 0, 0.0165, 0.0466 and 0.0924, in other

150
words practically the same range as the equivalent DVM runs (see Chapter
4).
X-wise mesh interval
Solution of the Poisson equation by means of the Fast Fourier Transform
permitted variable mesh spacing only in the y-direction, and not in the
x-direction. The x-wise mesh interval in the transform plane is
determined by two factors: by how many surface mesh points are required
in the rounding portion of the body, and by what the total mesh extent
in the x-wise direction is required to be. These are conflicting
requirements, in that as much definition in the rounding portion as
possible is required, and this is detrimental to the total mesh extent,
which is always required to be as great as is practicably possible. If
such a conflict exists, the only way of resolving it is to increase the
total number of mesh points in the x-direction. In the present work 64
points were used in the x-direction and 32 in the y-direction. It was
found that provided the calculations were not continued past about 2
time cycles the mesh was sufficiently large that vortex pairs did not
convect too near the inlet or outlet boundaries. The number of mesh
points within the edge rounding itself varied between 5 and 15 depending
on the flow length scale. For some test cases, several runs were
performed at different mesh spacings to ensure that an optimum mesh
spacing was being employed.
Y-wise mesh interval
Because of the facility in the program of varying the y-mesh spacing in
the transform plane, it was possible to crowd points near the surface,
151
as described earlier. The mesh was chosen to vary so that in the
transform plane the mesh interval increased at each station out by about
10%. Because of the nature of the transformation, this implies an
increase of somewhat more than 10% in the physical plane. The aim in
choosing the spacing itself (or, in other words, the value of "b" in
Equation 6.5.1) was to place an adequate number of points in the
boundary layer. The spacing was varied in order to find its optimum
value for each case.
Freestream velocity in the transform plane
Freestream velocity amplitude in the transform plane was varied to yield
the desired value of flow length scale and hence the desired value of
non-dimensional radius (r/l_z ).
Time step size
Time step size, At, was optimised by a process of numerical
experimentation. Too large a time step size leads to instabilities in
the solution and inaccuracies in the vortex representation. Too small a
time step leads to uneconomic computations. The same time step size was
employed for all the higher Stokes number calculations, leading to 1200
time steps per cycle. For the lower Stokes number calculation with
coarse mesh a much larger time step was possible, which led to 200 time
steps per cycle. In all cases the period of the oscillation, T = 10.
152
7.3.4 Flow results
The table below provides a list of the major calculations (except test
runs) which were carried out using the HMVDM, with corresponding figure
numbers for plots of vortex position and, in some cases, plots of
velocity vectors.
Run Mesh r / L z At St Figure numbers
(1) A 0.0 0.00833 28300 7.9,7.10
(2) A 0.0 0.00833 18400 7.11,7.12
(3) B 0.0165 0.00833 36700 7.13,7.14
(4) C 0.0466 0.00833 36700 7.15,7.16
(5) D 0.0924 0.00833 36700 7.17,7.18
(6) C 0.0466 0.00833 4600 7.19
(7) B 0.0466 0.05000 4600 7.20
In the plots of vortex positions it should be noted that the occurrence
of large concentrations of symbols does not necessarily imply the
presence of large concentrations of vorticity, since the process of
diffusion ma>' reduce the strength of the vortex cluster without
substantially altering its plotted appearance. Also, the size of the
vortex symbols plotted varies as the cube root of the vortex strength,
to enhance the visibility of areas of weak vorticity. There will be
occasions when many weak vortices in a plot will partially obscure the
153
presence of stronger areas of vorticity. For the above two reasons,
plots of vortex positions can not be seen as analogous to experimental
flow visualisations with either particles or dye. In this respect the
plots of velocity vectors give a more complete picture of the flow
situation. These are only presented here for cases where freestream
velocity is zero, as is indicated by the fact that many regions of the
velocity vector plots contain no plotted vector symbols.
RUN 1
Figures 7.9 and 7.10 show the first two time cycles of a square edge
calculation (Run 1). Each position plot is separated by a quarter of a
time cycle, and each vector plot by half a time cycle. Experimental flow
patterns for the same type of configuration are shown in Figure 2.4.
Experimental evidence would suggest that, after an initial period,
usually only one strong pair of vortices would be shed per cycle,
together with one weaker pair (the "residual vorticity"). In these plots
two pairs of nearly equal strength are being shed per cycle (ignoring
the first half cycle), a situation more reminiscent of shedding off
rounded edges (see Chapter 2). However, this is quite feasible. Two
pairs of nearly equal strength have been observed in experiment; also,
in the present calculation it is possible that if the computation were
extended a preferred shedding direction might eventually be established.
Both of the first two pairs which have been shed approach the body side
due to the greater vortex strength in each case of the vortex nearer the
body side. The first of the two pairs remains intact even by Figure 7.9h
and produces substantial secondary separations (Figures 7.10c and d),

154
with the two clusters in the pair separating under the influence of
their own images. On the other hand, the second vortex of the second of
the two pairs is drawn back into the main vortex shedding pattern to
form a part of the first vortex in the third pair. The rather intricate
pairing behaviour is typical also of experimental observations.
RUN 2
Run 2 was carried out at a slightly lower Stokes number than Run 1 but
with exactly the same mesh. The lower Stokes number is achieved through
a lowering of the flow length scale. Thus we should expect smaller
vortex structures to appear; the effect of the lower Stokes number
itself should be minimal. Figures 7.11 and 7.12 contain plots of vortex
position and velocity vectors for this case. Flow patterns are similar
to those observed in Run 1, but differences appear noticeably at the end
of i\ time cycles (Figures 7.Ilf and 7.12c). The second vortex to form
in the first pair collides with the body so near the shedding edge that
it immediately becomes drawn back into the shedding process and combines
with a portion of the second part of the second pair to form a third
weak pairing which convects away almost along the edge bisector. This
pair is so weak that by the end of the second cycle (Figure 7.12d) its
presence is scarcely noticeable, although from the position plot (Figure
7.11h) it might have appeared that its strength was quite considerable.
However, the third pair is strong enough also to interfere with the next
vortex pair to be shed by releasing some of its vorticity to be reformed
into part of this pair (the pair to the left of the shedding edge in the
above two figures). In Figure 7.H i yet another pair is seen forming,
but the final figure (7.11j ) might seem to suggest that both these last
155
two pairs have reformed into one larger pairing. Figure 7.12e confirms
the presence of a strong vortex pair at the edge.
RUN 3
The flow results of Run 3 (Figures 7.13 and 7.14) are fairly similar to
those of the previous runs despite the presence here of a small amount
of edge rounding. The double shedding pattern is still observed; on this
occasion the two vortex clusters in each pair are of more equal strength
so that the pairs do not convect towards the body side but out into the
field. Convection occurs at such a rate that by Figure 7.13g the first
pair has reached the outward extent of the computational mesh. The
general structure of the vortex clusters is very similar to those for
the plain square edge - compare, for instance, Figures 7.13c and d with
7.9c and d.
RUN 4
The effect of a further increase in edge rounding may be seen in Figures
7.15 and 7.16 (Run 4). The nature of shedding of vorticity has begun to
change. Due to the effects of edge rounding, vorticity is shed from the
edge in small "packets" onto the main vortex core (Figures 7.15b and
7.15d, for example). This form of vortex structure was also noted in
calculations performed with the HMVDM for uniform flow past rounded
edges, and is also observed experimentally. Figure 2.8 shows flow
visualisation pictures of the rolling barge at roughly the same values
of (r/l_z ) and Stokes number. Photos 3, 7 and 8 show clearly the same
types of shedding behaviour noted for Figures 7.15b and 7.15d. A rough
correspondence between the computed flow results and the experimental

156
flow visualisations of Figure 2.8 would be:
Figure 7.15b Photo 3
7.15c 5/6
7.15d 7/8
7.15e 1/2
7.15f 3/4
7.15g 5/6
7.15h 7/8
Considerable similarity between computation and experiment exists. The
vortex pairing is much less well structured than for Run 3. The first
pair to develop (Figures 7.15d and 7.16b) quickly collides with the
wall, the right hand member being almost completely recombined into the
next pair (Figure 7.15e), while the left hand member creates a secondary
separation at the wall, vorticity from which is then convected away into
the field (Figure 7.16c). Figures 7.15e and f, together with Figure
7.16c, show the second main pairing convecting away in two separate
parts from the right of the shedding edge, whilst the last figures
indicate another pairing convecting to the left and colliding with the
wall. The seemingly confused state of Figure 7.15h is actually
deceptive, as can be seen from the equivalent velocity vector plot.
RUN 5
Run 5 (r/Lz = 0.0924) has its approximate counterpart in the
experimental flow visualisation pictures of Figure 2.7. A rough

157
correspondence between the figures might be as follows:
Figure 7.17a ... Photo 7
" 7.17b ... " 8
" 7.17c " 2/3
" 7.17d ... " 4
" 7.17e ... " 5/6
" 7.17f ... " 8/1
11 7.17g ... " 2/3
" 7.17h " 4
As is the case elsewhere, the location of the regions of concentration
of vorticity should be ascertained with reference to plots of velocity
vectors (Figures 7.18). Because of the large edge radius, vortex
shedding is weak and less structured than at lower (r/l_z ). Two vortex
pairs per cycle are shed, although the second vortex of a pair appears
as little more than a thickening of the boundary layer * see for example
Figure 7.17c.
RUNS 6 AND 7
Runs 6 and 7 (Figures 7.19 and 7.20) are both carried out at a lower
Stokes number than the previous runs; the degree of edge rounding is the
same as that for Run 4. Both figures present the first time cycle of
flow. Run 6 is carried out on an identical mesh to that for Run 4 and
has the obvious difference of thicker boundary layer. Run 7 is a test of
sensitivity to mesh and time step size and is carried out on a

158
significantly larger, and therefore coarser, mesh than Run 6. Despite
this, the forms of the shedding in both Run 6 and Run 7 are very
similar.
7.3.5 Time dependent vortex force coefficient
Time dependent vortex force coefficient (Cpv ) is plotted against time in
cycles (t/T) in Figures 7.21 and 7.22.
Figure 7.21 shows a comparison between Runs 2 to 5, that is, between
different values of (r/l_z ). The general form of all the force traces is
similar. Increasing (r/l_z ) has the effect of reducing the force levels,
but introducing into the traces "wiggles" of increasing relative size.
The origin of the wiggles is thought to lie in the discretisation of the
vortex sheet into point vortices and to be the unavoidable consequence
of their occasional strong mutual interaction, especially near the
surface. Such an event happens once every 100 time steps or so and has a
negligible effect on force averages.
Force coefficient plots from Runs 4, 6 and 7 are plotted in Figure 7.22.
A comparison between Runs 4 and 6 reveals the differences due to an
almost tenfold reduction in Stokes number; a comparison between Runs 6
and 7 demonstrates the sensitivity of the HMVDM to changes in mesh and
time step size. A more detailed discussion of these differences may be
found in the next section. All the plots are of fairly similar
appearance, the two fine mesh plots (Runs 4 and 6) showing much more
detail. An instability in the computation can be seen starting to grow

159
at the end of the trace for Run 6. This prevented the run being extended
further. Apart from this instability, the force trace for the higher
Stokes number calculation shows a higher level of small-scale
disturbances. It is not clear whether this is a genuine or purely
numerical effect.
Figures 7.23, 7.25, and 7.27 present plots of drag coefficient for the
isolated edge (C^ - referred to here and in the work of Graham, 1980, as
the coefficient "A") against time for various computer runs. In fact it
is the half-coefficient, A2, which is actually shown plotted at every

2
half time cycle (Equation 7.3.3) and represents the integration of
vortex force coefficient over the previous half-cycle multiplied by two
to yield a coefficient compatible with a full-cycle coefficient. Figure
7.24 shows a plot of the half-coefficient, B; (related to the inertia

2
coefficient, C^), which has been similarly derived.
The method of calculation of "vortex" force takes account not only of
the component due to vortex shedding, but also of that due to skin
friction. However, there are two pieces of evidence to suggest that the
skin friction component is negligible. Firstly, a rough assessment of
the skin friction drag, in the absence of flow separation, was made
using oscillatory boundary layer theory (Section 7.1), and it was found
to be less than 5% of the vortex drag in all cases. Second, comparison
of force and drag coefficient data for runs at different Stokes number
160
reveals no significant trend due to lowering the Stokes number and hence
increasing the relative skin friction contribution, although, since a
lower Stokes number would also affect the pressure drag via changes in
the separation point and amount of diffusion in the wake, it is possible
that an actual significant increase in skin friction contribution has
been partially obscured by a corresponding decrease in pressure drag.
Also, in cases where the flow has developed significantly to cover much
of the body, it becomes difficult to distinguish between drag due to
vortex shedding and that due to skin friction.
SQUARE EDGE
Figure 7.23 shows the half-coefficient A2 results from two sets of

2
square edge computer runs (Runs 1 and 2). Run 2 is at a slightly lower
Stokes number and was carried out to determine whether Run 1, which
could not be extended to longer times because of mesh restrictions, had
been continued long enough to provide an adequate assessment of the
coefficient A for fully developed flow, that is, away from starting
conditions. Assuming the effect of changing Stokes number to be very
slight indeed, it can be seen from Run 2 that drag appears to have
settled down by the third half-cycle. A similar plot for the
half-coefficient B2 is presented in Figure 7.24. This coefficient is a

2
most sensitive quantity and can be seen to vary quite markedly with
time.
Average values of the coefficients A and B, calculated as described in
Section 7.3.3, are 1.398 and 0.017 for Run 2. The same results as
161
computed by the Discrete Vortex Method are 1.995 and 0.051. Most
interest lies in A, which, if we are to take the HMVDM results as
reliable, the DVM has overestimated by some 45%. This overestimation is
considered quite normal. Experiment (see Section 4.2.4) gives A as
approximately 1.4, which would make the HMVDM result seem most
respectable. The result for B is bracketed by the DVM results of on the
one hand the present work and on the other hand Graham (1980) who
computed a value of -0.04.
SQUARE EDGE WITH ROUNDING
Figure 7.25 shows the variation of half-coefficient Ai with edge
rounding over two time cycles. At higher edge non-dimensional radii drag
appears to settle down more quickly. This is especially clear at r/l_z =
0.0924. Average values of A were derived from the above figure using the
method of Section 7.3.3, save that at r/l_z = 0.0165 only one data point
was used, that at t/T = 1.5. At later times in the computation the drag
predictions for this case were considered unreliable due to the
proximity of the vortex pairs to the mesh boundary. Figure 7.26 shows
the variation of these averaged values of coefficient A with
non-dimensional radius. This figure is the HMVDM counterpart of the DVM
plot, Figure 4.22. A comparison of these two figures, again assuming
that the HMVDM is reliable, shows that the DVM overestimates A by an
increasing percentage as r/l_z is increased. This is in accordance with
the problems encountered with the DVM at the higher edge radii and
described in Section 4.3.2.

162
No graph of the variation of coefficient B with non-dimensional edge
radius is presented since, except for the square edge, no discernible
trend may be found in the computed values. In all cases they are very
small and positive.
7.3.7 Effect of variation of Stokes number
The effect of a nearly ten-fold reduction in Stokes number, with all
other computational variables held identical, can be seen in Figure
7.27a by a comparison of the half-coefficient A^ for Runs 4 and 6. The

2
onset of computational instabilities prevented Run 6 from being taken
any further than 1 time cycle, yet on the two data points available for
comparison the effect of this change in Stokes number would appear to be
insignificant. Of course, one should expect larger differences for flows
at very different Stokes number, especially when one enters the
turbulent regime. It is beyond the scope of this work to consider the
effects of turbulence; in laminar flow, a change in Stokes number might
be expected to affect both skin friction and pressure drag, in the
latter case via changes in the separation point and in the amount of
diffusion and mixing in the wake.
7.3.8 Sensitivity to mesh density and time step size
It is necessary to establish the sensitivity of the HMVDM to changes in
mesh density and time step size (At). As was reported earlier,
163
considerable testing of the HMVDM was carried out to verify that it was
insensitive to such parameters. In this work the results of only one
such test are presented. Runs 6 and 7 (Figure 7.27b) differ in that the
mesh for Run 7 (Mesh B) is four times as coarse as that for Run 6 (Mesh
C) and also in that time step for the former is six times as large.
Considering the vast difference in the parameters, and the fact that the
vortex shedding processes presented here are not in any case
deterministic, the HMVDM seems fairly insensitive to mesh and time step
refinement. Also, comparison may legitimately be made with the higher
Stokes number run (4) for the reasons described in the previous section;
such a comparison demonstrates that though especially at time t/T = 1
the coarse mesh predicts a higher drag level than the finer mesh this
discrepancy has virtually disappeared at t/T = 2.
7.4 Application of the Hybrid Moving Vortex Diffusive Method Results to
the Roll Dancing of a Rectangular Cross-Section Barge with Rounded
Bilges
The above described drag coefficient results for the isolated edge with
edge rounding may be applied to the case of the rolling rectangular
barge described in Chapter 5 by the matching technique incorporated into
the computer program BMTIMP also described in that chapter. Graphs of
roll damping coefficient calculated using the Discrete Vortex Method
drag predictions may be found in Figures 5.2, 5.6 and 5.7 for the
standard barge Cases 3S, 3 and 2 respectively. Exactly equivalent graphs
using the HMVDM predictions are given in Figures 7.28, 7.29 and 7.30. In
164
all three cases the HMVDM damping predictions seem to perform well at
low roll velocity amplitude, although there is very considerable scatter
in the experimental data. At higher roll velocity amplitudes, where
comparison is made with experimental forced roll results, the HMVDM
predictions underestimate damping by up to 30% or so, whereas the DVM
predictions significantly overestimate damping throughout the roll
velocity amplitude range.
The possible sources of error in the DVM damping predictions and the
experimental results have been outlined in Section 5.3. Any discrepancy
between theory and experiment here will be due to inaccuracies either in
the HMVDM, the matching process or the experiments themselves. The HMVDM
isolated square edge drag result apparently compares so well with the
square edge experiment (see Section 7.3.6) that the discrepancy between
the theoretical and experimental roll dampings for the rectangular barge
with square bilges (Case 3S) cannot be attributed except in small
measure to inaccuracies in the HMVDM. The underestimation of experiment
by theory in this case at the larger roll velocity amplitudes might
point to a systematic error in the forced-roll experimental results
which could explain the underestimation also present for Cases 3 and 2.
165
CHAPTER 8
CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK
8.1 Conclusions Concerning the Discrete Vortex Method
8.1.1 Isolated edge results
Conclusions concerning the Discrete Vortex Method may be drawn from the
work at two different levels. First the predictions, especially the drag
predictions, of the Discrete Vortex Method can be assessed directly by
comparison with experimental results. Second the predictions of the
BMTIMP computer program, which requires the DVM drag results as input,
can also be assessed by comparison with experiment. This will in turn
provide some insight into the reliability of the Discrete Vortex Method,
but will also test other aspects of the roll damping prediction program
such as the applicability of the matching technique. Conclusions
concerning the BMTIMP computer program may be found in Section 8.3; this
section discusses the conclusions which may be drawn directly relating
to the Discrete Vortex Method results.
Plots of theoretical vortex positions in general follow well the trends
shown in the experimental flow visualisation studies, especially for
edges with keels, where all the essential features are present in both
the shedding and convection phases. In the case of small amounts of
bilge rounding, theoretical flow patterns correspond well to the type of
flow visualised for a square edge, as we should expect. At higher radii,

166
however, the program, in its present form at least, is unable fully to
simulate the double shedding pattern observed experimentally.
Theoretical vortex sizes for flow round a square edge with varying keel
span compare favourably with experimentally measured sizes given the
uncertainty in the experimental results.
Calculated drag coefficients for the square edge and the flat plate are
near Graham's calculated results (Graham, 1980). A well known feature
of the Discrete Vortex Method is that it overestimates drag - see for
example Sarpkaya (1975). This overestimation in the present case is
about 30% for the flat plate and 45% for the square edge. The magnitude
of these figures is in accordance with other workers' findings. The fact
that drag results for the square edge with bilge keel are bounded on
both sides, in the limit of vanishing keel span by the square edge and
in the limit of long keel by the flat plate, by values which have been
compared with experiment and those of other workers must give confidence
that the error bound for the bilge keel results lies within these
limits.
8.1.2 General hull cross-section results
The test computations performed using the DVM as developed for general
hull sections indicated that despite the complexities of the geometry
and the difficulties of implementing the representation of roll the
method was performing well. Drag predictions for a rolling cylinder with
bilge keels were within 16 to 24% of experiment.

167
The results of the computations for the FPV Sulisker were less easy to
validate against experiment because of the fairly high level of
uncertainty in the experimental dampings and the position of the roll
centre. The DVM overestimates experimental total damping by between
roughly 30 and 120%. These two figures represent approximate lower and
upper bounds on the error, with all comparisons falling within these
bounds. Given the fact that the DVM itself would be expected to
overpredict damping, and also that the assumptions of strip theory are
expected to lead to total dampings which are too high, the DVM
predictions appear reasonable. Comparison with the semi-empirical method
of Ikeda et al. (1978) was very poor, with the Ikeda method yielding an
overall damping a factor of four lower than the DVM. The Ikeda method
underpredicts experiment by roughly a factor of three. There appear to
be weaknesses in Ikeda's method when employed for geometries
substantial ly different from those against which it was originally
validated. Comparison with the method of Odabasi et al. (1985) for the
same case is good. One should expect a correctly implemented Discrete
Vortex Method to yield more accurate answers than either of the other
two methods because it would provide a better model of the physical flow
processes. In fact, it is believed that this is the most complete
theoretical model of vortex roll damping for ships yet to be
constructed.
Discrete Vortex Method calculations were performed for the Sulisker with
added bilge keels. Since the damping increment due to the keels was
relatively small, comparisons either with the other damping prediction
methods or with experiment for this case were inconclusive, although the
168
DVM predictions appeared at least feasible.
Damping computations with the DVM for the Sulisker have shown that there
is a substantial damping increment due to the addition of a small bar
keel. The method has also demonstrated the sensitivity of the damping to
changes in the roll centre. It is anticipated that it will be of great
value in determining the effect of even small changes in geometry or
flow parameters.
Computations using the Discrete Vortex Method for other ships would
provide further means of assessing the accuracy of its predictions. One
particularly useful set of tests would involve on the one hand
comparisons between computed and experimental overall vessel dampings,
and on the other hand comparisons between computed sectional dampings
and dampings derived from experiments involving equivalent
two-dimensional sections. This would not only test the accuracy of the
DVM itself, but also the implications of the use of strip theory.
One possible future development of the DVM is its extension to three
dimensions. Although this would require a major research effort, it
would be expected to yield more accurate roll damping data for
three-dimensional vessels than the present strip approach. It would also
enable the effects of forward speed to be investigated.

169
8.1.3 General conclusions
The above work has highlighted many interesting features of the Discrete
Vortex Method as formulated here for both isolated edges and general
ship hull cross-sections.
- The method is very economical even for many cycles of sinusoidal flow
around intricate geometries.
Sinusoidal flow, in that there is no mean convective velocity, and
therefore any mean convection of the vortices is produced through the
shedding and pairing of the vortices themselves, is among the most
difficult types of flow for which the Discrete Vortex Method is used,
and yet feasible and consistent results can still be achieved for a wide
variety of complex geometries especially when the edge is sharp.
- A variety of different quantities can be deduced from the calculations
without too much effort.
- As the edge becomes less "sharp", either by an increase in angle or by
the application of edge rounding, the method becomes less able to cope
for two distinct reasons:
1) The separation point becomes more difficult to define and hence
the Kutta condition becomes more difficult to satisfy in a
meaningful way.
170
2) The subsequent convection of vortices becomes more dominated by
viscous effects. Experiments indicate that for the less "sharp"
edges the vortex pairs do not convect out into the flow but convect
along the edge. Also diffusion takes place rapidly so that the
convection velocities reduce quickly. These effects cannot be
modelled adequately by the Discrete Vortex Method.
Because of the above deficiencies in the method adequate results cannot
be provided by the Discrete Vortex Method employed in oscillatory flow
around significantly rounded edges. In such circumstances a full
Navier-Stokes solution of the flow provides the most realistic
alternative.
8.2 Conclusions Concerning the Hybrid Moving Vortex Diffusive Method
The HMVDM provides a model of the unsteady two-dimensional Navier-Stokes
equations, excluding the effects of turbulence. How good a model of
vortex shedding this then makes depends partly on what Reynolds number
is being considered, partly on what the effects of
three-dimensionalities in the real flows are likely to be and partly on
the adequacy of the finite difference approximations inherent in the
method.
The quality of results from test computations including the oscillatory
boundary layer case was encouragingly good. The HMVDM was then applied
171
to vortex shedding from an isolated square edge with rounding. Plots of
vortex positions and velocity vectors demonstrated the ability of the
method to predict fundamental flow processes correctly. Even small
details of the shedding structure which were noted in experimental flow
visualisations were also apparent in the theoretical predictions. Drag
results displayed virtual insensitivity to mesh refinement, time step
size or Stokes number within reasonable limits. The one direct
comparison with an experimental drag result which was possible - that
for the plain square edge - was excellent, the two results differing by
just 0.14%, although there is a degree of uncertainty in experiment.
Comparison with the DVM drag result for the plain square edge showed
that the DVM result exceeded the HMVDM by about 45%, which would be
expected considering the limitations of the DVM. There was worsening
agreement at higher edge radii, due presumably to inadequacies in the
DVM model. In the limit of infinite edge radius, we retrieve effectively
the problem of the oscillatory boundary layer, which was predicted very
well by the HMVDM, the maximum value of skin friction force being
overestimated by 4% compared with the analytical.
The HMVDM drag results have applicability, via matching processes, not
only to the prediction of the roll damping due to vortex shedding of
barges with rounded bilges, but also to other related problems such as
the prediction of damping from a square cylinder with corner rounding.
The hybrid nature of the HMVDM appears to encompass the best features of
both mesh and particle methods. On the one hand, the plots of velocity
vectors formed on the mesh indicate an entirely satisfactory mesh

172
representation, whilst on the other hand the plots of vortex positions
provide an extremely detailed picture of the vortex shedding processes.
This method would provide an ideal means of analysing in a detailed way
the processes of vortex shedding and making more detailed comparisons
with experiment than was feasible here. Analysis would not need to be
confined to oscillatory flow around isolated edges. It is considered
that the method is capable of extension to other geometries and flow
configurations.
A particular area where more work would be valuable is in extending the
HMVDM to cope with complete ship cross-sections in the same way that the
DVM was extended. This would be a major exercise, but would provide a
useful tool for the assessment of the damping contribution from sections
of ships where no geometrically definable shedding point existed.
One limitation with the present HMVDM calculations was that it was not
possible to continue them for much more than two cycles of oscillation
because of the excessive computational expense involved. Although the
accuracy of the results was not believed to have been seriously impaired
by this, an important area of future work would be to devise ways of
continuing the calculations much longer, in order to assess long term
computational trends.
8.3 Conclusions Concerning the BMTIMP Computer Program
Comparisons with experiment were made for a rectangular cross-section

173
barge with square bilges using both the DVM and the HMVDM drag results
as input to BMTIMP. In general, experimental roll damping values were
bounded by on the one hand the DVM prediction and on the other the
HMVDM. Since we can lay considerable confidence in the HMVDM drag
prediction, this might suggest either that there is a small systematic
error in experiment or that the matching process in BMTIMP has caused
damping to be slightly underpredicted.
The predictions using the HMVDM for barges with rounded bilges
underestimate experiment by up to 30% at the higher roll amplitudes.
Predictions at lower roll amplitudes are often well within the scatter
of experimental points. This again might suggest some doubt about the
accuracy either of BMTIMP or the experimental results themselves.
There are no applicable experimental data available for the rectangular
barge with bilge keels, but the experiments by Ikeda, Komatsu, Himeno
and Tanaka (1977) for barges with bilge rounding and bilge keels compare
extremely well with the BMTIMP theoretical predictions using drag values
predicted by the DVM.
Use of the modified BMTIMP program in forced-roll tests in conjunction
with Discrete Vortex Method and Hybrid Moving Vortex Diffusive Method
isolated edge drag predictions has yielded the following important
findings.
- The modified version of BMTIMP is capable of the prediction, to useful
accuracy, of the roll damping of barges of rectangular cross-section

174
with a variety of bilge geometries.
- Even very small bilge keels produce a significant damping increment
over the square edge or rounded bilge. Particularly noteworthy is the
increase in drag if a "bounded" keel (i.e. of span 0.414 r) is added to
a rounded bilge.
- The trend of the damping versus roll amplitude curve for both the edge
with keel and the rounded bilge is that the damping at low roll
amplitude diverges from that for the square bilge. At the higher
amplitudes, however, the two curves run parallel.
Further experimental evidence would be useful in establishing the
validity of the current modifications in matching to the BMTIMP barge
rolling program and to fix the domain within which the assumptions
underlying the modifications may be considered to be valid. This
experimentation could include forced--roll tests on a rectangular
cross-section barge with different sizes of bilge keel and different
degrees of bilge rounding. The experiments could also include tests on a
freely floating barge in regular beam waves, since the program has also
been modified to cope with this condition.

175
REFERENCES
ABERNATHY, F. H. and KRONAUER, R. E. 1962
The formation of vortex streets.
J. Fluid Mech., 13, 1.
BAKER, G. R. 1979
The "Cloud in Cell" technique applied to the roll up of vortex sheets.
J. Comp. Phys., 31, 76.
BASUKI, J. 1983
Unsteady flows over aerofoils and cascades.
PhD. thesis, Imperial Coll., Aero. Dept..
BAYLISS, A., GUNZBURGER, M. and TURKEL, E. 1982
Boundary conditions for the numerical solution of elliptic equations in
exterior regions.
SIAM J. Appl. Math., 42, no. 2.
BEARMAN, P. W., DOWNIE, M. J., GRAHAM, 0. M. R. and OBASAJU, E. D. 1985
Forces on cylinders in viscous oscillatory flow at low
Keulegan-Carpenter numbers.
BORTHWICK, A. 1986
Comparison between two finite-difference schemes for computing the flow
around a cylinder.
Int. J. for Numerical Methods in Fluids, 6, 275.

176
BROOK, A. K. 1986
The role of simulation in determining the roll response of a vessel in
an irregular sea-way.
RINA International Conference on Ship Stability and Safety: The Safeship
Project, June 1986.
BROOK, A. K. 1987
Pr ivate commun icat io n .
BROWN, C. E. and MICHAEL, W. H. 1955
On slender delta wings with leading edge separation.
NACA Tech. Note 3430.
BROWN, D. T., EATOCK TAYLOR, R. and PATEL, M. H. 1983
Barge motions in random seas - a comparison of theory and experiment.
CHORIN, A. J. 1973
Numerical study of slightly viscous flow.
J. Fluid Mech., 57, pt. 4, 785.
CHORIN, A. J. and BERNARD, P. S. 1973
Discretisation of a vortex sheet, with an example of roll-up.
J. Comp. Phys., 13, 423.
CHRISTIANSEN, J. P. 1973
Numerical simulation of hydrodynamics by the method of point vortices.
J. Comp. Phys., 13, 363.

177
CLEMENTS, R. R. 1973
An inviscid model of two-dimensional vortex shedding.
CLEMENTS, R. R. and MAULL, D. J. 1975
The representation of sheets of vorticity by discrete vortices.
Prog. Aero. Sci., 16, no. 2, 129.
COOLEY, J. W. and TUKEY, J. W. 1965
An algorithm for the machine calculation of complex Fourier series.
Math. Comp., 19, 297.
COZENS, P. D., BEARMAN, P. W. and GRAHAM, J. M. R. 1986
The prediction of roll damping of barges due to vortex shedding: effects
of bilge keels and cross-section shape.
BMT Ltd. Technical Memorandum no. 125.
DALZELL, J. F. 1978
A note on the form of ship roll damping.
J. Ship Research, 22, no. 3, 178.
DAVIS, R.T. 1983
Numerical methods for coordinate generation based on a mapping
technique.
In "Computational methods for turbulent, transonic, and viscous flows",
edited by Essers, J. A., Springer-Verlag, Berlin.

178
DEFFENBAUGH, F. D. and MARSHALL, F. J. 1976
Time development of the flow about an impulsively started cylinder.
AIAA J., 14, no. 7, 908.
DENISE, J-P. F. 1983
On the roll motion of barges.
RINA Supplementary Paper, 125, 255.
DOWNIE, M. J., BEARMAN, P. W. and GRAHAM, J. M. R. 1984
Prediction of the roll damping of barges including the effects of vortex
shedding.
NMI Ltd., Report R188.
DOWNIE, M. J., BEARMAN, P. W. and GRAHAM, 0. M. R. 1985
Prediction of roll damping for a barge in regular waves.
NMI Ltd., Technical Memorandum No. 115.
FALTINSEN, 0. M. and SORTLAND, B. 1987
Slow drift eddy making damping of a ship.
Appl. Ocean Res., 9, no. 1.
FARMER, C. L. 1985a
A moving point method for arbitrary Peclet number multi-dimensional
convection-diffusion equations.
IMA J. Num. Ana., 5, 465.

179
FARMER, C. L. 1985b
Moving point techniques.
Proc. NATO Advanced Study Institute on Fundamentals of Transport
Phenomena in Porous Media, pub. Martinus Nijhoff, Newark, Delaware.
FINK, P. T. and SOH, W. K. 1974
Calculation of vortex sheets in unsteady flow and applications in ship
hydrodynamics.
10th Symp. Naval Hydrodynamics, Cambridge, Mass.
FRANK, W. 1967
Oscillation of cylinders in or below the free surface of deep fluids.
NSRDC, Washington, D. C., Rep. 2375.
FROUDE, W. 1874
On the resistance in rolling of ships.
Naval Science, 3, 312.
GARDER, Jr., A. 0., PEACEMAN, D. W. and POZZI, Jr., A.L. 1964
Numerical calculation of multi-dimensional miscible displacement by the
method of characteristics.
SPEJ, 26.
GAWN, R. W. L. 1940
Rolling experiments with ships and models in still water.
Inst. Nav. Arch., 82.

180
GENTRY, R. A., MARTIN, R. E. and DALY, B. J. 1966
An Eulerian differencing method for unsteady compressible flow problems.
J. Comp. Phys., 1, 87.
GERRARD, J. H. 1967
Numerical computation of the magnitude and frequency of the lift on a
circular cylinder.
Phil. Trans. Roy. Soc., 261, no. 1118, 137.
GRAHAM, J. M. R. 1980
The forces on sharp-edged cylinders in oscillatory flow at low
Keulegan-Carpenter numbers.
J. Fluid Mech., 97, 1, 331.
HAMA, F. R. and BURKE, E. R. 1960
On the rolling up of a vortex sheet.
Univ. of Maryland, Tech. Note BN 220.
HARLOW, F. H. 1964
The particle-in-cell computing method for fluid dynamics.
Methods in Comp. Phys., 3, 319.
IKEDA, Y. and HIMENO, Y. 1981
Calculation of vortex-shedding flow around oscillating cylinder and
Lewis-form cylinders.
Proceedings of the 3rd International Conference on Numerical Ship
Hydrodynamics, Paris.
181
IKEDA, Y., HIMENO, Y. and TANAKA, N. 1977
On eddy making component of roll damping force on naked hull.
J. Society of Naval Architects of Japan, 142.
IKEDA, Y., HIMENO, Y. and TANAKA, N. 1978
A prediction method for ship roll damping.
Dept, of Nav. Arch., Univ. of Osaka, Japan, Rep. 00405.
IKEDA, Y., KOMATSU, K., HIMENO, Y. and TANAKA, N. 1977
On roll damping force of ship - effect of hull surface pressure created
by bilge keels.
J. Kansai Society of Naval Architects, no. 165.
JAROCH, M. 1986
An introduction to the method of discrete vortices and the application
to the problem of modelling the flow past a normal flat plate with long
wake-splitter plate.
Institutsbericht IB 01/86, Hermann-Fottinger-Institut fur Thermo- und
Fluiddynamik, Technische Universitat, Berlin.
KADEN, H. 1931
Aufwicklung einer unstabilen unstetigkeitsflache.
Ing. Arch., 2, 140.
KAMEMOTO, K. and BEARMAN, P. W. 1978
The importance of time step size and initial vortex position in
modelling flows with discrete vortices.
Imperial Coll. Aero. Tech. Note 78-108.

182
KATO, H. 1958
On the frictional resistance to the rolling of ships.
Soc. of Nav. Arch, of Japan J.s 102, 115.
KATO, H. 1966
Effect of bilge keels on the rolling of ships.
Memoirs of the Defense Academy, Japan, 4, no. 3, 369.
KIYA, M. and ARIE, M. 1977
A contribution to an inviscid vortex-shedding model for an inclined flat
plate in uniform flow.
J. Fluid Mech., 82, part 2, 223.
KIYA, M., ARIE, M. and HARIGANE, K. 1979
Unsteady separated flow behind a normal plate calculated by a
discrete-vortex model.
Memoirs of the faculty of Engng., Hokkaido Univ., 15, 199.
KORVIN-KROUKOVSKY, B. V. and JACOBS, W. R. 1957
Pitching and heaving motions of a ship in regular waves.
Trans. SNAME, 65, 590.
KUWAHARA, K. 1973
Numerical study of flow past an inclined flat plate by an inviscid
model.
J. Phys. Soc. Japan, 35, no. 5, 1545.

183
LECOINTE, Y. and PIQUET, J. 1986
Periodic and multiple periodic behaviour of locked in vortex shedding.
16th O.N.R. Symposium on Naval Hydrodynamics, Berkeley, Univ. of
California.
LEONARD, A. 1980
Vortex methods for flow simulation.
J. Comp. Phys., 37, 289.
LEWIS, R. I. and P0RTH0USE, D. T. C. 1983
Recent advances in the theoretical simulation of real fluid flows.
Proc. North East Coast Inst. Engrs. and Shipbuilders, March Ed..
MERCIER, J. A. 1965
Scale effect on roll-damping forces at zero forward speed.
Davidson Laboratory Report 1057.
MILINAZZO, F. and SAFFMAN, P. G. 1977
The calculation of large Reynolds number two-dimensional flow using
discrete vortices with random walk.
J. Comp. Phys., 23, no.4, 380.
MILNE-THOMSON, L. M. 1968
Theoretical hydrodynamics.
Macmillan & Co., 5th ed..

184
MOORE, D. W. 1974
A numerical study of the roll-up of a finite vortex sheet.
J. Fluid Mech., 63, part 2, 225.
NAYLOR, P. J. 1982
A discrete vortex model for bluff bodies in oscillatory flow.
PhD. thesis, Imperial Coll., Aero. Dept..
NOBLE DENTON and ASSOCIATES 1985
Report on the prediction of roll damping of transportation barges.
Final phase - barge motions research programme for Noble Denton and
Associates.
Project no. 20009 (Commercial-in-Confidence).
ODABASI, A. Y., INCECIK, A. and EDWARDS, H. M. 1985
Analytical investigation of roll damping estimation for large amplitude
motions.
BSRA Contract Report W1278 (Commercial-in-Confidence).
OGILVIE, T. F. and TUCK, E. 0. 1969
A rational strip theory of ship motion.
Dept. Nav. Arch., Univ. of Michigan, Rep. 013.
PATEL, M. H. and BROWN, D. T. 1986
On predictions of resonant roll motions for flat-bottomed barges.
RINA Supplementary Papers, 128, 173.

185
PEACE, A. J. and RILEY, N. 1983
A viscous vortex pair in ground effect.
PEACEMAN, D. W. and RACHFORD, H. H., Jr. 1955
The numerical solution of parabolic and elliptic differential equations.
J. Soc. Indust. Appl. Math., 3, no. 1, 28.
PRICE, H. S., CAVENDISH, J. C. and VARGA, R. S. 1968
Numerical methods of higher-order accuracy for diffusion-convection
equations.
SPEJ, 8, 293.
PULLIN, D. I. 1978
The large-scale structure of unsteady self-similar rolled-up vortex
sheets.
J. Fluid Mech., 88, pt. 3, 401.
RAVIART, P. A. 1986
Particle numerical models in fluid dynamics.
In "Numerical Methods for Fluid Dynamics II", edited by Morton, K. W.
and Baines, M. J., OUP, Oxford.
RICHTMYER, R. D. 1957
Difference methods for initial-value problems.
Interscience Publishers, Inc., New York.

186
ROACHE, P. J. 1970
Sufficiency conditions for a commonly used downstream boundary condition
on stream function.
J. Comp. Phys., 6, no. 2, 317.
ROACHE, P. J. 1972
Computational fluid dynamics.
Hermosa, Albuquerque.
ROBINSON, R. W. and STODDART, A. W. 1986
An engineering assessment of the role of non-linearities in
transportation barge roll response.
RINA paper no. 8.
ROSENHEAD, L. 1931
The formation of vortices from a surface of discontinuity.
Proc. Roy. Soc. London, Ser. A, 134, 170.
SAKATA, H., ADACHI, T., SAITO, T. and INAMURO, T. 1983
A numerical analysis of the flow around structures by the Discrete
Vortex Method.
Mitsubishi Heavy Industries, Ltd., Technical Review, October Ed..
SALVADORI, M., and BARON, M. 1961
Numerical methods in Engineering.
Prentice-Hall Inc., Englewood Cliffs, New Jersey.

187
SALVESEN, N., TUCK, E. 0. and FALTINSEN, 0. 1970
Ship motions and sea loads.
Trans. Soc. Nav. Arch. Marine Engrs., 78, 421.
SARPKAYA, T. 1968
An analytical study of separated flow about circular cylinders.
Trans ASME J. Basic Engng, 90, 511.
SARPKAYA, T. 1975
An inviscid model of two-dimensional vortex shedding for transient and
asymptotically steady separated flow over an inclined plate.
J. Fluid Mech., 68, 1, 109.
SARPKAYA, T. and SCHOAFF, R. L. 1979
Inviscid model of two-dimensional vortex shedding by a circular
cylinder.
AIAA J., 17, no. 11, 1193.
SCHLICHTING, H. 1960
Boundary layer theory.
McGraw-Hi11.
SMITH, W. E. and SALVESEN, N. 1970
Comparison of ship-motion theory and experiment for destroyer with large
bulb.
J. Ship Res., 14, no. 1.

188
SOH, W. K. and FINK, P. T. 1971
On potential flow modelling of the action of ships' bilge keels.
4th Aust. Conf. Hyd. Fluid Mech., Monash Univ..
SPOUGE, J. R. 1987
Private communication.
SPOUGE, J. R. and IRELAND, N. 1986
An experimental investigation of large amplitude rolling motion.
RINA International Conference on Ship Stability and Safety: The Safeship
Project, June 1986.
STANSBY, P. K. 1977
An inviscid model of vortex shedding from a circular cylinder in steady
and oscillatory far flows.
Proc. Instn. Civ. Engrs., part 2, 63, 865.
STANSBY, P. K. and DIXON, A. G. 1983
Simulation of flows around cylinders by a Lagrangian vortex scheme.
Applied Ocean Research, 5, no. 3, 167.
THOMAN, D. C. and SZEWCZYK, A. A. 1969
Time-dependent viscous flow over a circular cylinder.
Physics of Fluids Supplement II, 76.
VARGA, R. S. 1962
Matrix iterative analysis.
Prentice-Hall Inc., Englewood Cliffs, New Jersey.

189
WESTWATER, F. L. 1935
The rolling up of the surface of discontinuity behind an aerofoil of
finite span.
R and M 1692, Aero. Res. Council, GB.
WHITE, SIR W. 1895
Notes on further experiences with first class battleships.
Trans. Inst. Nav. Arch., 36, 127.
WOODS, L. C. 1954
A note on the numerical solution of fourth order differential equations.
Aero. Q., 5, pt. 3, 176.
WU, J. C. 1981
Theory for aerodynamic force and moment in viscous flows.
AIAA J., 19, no. 4.

190
Bilge keel spans tested Bilge radii tested

(mm) (mm)
0
1 (chamfered) square edge
2 13
6 25
13
25
76
Figure 2.1 Details of the model barge.

variable
amplitude
Figure 2.2 Schematic representation of experimental apparatus.

192
Figure 2.3 Flow visualization of vortex shedding from a rolling barge

with rounded edges and bilge keels. Bilge radius, r = 13 mm;
keel span, a = 13 mm; roll angular frequency, w = 1.5 rad/sec;
and roll amplitude, fi = 1 1 . 4 degrees.
193
Figure 2.4 Visualization of flow round a rolling barge at successive

stages through one cycle of vortex shedding, using tracer
particles. Barge with square edges; roll angular frequency,
a) = 4.0 rad/sec; and roll amplitude, n = 9.9 degrees.
194

stages through one cycle of vortex shedding, using dye.
Barge with square edges and bilge keels; keel span,
a = 25 mm; roll angular frequency, w = 6.0 rad/sec; and
roll amplitude, n = 5.7 degrees.
CO
cn

stages through one cycle of vortex shedding, using tracer
particles. Barge with square edges and bilge keels;
keel span, a = 76 mm; roll angular frequency, w = 4.0 rad/sec;
and roll amplitude, q = 5.7 degrees.
196

Barge with rounded edges; bilge radius, r = 25 mm; roll
angular frequency, go = 1.5 rad/sec; and roll amplitude,
n = 1 1 . 4 degrees.
197
2 6
•\ Ig,
a B

Barge with rounded edges; bilge radius, r = 13 mm; roll
angular frequency, oo = 1.5 rad/sec; and roll amplitude,
fj = 1 1 . 4 degrees.
198
a 8

Barge with rounded edges and bilge keels; bilge radius,
r = 13 mm; keel span, a = 6 mm; roll angular frequency,
199
1 5
2 6
3 7
B
4

Barge with rounded edges and bilge keels; bilge radius,
r 13 mm; keel span, a = 13 mm; roll angular frequency,
w = 1.5 rad/sec; and roll amplitude, n = 11.4 degrees.
~
200
Figure 2.11 Flow visu al ization using twin dye sources of the vortex
shedding fr om a rolling barge with square edges and
bilge keels. Keel span, a = 76 mm; roll angular frequency,
201
Physical plane Transform plane
Figure 3.1 General numerical Schwarz-Christoffel

transformation.
a
1
(O+Oi)
Transform plane
Figure 3.2 Schwarz-Christoffel transformation applied

to the plain isolated edge.
202
(0 + iVb)
a - (O+Oi)
a ^ (0- iVb)
1^
Transform plane
F i g u r e 3.3 S c h w a r z - C h r i s t o f f e l t r a n s f o r m a t i o n a p pl ie d
to the isolated edge w i t h bi l g e keel.
F i g u r e 3.4 Ma i n fe at ur es of the D i s c r e t e Vort ex Method

as a p p l i e d to bluff bodies.
203
F i g u r e 3.5 Schwarz-Christoffel transformation applied

to a ba si c hull c r o s s - s e c t i o n .
F i g u r e 3.6 R e p r e s e n t a t i o n of roll in the physical plane.

204
/ x (u-iv) at )
F i g u r e 3.7 R e p r e s e n t a t i o n of roll in the

t r an sf or m plane.
F i g u r e 3.8 R e p r e s e n t a t i o n of roll for a ship section

using roll and sw ay components.
PRESSURE COEFFICIENT DISTRIBUTION
PL OT OF V O R T E X POSITIONS
□ LOWER SURF AC E — 11.0

Y
o UPPER S URF AC E
— 10-0 PRESSURE
COEFFICIENT
—9 .0
X
— 8.0 Y
—7 .0
— 6.0 X
z
—5 .0
PO
o
en
(VORTEX LENGTH S C A L E ) * 0-5
F i g u r e 4.1 P r e s s u r e d i s t r i b u t i o n a nd v o r t e x positions. Non-dimensional keel span,

a / L z = 2.2; t ' / T = 0.0,
PRESSURE COEFFICIENT DISTRIBUTION 12.0 PLOT OF VORTEX POSITIONS
LOWER SURF AC E C1-ll .0

UPPER S URF ACE
PRESSURE
COEFFICIENT X
a
□ — 9 .0
o)
□
0 8.0
—
Y
□ 0
—7 .0 X
□
□ o
^•-6 •0
□
□ — 5.0 °ao *o
o j
car Bno _
□
□ D *«o g
— 4.0 aODnn
□
o ro
0 -3 .0 o
cr>
0
— 2.0
o — 1 .0
Q
0.50 l 0.0
o
o cbi -o
o (VORTEX LENGTH SCRLE) * 0-5
0 2.0
o o
O O o
h3.0
■4.0
4 D I S T A N C E FROM S H E D D I N G EDGE
(IN TERMS OF VORTEX LENGTH SCRLE) L5 .0
Figure 4.2 Pressure distribution and vortex positions. Non-dimensional keel span,
a/Lz = 2.2; t'/T = 0.25.
PRESSURE COEFFICIENT DISTRIBUTION r-12.0 PLOT OF VORTEX POSITIONS
LOWER S U R F AC E — 11.0
UPPE R SURF AC E
— 10.0 P R E SS U R E Y
COEFFICIENT
-9 .0
b-8.0 X Y
X
-7 .0
<£-6.0
a
[P -5 .0
-4 .0
o
o
O □ — 3.0
207
o
o
o □ 2-0
o -
o o o □
C3
o -1.0
O O
□
0.50 l ■0.0
□
□
-1 -0
□ [3
(VORTEX LENGTH S C A L E ) * Q.5
□ -2.0
□
□
□ □ □ □ -3.0
C14 . 0
D I S T A N C E FROM S H E D D I N G EDGE
(IN TERMS OF VORTEX LENGTH SCRLE) -5 .0
a/L, = 2.2; t'/T = 0.50 .
PRESSURE COEFFICIENT DISTRIBUTION PLOT OF VORTEX POSITIONS
r-12-0
□ LOWER S U R F R C E — 11.0
Y
o UPPER SURF RC E
— 10.0 PRESSURE
COEFFICIENT
— 9.0 X
Y
(£-8.0
-7 .0
n□ □o □
_
© o o13 n B®
om□»
oo — 6.0
o
o
o — 5.0
o
o
o CL-4.0
O
O
0 —3.0
□
208
— 2-0
□
— 1.0
□
0.50 l EH3.0
□ -1 .0
□ -2.0
(VORTEX LENGTH SCRLE) * 0-5
□
[13 .0
□
□ -4.0
□ □ □
□ □ □ □ □
05.0
------------ 0 1STf i NCE FROM S H E D D I N G E D GE
(IN TERMS OF VORTEX LENGTH SCR L E ) -6.0
a/L = 2.2; t'/T = 0.75 ,
209
Figure 4.5c t'/T = 0.5
Figures 4.5 Velocity vectors at successive stages through one time cycle. Non-dimensional
edge radius, r/l_z = 1.0; non-dimensional keel span, a/Lz = 1.0.
210
Beginning of time cycle
Half way through time cycle

OOOo
o O
o
Figure 4.6 Comparison between theoretical vortex positions and experimental flow visualization;
rounded edges with bilge keels, r/[_z = 0.039, a/l_z = 0.039.
211
Figure 4.7 Comparison between experimental and theoretical vortex

sizes plotted as a function of non-dimensional keel
span, a/Lz .
O experiment; — *— theory
212
Figure 4.8a Trace of time dependent vortex force coefficient

(Cp ) against time in cycles (t/T).
Non-dimensional keel span, a/l_z = 0.2.
Figure 4.8b Trace of time dependent vortex force coefficient

(Cpv ) against time in cycles (t/T). Square edge.
SQ U ARE ED G E ROUNDED EDGE WITH
W ITH KEEL BOUNDED KEEL
213
■©— square edge with keel; Q— rounded edge with bounded keel (a/r = /2 - 1).
SQ U A RE EDGE ROUNDED ED G E WITH
WITH K E E L BO UNDED K E E L
214
■©— square edge with keel; ■g— rounded edge with bounded keel (a/r = /2 1 ).
215
Figure 4.12 t'/T = 0.25
Figure 4.14 t'/T = 0.75
Figures 4.11 - 4.14 show computed vortex positions for a rounded edge at
successive stages through a time cycle, Non-dimensional edge
radius, r/Lz = 0.039.
216
Figure 4.16 t'/T = 0.5
Figures 4.15 - 4.17 show computed vortex positions for a rounded edge at
successive stages through a time cycle. Non-dimensional edge
radius, r/L = 0.109.
PRESSURE C OEFFICIENT- D I S T R I B U T I O N — 9 .0 PLOT OF VORTEX POSITIONS
LOWER S URF A C E
— 8.0
U PPE R SURF AC E
P R E S SU RE
— 7. 0 COEFFICIENT
X
- 6.0
— 5. 0
Y
-4.0
— 3- 0
O ©
O
O O
O 2.0
Q -
217
O O — 1.0
O
Y *
□ Q
Figure 4.18 Pressure distribution and vortex positions. Non-dimensional edge radius,
r/Lz = 0.086; t'/T = 0.0.
PRESSURE C O E F F I C I E N T D I S T R I B U T I O N — 9.0 PLOT OF VORTEX POSITIONS
LOWER SURF AC E
U PPE R S URF A C E
P R E SS U R E
COEFFICIENT
— 6.0
— 5.0
— 4.0
ro
oo
(VORTEX LENGTH SCALE; * 0-5
Figure 4.19 P r e s s u r e d i s t r i b u t i o n and v o r t e x positions. Non-dimensional e d g e radius,

r/l_z = 0 . 0 8 6 ; t ’/T = 0.25.
PRESSURE COEFFICIENT DISTRIBUTION — 9.0 PLOT OF VORTEX POSITIONS
LOWER S U R F AC E
— 8.0
UPPE R S U R F AC E
PRESSURE
_7 . 0 COEFFICIENT
— 6.0
— 5.0
— 4.0
— 3.0
— 2.0
219
— 1 .0
Figure 4.20 Pressure distribution and vortex positions. Non-dimensional edge radius,
r/Lz = 0.086; t'/T = 0.5.
PRESSURE COEFFICIENT DISTRIBUTION — 9.0
PLOT OF VORTEX POSITIONS
□ LOWER S U R F A C E
— 8.0
o UPPER S U R F AC E
P R E SS U R E
COEFFICIENT
O O
220
-2.0
9-3.0
[&4 .0
(VORTEX LENGTH SCREE) * 0-5

------ DI STANCE FROM SHEDDING EDGE
(IN TERMS OF VORTEX LENGTH SCRLE! •5.0
F i g u r e 4.21 Pressure distribution and vortex positions. Non-dimensional edge radius

r / L z = 0 . 0 8 6 ; t ’/T = 0.75.
221
I_________ i_________ i_________ i_________ i_________ i

0 0*02 0-04 0-06 0-08 0*10
r / Lz
Figure 4.22 Plot of coefficient A against non-dimensional

edge radius, r/L , for rounded edges, calculated
using the Discrete Vortex Method.
222
f)1 = su rge, ij2 = s w a y , t \3 = h e a v e
t \4 = roll , i)5 = pitch, \ = y a w
Figure 5.1 Schematic representation of barge motions

(reproduced from Downie, Bearman and
Graham, 1984).
FULL SCALE SECTION
CASE 3S DIMENSIONS
G ------- THEORY SHARP CORNER _1

JJ-Tm
+ MODEL TEST TRANSIENT RESULTS
27-Am
o MODEL TEST FORCED ROLLING RESULTS
223
Figure 5.2 Variation of the roll damping coefficient with roll velocity amplitude
cor| for barge Case 3S (model scale values). Theory calculated using the
Discrete Vortex Method.
B /A
V/fU x 100°/o
Bc
220
200
180
160
140
224
120
100
80
60
n a / qZ>/3
40
0 001 0*01 0*1 1-0 10
Figure 5.3 Variation of vortex damping coefficient, per unit of roll, as a percentage
of critical damping coefficient, with roll amplitude/(keel span)^^.
Forced roll test.
( Nm / Rad / S e c )
VORTEX DAMPING COEFFICIENT B„
225
Figure 5.4 Variation of By with for keels of different spans (model scale).
Full scale keel spans: . ,— Q — . 2m; — e — 1m; 0.5m; — 0.2m;

— ^— 0.1m; — X— no keel - equivalent to Case 3S.
ro
ro
cn
Figure 5.5 Percentage gain in-vortex damping over square edge case for'keels of different
spans, as a function of roll velocity amplitude (model scale).
Full scale keel spans: ■e— 2m; ■B— lm; ■A— 0.5m; — §— 0.2m,
— 0.1m.
CASE 3 F U L L S C A L E SECTION
-----------------ROBINSONS EMPIRICAL CORRECTION D IM E N S IO N S
-------- 0 ------ THEORY SHARP CORNER BILGE RADIUS 0-48m
— -o ------ THEORY ROUNDED BILGE
2 -7m
o
MODEL TEST FORCED ROLLING RESULTS
274m !
no
ro
OJrj^ (Rad/Sec)
ojtV for barge Case 3 (model scale values). Theory calculated using the
k
CASE 2 FU LL S C A L E SECTION
----------------- ROBINSONS EMPIRICAL CORRECTION D IM E N S IO N S
-------- 0 ------ THEORY SHARP CORNER BILGE RADIUS 0A8m
-------- 0------- THEORY ROUNDED BILGE
X M 2-9m
27-Am l
228
can for barge Case 2 (model scale v a l u e s ) -. Theory calculated using the
4
( Nm / Rad /Sec)
229
GO (Rad /Sec )
Fi g u r e 5.8 V a r i a t i o n of B v wi th roll an gu la r f r e q u e n c y w for Ikeda's barge model B;
(beam 0.28m, d r a f t 0.1T2m, keel span a = 0.5 cm, bilge radius r = 1cm).
U
<b ----- THEORY
U1 0-227 Rad
O MODEL TEST
(Nm/Rad/
ro
co
o
Fi g u r e 5.9 V a r i a t i o n o f B v wi th roll a n g u l a r f r eq ue nc y oj for Ikeda's barge model C;
(beam 0.28m, d r a f t 0.112m, keel span a = 1cm, bilge radius r = 2 cm).

VORTEX DAMPING COEFFICIENT Bw ( N m / R a d / S e c )
F i g u r e 5.10 V a r i a t i o n of Bv wi t h roll a m p l i t u d e h4 for Ikeda's barge model C;
(beam 0.28m, d r a f t 0.112m, keel span a = 1 cm, bilge radius r = 2 c m ) .

232
F i g u r e 5.11 D i a g r a m of model cylinder.

(
233
TAWS
i i
0
rri'T
5
1 1 t
V
□
1S
SCAIE H HETPfS
Figure 5.12 F.P.V. Sulisker general arrangement.

(
234
F i g u r e 5.13 F.P.V. Sulisker cross-sectional data.
235
Figure 5.14 Computed vortex positions for a selection of sections from the F.P.V. Sulisker.
Run 1, t/T = 4.5 .
S e c t i o n 15 Sect i o n 17
S e c t i o n 13
236
Section 21 Section 22
Section 19
*,*
F i g u r e 5.14 (cont'd)
237
Figure 5.15 V e l o c i t y v e c t o r s at s u c c e s s i v e s t a g e s t h r o u g h o ne t i m e c y c l e .
R un 1; S e c t i o n 17 of t h e F.P.V. S u l i s k e r .
2 38
Figure 5.16 V e l o c i t y v e c t o r s at s u c c e s s i v e s t a g e s t h r o u g h o n e t i m e c y cle.

Run 3; S e c t i o n 17 o f t h e F.P. V . S u l i s k e r .
239
Figure 5.17 V e l o c i t y v e c t o r s at s u c c e s s i v e s t a g e s t h r o u g h o ne t i m e cycle.

R u n 5; S e c t i o n 17 o f t h e F . P.V. S u l i s k e r .
2 40
F i g u r e 5.18 V e l o c i t y vectors at s u c c e s s i v e stages thro u g h one time cycle.

Run 7; S ection 14 of the F.P.V. Sulisker, with bilge keels.
NON-D IM ENSIONAL SECTION O R D IN A T E
(tonnes m l rad )
Figure 5.19 Distributions of sectional vortex rol 1 damping moment coefficient, .c d

for Discrete Vortex Method Runs 1-3 and 5 and for the BMT method. F.P.V.
Sulisker.
242
Ax m-1. Ax m
X_ 9
m-Z
Xm 1
m-1
Xm
m
Xm
m+1
- xm 0
m+2
Figure 6.1 Definition of mesh spacing employed

in the derivation of finite difference
formulae.
243
Figure 6.2 Flow chart of the Hybrid Moving Vortex

Diffusive Method.
244
Physical plane Transform plane
Figure 6.3 Application of the Schwarz-Christoffel transformation to the generation of

computational meshes. In the transform plane the x-mesh spacing is regular;
the y-mesh spacing expands exponentially outwards.
245
1.2 r..2
A y..i
t.i r..,
Ay
K
\ ^ -\* \ \— ^ ^ — 'n— ^ .
Typical velocity
profile near the
wal 1
Figure 6.4 The no-slip vorticity boundary condition.
O u te r (upper) boundary
cons ta n t j = MY
t
In le t Outlet
(upstream) (downstream)
Stream
boundary boundary
Ay, direction
y T T
X 1 =0
A \ \ \ \ \ \ \ ^~ V \ \ "V V V \ \ \ \ j=1
*“ Wall boundary i=NX
Figure 6.5 Boundary conditions on stream function.

246
\ = V W \
Figure 6.6 The velocity interpolation scheme.
A I = A 1 + A 2 + A 3+ A 4
Figure 6.7 Redistribution of point vorticity onto

the mesh system.
---------------------i--------------------- --------------------- 1--------------------- -
(Kj+1) | (i,j+1) j (i+1 j+1)
I 1
I I l 1
I I ___________ I 1
i i
i
x tu ) !
( K j) i /£ . H I
k \/ 1 k \ \ \ v
! ^ y 7 i
l i
i 1
i i I 1
i i I 1
1
(ij-1 ) j (i+1,j-1) .(i - 1 , j - 1 ) |
____________ 1____________
KEY : Inclusion zone
Figure 6.8a First Figure 6.8b Second

amalgamation of amalgamation of
like-signed vortices like-signed vortices
Figure 6.8 Definition of the inclusion zone for the various stages of the
vortex amalgamation process.
2 48
Figure 7.1a Oscillatory boundary layer over a

flat plate: computational mesh.
Figure 7.1b Oscillatory boundary layer over a

flat plate: velocity vectors over
the computational domain at the end
of the second time cycle.
249
□ Computational re su lt
----- Analytical re su lt
Figure 7.1c Oscillatory boundary layer over a flat plate:

detail of velocity distribution at one x-station
at the end of the second time cycle.
250
KEY : ---------Computational result; --------- Analytical result
Figure 7.2 Oscillatory boundary layer over a flat plate: comparison between analytical and
computed wall shear stress".
(
Figure 7.3 Plain square edge in uniform flow: computational mesh.

(
Key: ..... vortex sheet; — *— streamlines
Figure 7.4 Plain square edge in uniform flow: Comparison of flow character!’sties as
calculated by the Discrete Vortex Method (reproduced from Pull in, 1978)
and the Hybrid Moving Vortex Diffusive Method.
< (
253
Figure 7.5 Computational Mesh A.
254
Figure 7.6 Computational Mesh B .
255
Figure 7.7 Computational Mesh C ,
1
256
Figure 7.8 Computational Mesh D .
257
Figures 7.9 Vortex positions at successive stages through two time cycles.
Run 1; Square edge.
258
Figures 7.9 (cont'd)

259
Figures 7.10 Velocity vectors at successive stages through two time cycles.
Run 1; Square edge.
260

261
Figures 7.11 Vortex positions at successive stages through two and a half time
cycles. Run 2; Square edge.
262

263
Figures 7.11 ( c o n t ’d)

264
Figures 7.12 Velocity vectors at successive stages through two and a half time
cycles. Run 2 ; Square edge.
265

266
F i g u r e s 7.12 (cont'd)
267
Figures 7.13 Vortex positions at successive stages through one and three quarter
LO
time cycles. Run 3; non-dimensional edge radius, r/Lz = 0.016

268
Figures 7.13 ( c on t' d)

269
Figures 7.14 Velocity vectors at successive stages through one and three quarter
time cycles. Run 3; non-dimensional edge radius, r/L = 0.0165.
270

271
Run 4; non-dimensional edge radius, r/l_z = 0.0466.
212

273
F i g u r e s 7.16 V e l o c i t y ve ct or s at s u c c e s s i v e stages t h ro ug h two t i me cycles.

Run 4; n o n - d i m e n s i o n a l edge radius, r/L = 0.0466.
274

275
276

277
Figures 7.18 Velocity vectors at successive stages through two time cycles.
Run 5; non-dimensional edge radius, r/Lz = 0.0924.
278

279
Figures 7.19 Vortex positions at successive stages through one time cycle.
280
Figures 7.20 Vortex positions at successive stages through one time cycle.
Run 7; non-dimensional edge radius, r/L = 0.0466.
281
O
CM
CM
Figure 7.21 Traces of time depe nd en t vortex force coef fi ci en t

(Cpv ) ag a i n s t time in cycles (t/T) for d i ff er en t
values of non-dimensional edge radius, r/L .
282
F i g u r e 7. 2 2 Tr a c e s of time d e p e n d e n t vort ex force c o e f f i c i e n t

(Cpv ) a g a i n s t time in cycles (t/T) for d i f f e r e n t
S t o k e s numbers, m e s h e s and time step sizes.
283
Fi g u r e 7.23 H a l f - c o e f f i c i e n t , A^, versus time in cycl es (t/T)
C o m p a r i s o n b e t w e e n Runs 1 and 2.
0.0300
STOKES NUMBER
□ 28300, RUN 1
o 18400, RUN 2
T-frVl
m 0.0225 a 0
o
L l_
L l_ 0.0150
LJ
O
O
I 0
Li_
_I
< 0.0075
X
Q
□
0.0000
0 .0 0.5 1 .0 1 .5 2 .0 2.5
TIME (IN CYCLES)

F i g u r e 7. 2 4 Half-coefficient, B 1# versus time in cycl es (t/T)
2
C o m p a r i s o n b e t w e e n Runs 1 and 2.
284
2.0 J.
NON-DIMENSIONAL RADIUS
□ 0.0. RUN 2
o 0.0165. RUN 3
A 0.0466. RUN 4
+ 0.0924. RUN 3
HALF-COEFFICIENT A.
1 .5 -
l .0 -
0.5 -
+ ---------- + ---------- +
0.0 - — i—
0. 0.5 1 .0 I .5 2 .0 2.5
TIME (IN CYCLES)
Figure .25 Half-coefficient, , versus time in cycles (t/T).

Comparison between Runs 2 to 5.
285
Figure 7.26 Plot of coefficient A against non-dimensional edge

radius, r/l_z , for rounded edges, calculated using
the Hybrid Moving Vortex Diffusive Method.
286
0-0
STOKES NUMBER
a 36700. RUN 4
< 0.6 © 4600. RUN e
LlJ
O
Ll
Ll 0.4 -
LlJ
0
o
1
Ll
_I
< 0.2 -
0.0
0.0 0.5 1 .0 1 .5 2 .0 2-5
TIME (IN CYCLES)
Figure 7.27a Half-coefficient, A 1S versus time in cycles (t/T)

2
Comparison between Runs 4 and 6.
0.8
MESH TYPE AND TIME STEP SIZE
□ C. d T -0 .0 0 0 3 3 , RUN 4
© C. A T -0 .0 0 8 3 3 , RUN 6
A B, 4 T -0 .0 5 . RUN 7
< 0.6 -
f-
z
UJ
0
Ll
Ll
Ll I
O
O
1
Ll
_l
<
X
0 .0 - l - T —r - .—T ■ .................................................I ■■ . . - I - I ............................ .......... ...................

0 .0 0.5 1 -0 1.5 2.0 2 .5
TIME (IN CYCLES)

Figure 7.27b Half-coefficient, A,, versus time in cycles (t/T)
2
Comparison between Runs 4, 6 and 7.
FULL S C A L E SECTION
C A SE 3 S D IM E N S IO N S
0 ------- THEORY SHARP CORNER

I 2*7m
27-Am
287
con for barge Case 3S (model scale values). Theory calculated using the
Hybrid Moving Vortex Diffusive Method.
CASE 3 F U L L S C A L E SECTIO N
-----------------ROBINSONS EMPIRICAL CORRECTION D IM EN S IO N S
-------- 0 ------ THEORY SHARP CORNER BILGE RADIUS 0-48m
-------O-------- THEORY ROUNDED BILGE
2-7m
27-4m
288
wri for barge Case 3 (model scale values). Theory calculated using the
CASE 2 FU LL S C A L E SECTION
-----------------ROBINSONS EMPIRICAL CORRECTION D IM E N S IO N S
-------- 0 ------ THEORY SHARP CORNER BILGE RADIUS 0 A8m
-------- O— THEORY ROUNDED BILGE J2-9m
27-4m ~7
con for barge Case 2 (model scale values). Theory calculated using the

Cozens PD 1987 PHD Thesis

Uploaded by

Copyright:

Available Formats

Cozens PD 1987 PHD Thesis

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Cozens PD 1987 PHD Thesis

Uploaded by

Copyright:

Available Formats

NUMERICAL MODELLING OF THE

ROLL DAMPING OF SHIPS DUE TO

Paul D e n n i s C o z e n s , M.A. (Cantab.)

Flow visualization of the vortex shedding from the bilges of a rolling

rectangular cross-section barge with bilge keels

The second model is an extension of the basic Discrete Vortex Method to

Finally, a more advanced computer model of the full two-dimensional unsteady

Flow visualisation experiments were undertaken to study the vortex shedding

It is with a deep sense of gratitude that I write these

acknowledgements. The very considerable rewards of the research work

itself have been added to immeasurably by the kindness of all those at

Imperial College who have made it possible. In the first instance, I

count myself most fortunate that I should have been supervised by

Professor P. W. Bearman and Dr. J. M. R. Graham who have both been an

unquenchable source of ideas, inspiration and guidance.

I am most grateful to BMT Ltd. and the Marine Technology Directorate of

for their advice and support.

I am indebted to Mr. J. O'Leary for his patient assistance in setting up

the experiments, Mrs. R. Fairhurst for typing many parts of my thesis,

and Mr. J. D. Powell for proof reading.

My thanks also to Dr. M. J. Downie, Dr. J. P. Felix and Mr. P. Kalkanis

Last, but by no means least, my heartfelt thanks to my wife, Elizabeth,

1.1 The Prediction of the Motions of a Vessel in a Seaway 15

1.1.1 The prediction of vessel response in waves using

1.1.2 The roll response of a vessel at zero forward speed 16

1.2 The Discrete Vortex Method 19

1.2.1 Early development of the Discrete Vortex Method 19

1.2.2 The application of the Cloud-in-Cell technique to the

Discrete Vortex Method 24

1.3 The Solution of the Full Two-Dimensional Unsteady

1.3.1 Eulerian solutions of the Navier-Stokes equations 27

1.3.2 Lagrangian solutions of the Navier-Stokes equations 29

1.3.3 Comparison of Eulerian and Lagrangian schemes 31

1.3.4 Hybrid Lagrangian/Eulerian solutions of the

1.4 The Application of Vortex Methods to Vortex Shedding

from Vessel Hulls 33

2. DESCRIPTION OF THE PROCESSES OF VORTEX SHEDDING

2.1 Experimental Method 38

2.1.1 Experimental apparatus 38

2.1.2 Test procedure 38

2.1.3 Flow visualization techniques 40

2.2 Photographs of the Vortex Shedding Processes 41

3. THE DISCRETE VORTEX METHOD 46

3.1.1 Plain isolated edge transformation 47

3.1.2 Plain isolated edge with bilge keel transformation 48

3.1.3 Full numerical solution of the Schwarz-Christoffel

3.2 Basic Concepts of the Discrete Vortex Method 54

3.3 Treatment of the Nascent Vortex 56

3.3.1 Satisfying the Kutta Condition 56

3.3.2 Determination of the position of the nascent vortex 57

3.4 A Scheme for Amalgamating Vortices 59

3.5 Starting the Flow Pattern 61

3.6 Decaying the Vortex Pairs 61

3.7 Calculation of Pressure Coefficient 63

3.8 Calculation of Vortex Force 64

3.9 Extension of the Discrete Vortex Method to Model the

Vortex Shedding from Rolling Ships 64

3.9.1 Vessel in pure sway 64

3.9.2 The representation of roll 65