Computational Methods For Hydrofoil Design - A Composite Analysis Using Panel Code and RANS

Download as pdf or txt
Download as pdf or txt
You are on page 1of 72

See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/271847208

Computational Methods for Hydrofoil Design - A composite analysis using


panel code and RANS

Thesis · January 2015

CITATION READS

1 20,596

1 author:

Gilberto Besana
University of Southampton
1 PUBLICATION 1 CITATION

SEE PROFILE

All content following this page was uploaded by Gilberto Besana on 06 February 2015.

The user has requested enhancement of the downloaded file.


University of Southampton
Faculty of Engineering and the Environment

3rd Year Individual Project

Computational Methods for Hydrofoil Design


A composite analysis using panel code and RANS

Gilberto Besana MEng Ship Science


Professor Stephen R. Turnock Supervisor

January 2015
Abstract
A composite investigation with panel codes and Reynolds-Averaged Navier-Stokes (RANS)
method is conducted on 3D and 2D hydrofoil geometries and NACA sections to quantify the
forces and the performance of these structures. This methodology is aimed to prove that it is
possible to combine different numerical methods for a time-effective hydrofoil design approach.
This combination will use panel methods and thinship theory to investigate 3D NACRA F20
C-foil geometries. RANS, with k - ω SST as turbulent model, is used to analyse 2D NACA 0012
sections. The results will be compared between the two methods and with the experimental
data available. Results proved to be satisfactory being able to encompass in the study many
different performance criteria of hydrofoil design. The data between the codes showed similar
trends although some differences were encountered. A design methodology has been presented
to show how each code is adapt to study which part of the physics involved with this flow
problem. Time proved to be a major restriction but the main objectives of the project have
been met. A consideration on future works was proposed based on the outcomes of the research.

1
Aknowledgements
I would like to thank Professor Stephen R. Turnock for his constant support and valuable advice
throughout the project. His shared knowledge with myself helped carry out the work exposed
in this paper.
Furthermore I’d like to mention two other people who have supported my work. The first is
Dr. Dominic Taunton, who wrote and made available to me Wavelstl, the Matlab code utilised
throughout my project. Secondly Artur Lidtke, a postgraduate student, who gave advice and
guidance for OpenFOAM as well as sharing with me a python script which generates 2D airfoil
meshes for OpenFOAM.

2
Contents
Abstract 1

Aknowledgments 2

List of Figures 5

List of Tables 7

Nomenclature 8

1 Introduction 9
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2 Aims and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Hydrofoils and Foil Sailing Background . . . . . . . . . . . . . . . . . . . . . . . 11
1.4 Development of Computational Fluid Dynamics . . . . . . . . . . . . . . . . . . . 11
1.5 Risk Assessment and Ethical Approval . . . . . . . . . . . . . . . . . . . . . . . . 12
1.6 Layout of Project . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2 Hydrofoils 13
2.1 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.1 Lift Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1.2 Drag Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1.3 Parameters affecting Lift and Drag . . . . . . . . . . . . . . . . . . . . . . 15
2.1.4 Hydrofoil Design and Performance . . . . . . . . . . . . . . . . . . . . . . 16

3 Governing Principles of CFD 18

4 Potential Flow Theory Analysis 20


4.1 Palisupan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.1.1 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.1.2 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1.3 NACRA F20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.1.4 NACA 0012 Hydrofoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 Wavelstl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2.1 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2.2 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.2.3 NACRA F20 foils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.2.4 NACA 0012 horizontal foil . . . . . . . . . . . . . . . . . . . . . . . . . . 25

5 RANS Analysis 26
5.1 OpenFOAM Pre-processing & Guidelines . . . . . . . . . . . . . . . . . . . . . . 26
5.2 Flow around foil section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5.2.1 Mesh grid independence analysis . . . . . . . . . . . . . . . . . . . . . . . 28
5.2.2 Varying angle of attack analysis . . . . . . . . . . . . . . . . . . . . . . . . 30
5.3 Free surface flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.3.1 Varying submergence depth . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3
6 Results 33
6.1 NACA 0012 horizontal foil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
6.1.1 Investigation of Panel size and distribution . . . . . . . . . . . . . . . . . 33
6.1.2 Angle of attack variation . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
6.1.3 Free surface study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.2 NACRA F20 hydrofoils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
6.2.1 Investigation of Panel size and distribution . . . . . . . . . . . . . . . . . 40
6.2.2 Angle variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
6.2.3 Free surface deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

7 Conclusion 44

8 Recommendation for future work 44

References 45

Appendices 47
A Hydrofoils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
B Plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
C Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
D Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

4
List of Figures
1.1 Drag reduction employing hydrofoils . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.2 Early and recent developments of vessels adopting hydrofoil configuration . . . . 11
2.1 Schematic of physics of a submerged hydrofoil under the free surface . . . . . . . 13
2.2 Lift and drag force vectors for a hydrofoil catamaran . . . . . . . . . . . . . . . . 14
2.3 Breakdown of drag components . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.1 Sectional point order definition for Palisupan . . . . . . . . . . . . . . . . . . . . 21
5.1 Lift coefficient variation at increasing mesh refinement . . . . . . . . . . . . . . . 29
5.2 Drag coefficient variation at increasing mesh refinement . . . . . . . . . . . . . . 29
5.3 Lift/Drag ratio curve over angle of attack showing the maximum value at the
given angle obtained. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
6.1 Lift [top] and Drag [bottom] coefficient panel number convergence for Palisupan . 33
6.2 Lift coefficient at varying AoA for a NACA0012 foil . . . . . . . . . . . . . . . . 34
6.3 Drag coefficient at varying AoA for a NACA0012 foil . . . . . . . . . . . . . . . . 34
6.4 CL 2 versus CD plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6.5 CP distribution comparison at 0◦ AoA . . . . . . . . . . . . . . . . . . . . . . . . 35
6.6 CP distribution comparison at 10◦ AoA . . . . . . . . . . . . . . . . . . . . . . . 36
6.7 CP distribution comparison at 15◦ AoA . . . . . . . . . . . . . . . . . . . . . . . 36
6.8 Influence of foil depth from free surface on Drag coefficient studied with OpenFOAM 37
6.9 Influence of foil depth from free surface on Lift coefficient studied with OpenFOAM 37
6.10 OpenFOAM CP distribution for a range of depth . . . . . . . . . . . . . . . . . . 38
6.11 Free surface plots at varying immersion depths . . . . . . . . . . . . . . . . . . . 38
6.12 Wavelstl free surface deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6.13 Free surface plots at immersion depths = c . . . . . . . . . . . . . . . . . . . . . 39
6.14 Drag coefficient panel number convergence for Palisupan . . . . . . . . . . . . . . 40
6.15 Lift coefficient panel number convergence for Palisupan . . . . . . . . . . . . . . 40
6.16 Drag force of NACRA F20 foils at a range of angles for yaw and pitch analysed
with Palisupan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6.17 Lift force of NACRA F20 foils at a range of angles for yaw and pitch analysed
with Palisupan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6.18 Lift force of NACRA F20 foils at a range of angles of heel analysed with Palisupan 41
6.19 Wave resistance variation respect to depth . . . . . . . . . . . . . . . . . . . . . . 42
6.20 Diagram of the proposed methodology for hydrofoil design . . . . . . . . . . . . . 43
A.1 F20 NACRA C-foil geometry and coordinate system . . . . . . . . . . . . . . . . 48
A.2 Foil configuration and NACRA F20 sailing . . . . . . . . . . . . . . . . . . . . . . 49
A.3 Fully submerged and surface piercing hydrofoil example . . . . . . . . . . . . . . 50
B.4 Detail of Figure 6.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
B.5 Drag force of NACRA F20 foils at a range of angles of yaw analysed with Palisupan 51
B.6 Lift force of NACRA F20 foils at a range of angles of pitch analysed with Palisupan 52
B.7 Wavelstl total resistance experienced by the foil at a range of angles of yaw
investigated with Wavelstl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
C.8 Mesh and coordinate system for first OpenFOAM analysis . . . . . . . . . . . . . 53
C.9 Free surface OpenFOAM mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
C.10 Detail of mesh around hydrofoil - OpenFOAM . . . . . . . . . . . . . . . . . . . 55
C.11 Division domain of the air-water interface with use of the setFIelds utility . . . . 55
C.12 Detail of region with geometry issues . . . . . . . . . . . . . . . . . . . . . . . . . 56
C.13 Improved geometry at the foil nose . . . . . . . . . . . . . . . . . . . . . . . . . . 56
C.14 Velocity magnitude field plot for Re = 6 million at 0◦ angle of attack. . . . . . . 57

5
C.15 Pressure field plot for Re = 6 million at 0◦ angle of attack. . . . . . . . . . . . . 57
C.16 Velocity magnitude field plot for Re = 6 million at 1◦ angle of attack. . . . . . . 58
C.17 Pressure field plot for Re = 6 million at 1◦ angle of attack. . . . . . . . . . . . . 58
C.18 Velocity magnitude field plot for Re = 6 million at 2◦ angle of attack. . . . . . . 59
C.19 Pressure field plot for Re = 6 million at 2◦ angle of attack. . . . . . . . . . . . . 59
C.20 Velocity magnitude field plot for Re = 6 million at 3◦ angle of attack. . . . . . . 60
C.21 Pressure field plot for Re = 6 million at 3◦ angle of attack. . . . . . . . . . . . . 60
C.22 Velocity magnitude field plot for Re = 6 million at 4◦ angle of attack. . . . . . . 61
C.23 Pressure field plot for Re = 6 million at 4◦ angle of attack. . . . . . . . . . . . . 61
C.24 Velocity magnitude field plot for Re = 6 million at 5◦ angle of attack. . . . . . . 62
C.25 Pressure field plot for Re = 6 million at 5◦ angle of attack. . . . . . . . . . . . . 62
C.26 Palisupan pressure coefficient distribution for the F20 c-foils obtained through
Paraview at an angle of yaw of -3◦ . . . . . . . . . . . . . . . . . . . . . . . . . . 63
C.27 Palisupan pressure coefficient distribution back view of F20 c-foils obtained through
Paraview at an angle of yaw of 0◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
C.28 Palisupan pressure coefficient distribution front view of F20 c-foils obtained
through Paraview at an angle of yaw of 0◦ . . . . . . . . . . . . . . . . . . . . . . 64
C.29 NACRA F20 free surface deformation at 10 m/s . . . . . . . . . . . . . . . . . . 64
C.30 NACRA F20 free surface deformation at 5 m/s . . . . . . . . . . . . . . . . . . . 65
C.31 NACA 0012 nose CP distribution at 5 AoA with Palisupan . . . . . . . . . . . . 65
C.32 NACA 0012 top surface CP distribution at 5 AoA with Palisupan . . . . . . . . . 66
C.33 OpenFOAM lower surface CP for NACA0012 at 5 AoA . . . . . . . . . . . . . . 66
C.34 OpenFOAM top surface CP for NACA0012 at 5 AoA . . . . . . . . . . . . . . . 67
D.35 Screenshot from Battle of Boats - how technology won the America’s Cup . . . . 70

6
List of Tables
4.1 Particulars of the F20 NACRA geometry and run conditions . . . . . . . . . . . 23
4.2 Particulars of NACA 0012 horizontal foil geometry and run conditions . . . . . . 24
5.1 Principal dimension of the NACA0012 used in OpenFOAM . . . . . . . . . . . . 28
5.2 Lift and drag coefficient difference at varying mesh size including percentage
difference compared to the most refined mesh. . . . . . . . . . . . . . . . . . . . . 30
D.1 OpenFOAM results force breakdown computed for the free surface case using
LTSInterFoam for the range of depth analysed. . . . . . . . . . . . . . . . . . . . 68
D.2 Summary of boundary conditions for flow past hydrofoil OpenFOAM analysis . . 68
D.3 Summary of boundary conditions for free surface OpenFOAM analysis. . . . . . 69

7
Nomenclature
Symbol Description Unit

A Cross section area m2



AoA Angle of attack
AR Aspect Ratio -
c Chord length m
CD Drag coefficient -
Cf Skin friction coefficient -
CL Lift coefficient -
Cp Pressure coefficient -
CW Wave resistance coefficient -
D Drag force N
FX,Y,Z Force acting in x, y, z direction N
g Gravitational acceleration m s−2
k Turbulent kinetic energy kg m2 s−2
L Lift force N
l Length m
l/c Length to chord ratio -
M Mass kg
N Number of panels -
n Number of cells -
P Pressure Pa
Re Reynolds number -
Rw Wave resistance N
S Wetted surface area m2
t Time s
t Thickness m
V Inlet velocity m s−1
x, y, z Cartesian coordinates m
 Dissipation rate of turbulent kinetic energy m2 s−3
ν Kinematic viscosity m2 s−1
ρ Fluid density kg3 m−1
ζ Wave elevation m

8
1 Introduction
Hydrofoils are lifting surfaces that have been adopted for many decades in the marine industry
such as in sailing, defence and passenger craft. There are many available shapes for these
structures (L-foils, C-foils, horizontal foils) that are attached to the hull of a vessel in different
configurations (Figure A.3). As the forward speed of the craft increases they generate a lift force
that makes the hull rise out of the water. In this way wave and skin-friction drag are reduced,
enabling to achieve significantly higher speeds with the same amount of power supplied.
Recently this technology has been adopted by many sailing classes at professional racing
level. One of the reasons behind this is that in the past decade extensive research has been done
on foils applied to sailing boats. A substantial contribution came from CFD (Computational
Fluid Dynamics) becoming more accessible thanks to the significant progress of computing
power available nowadays.
It’s important to express that the application hydrofoil are destined is the first thing to
keep in mind while developing the design. Hydrofoils developed for sailing boats have particular
requirements and different operability conditions. Their range of speed varies from 10 to 20 knots
for dinghy classes up to almost 50 knots how the America’s Cup AC72 class has demonstrated
to be possible.
This project will use different codes to analyse the forces acting on foils having particular
concern for free surface effects, the interface region between air and water.

1.1 Background
The motivation of this project is to provide a combined approach for hydrofoil design. In present
days the wide availability of tools to analyse hydrofoil performance makes a composite approach
sensible for several reason:

• Computational costs can be reduced and thus time can be saved if the most appropriate
tool is chosen at the correct stage of the design spiral

• The complexity of physics makes specific tools more appropriate than others when ana-
lysing particular aspects of the design or the physics involved. For example surface panel
method is very good at capturing well lift force but doesn’t have the capacity to capture
forces arising due to viscous effects that can be analysed well with RANS codes.

Figure 1.1 shows the increased performance that hydrofoils can generate. Figure D.35 shows
how much foils can improve sailing performance if adapted properly. This was during the design
phase of the AC72, the catamaran class for the past America’s Cup.

9
Figure 1.1: Plot showing how hydrofoils aid in reducing drag at higher speeds

1.2 Aims and Objectives


The aim of this project is to understand the physics and performance of hydrofoil(s) interacting
with the fluid using different computational methods for a composite design approach.
In first place two codes, developed by University academics, will be used. These generate
results with a relatively small computational effort, being based on potential flow theory and
thinship theory. In parallel to this a second analysis will be made with an open-source code,
using RANS equations, that can achieve more accurate results at the expense of time and effort
to set up the computational case.
One key aspect that is considered is the free surface, which is the interface between water
and air. In this region, where the foil pierces the water or sits shallowly submerged under it,
additional resistance is encountered due to the wave pattern generated by the foil called wave
drag.
It is important to understand, at a concept design stage, how the geometry of the foils
will affect the overall performance. Specially in professional sailing, where every fraction of
time saved is essential, intensive research is carried out to optimise the design and enhance the
performance of these structures. To do this CFD in present days has become the most cost-
effective tool for developing a winning design because it has reached very powerful capabilities
of estimating and predicting accurately physical fluid flows.
The complexity and diversity of conditions that foils experience when operating on a sailing
boats make CFD the best tool given the ease of simulation of many different configurations and
cases.
F20 NACRA C-foils are used for the 3D analysis using Palisupan and Wavelstl (Fig. A.1).
Palisupan does not incorporate the free surface but can generate lift induced drag and pressure
forces. Different geometry variations are easily tested such as change in pitch/yaw angles.
Wavelstl is a Matlab code that computes the deformed free surface shape and resistance forces.

10
In OpenFOAM a 2D analysis is going to be performed on a NACA 0012 section profile.
Initially the flow over the foil at different angles of attack will be analysed capturing forces such
as lift and drag. Because of the time restriction of the project (one semester) a secondary aim
is to arrive at studying the foil submerged under the free surface at finite depths.
To conclude, a design methodology will be outlined stating how and when its most suitable
to use which code for analysing hydrofoil performance.

1.3 Hydrofoils and Foil Sailing Background


The conception of employing hydrofoils on hulls dates back to the early 1900s. Examples like
Italian inventor Forlanini’s ladder system (Figure 1.2a) employed to elevate the hull out of the
water or later American Alexander Graham Bell’s hydrofoil boat show how this technology isn’t
recent. The main difference between these early concepts and today’s evolution of foils is that
with the aid of tools like CFD and experimental testing the optimisation of the design is much
faster and of higher quality.

(a) Forlanini’s hydrofoil concept on Lago Maggiore, (b) Luna Rossa’s AC45 fully foiling 2014, Photo
Italy. Photograph originally published on ”Rivista copyright Carlo Borlenghi / Luna Rossa
TCI”, 1910

Figure 1.2: Early and recent developments of vessels adopting hydrofoil configuration

Lifting foils have revolutionised also sailing at a professional level. In the past edition of the
oldest sailing trophy in the world, the America’s Cup, thanks to a loophole in the rules, the
use of lifting foils was introduced to enhance performance of the catamaran class (Figure 1.2b).
Also the new AC62 design, which will be the class for the next edition, has been developed
around the ability of the catamarans to fully foil. Other sailing classes, both monohull and
multihull, have seen new rules to enable the use of foils to maximise their performance like the
Moth Class, the A-class catamaran or the NACRA class.

1.4 Development of Computational Fluid Dynamics


CFD is a computer based simulation of fluid flow that uses numerical methods to solve and
analyse the equations involved. The use of CFD for industry application started in the 1960s
for aerospace design and research. The limited computational power of the time restricted its
use only to high end R&D industrial products. The reason is that the extremely complex nature
of the solution was not affordable for many industries. With the availability of cost-effective

11
super-computers from the 1990s onwards and the development of fast and reliable CFD packages
interest in this method gained popularity (Versteeg and Malalasekera, 1995).
Different CFD methods are now available and are a balance between computational cost
and accuracy of results. A general subdivision, as described by Molland and Turnock (2007),
can be to group them in four different categories of ascending complexity:

1. lifting-line methods;

2. boundary element (surface panel) methods;

3. Reynolds-averaged Navier Stokes (RANS) methods;

4. Large Eddy Simulations (LES) and Direct Numerical Simulation (DNS) methods.

The first method is a model that predicts the lift distribution of a given 3D geometry using
mathematical approximations. The limitations of this method are that it is based on potential
flow theory so it does not incorporate viscous effects, it is limited to lifting surfaces of high aspect
ratio and neglects any flow effects from section camber of thickness (Molland and Turnock 2007).
The methods that will be used in this project are surface panel and RANS methods. The
former method has the characteristic of achieving good lift data with a comparatively small
computational effort. It’s a useful tool to use at a concept design stage to generate quickly a
sensible design. The latter method on the other hand, because of the increased computational
resources, is used to optimise the detailed design. This method will be used for the 2D analysis
of the project.

1.5 Risk Assessment and Ethical Approval


The nature of the project, which involves running simulations on a personal computer, does not
require any further risk assessment procedures as there is minimal to no risk in carrying out
this research project.
Furthermore, because no humans or animals will be involved in the project, there are no
ethical implications to be assessed.

1.6 Layout of Project


In the Hydrofoils chapter a general overview of hydrofoil theory will be given along with a selec-
tion of research papers that show what similar work has already been done. These papers will
mostly focus on the experimental testing of hydrofoils and what numerical analysis approaches
have been developed to analyse the flow around hydrofoils.
The following chapter will give a description about computational fluid dynamics; how it
works and the theory behind it.
The chapter on Potential Flow Theory will present the description of both Palisupan as well
as Wavelstl. A description of how they work and the methods adopted for the analysis will be
presented.
Next, the OpenFOAM chapter will present the basic functionality of the code and what
tools have been chosen in the computational analysis, the pre and post processing techniques.
The last chapter before concluding remarks and future work recommendations will present
and combine the results obtained from the three codes and how these can be utilised for a
cost-effective and optimised approach to hydrofoil design.

12
2 Hydrofoils
The analysis of flow past hydrofoils has been a relevant topic of research for many decades now.
It has been proved and demonstrated how the use of CFD has become a complement tool to
experimental testing in engineering practice, efficiently arriving to be almost a substitute to
expensive model testing (Zhang et al. 2006).
The first experiments were carried out in the 1950’s mainly by NASA at the Langley Research
Centre both in wind and water tunnel facilities (Kermeen 1956). Experiments were focused on
hydrodynamic and performance characteristic of foil sections in cavitating and sub-cavitating
flows. Acosta in 1973 wrote Hydrofoils and Hydrofoil Craft which is a paper that focuses on the
analysis of hydrofoil craft at the time, their performance focusing on cavitation and ventilation
problems. In the 1990’s with the significant progress of computational fluid dynamics hydrofoils
and craft configuration started to be investigated also with this method.
Research expanded also in developing tools to predict the behaviour and the motion of a
hydrofoil assisted boat numerically (Duncan and Matveev 2005). The focus was around the
general craft dynamics and motion prediction both in calm water and unsteady condition.

2.1 Theoretical Background

Figure 2.1: Schematic of physics of a hydrofoil submerged under the free surface. U is the inflow
velocity, c the chord of the foil, h the immersion depth, α the angle of attack and α0 the effective
zero angle of attack. As the flow of water travels around the hydrofoil two resultant forces re
produced which are Lift (L) and Drag (D). The coordinate system matches the OpenFOAM
axis system with Fx (negative direction) as the drag force and Fy as the lift force.

The purpose of lifting foils, explained earlier, is to produce a significant upward force that will
improve the performance of the vessel.

13
2.1.1 Lift Force
The lift (L) force acts perpendicular to the free stream flow, in the plane of the foil cross section.
In free stream, the lift generated by a two-dimensional foil is:
1
L= ρ V 2 A CL (1)
2
The lift coefficient CL is known to vary linearly with respect to the angle of attack. An approx-
imate formula to show the variation of CL is:

CL = 2π αT (2)

The angle αT is taken to be the AoA relative to the angle of zero lift. It is shown from formula
1 that the lift varies quadratically with velocity (V). This demonstrates that lifting surfaces
aren’t particularly suitable for low-speed craft but give distinctive advantages to crafts which
reach high speeds.
Physically lift can be expressed in different ways:

1. Difference in pressure: Water pressure on the lower surface of the foil is higher respect to
the upper surface. The difference in pressure creates the resultant upwards force. This
formulation comes from Bernoulli’s equation.

2. Newton’s 3rd Law: The hydrofoil pushes water downwards and this creates an opposite
upwards force on the hydrofoil.

Figure 2.2: Lift and drag force vectors for a hydrofoil catamaran

2.1.2 Drag Force


The drag (D) force instead acts in the same direction of the free stream flow. The general
equation for drag force is:
1
D = ρ V 2 A CD (3)
2
The drag force is considered to be dependent on may different factors and can be broken into
single components the most important of which are:

14
1. Viscous Drag: Arises due to the interaction between the body and the fluid, which has
viscous properties. There is a very thin region close to the body called boundary layer
that arises when a viscous fluid like water flows around a body. At the body the velocity
is equal to zero and gradually reaches the free stream velocity at a distance δ from the
body.
2. Induced Drag: Arises because of the lift force generated by the foil.
3. Wave and Spray Drag: This is generated by the creation of the wave pattern at the free
surface and the generation of spray (applies to surface piercing hydrofoils).

0.06

0.05
Drag coefficient CD

0.04

0.03

0.02
Induced drag Cdi
0.01
Profile drag Cd0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Angle of attack AoA [ ]

Figure 2.3: Breakdown of drag components

2.1.3 Parameters affecting Lift and Drag


The area close to the hydrofoil is a thin region called boundary layer where shear stresses are
present due to the viscosity of the fluid that flow around the body. Approaching the trailing
edge this layer thickens to about 2/3 % of the body length. Typically there are three regions
at the boundary layer:
1. Laminar flow region: the flow is smooth, orderly and steady
2. Transition region: the flow commences to break down
3. Turbulent region: flow becomes erratic with random motion and the boundary layer
thickens.
The Reynolds number, Re, is a non-dimensional quantity which determines the transition
between laminar and turbulent flow. It describes the nature of the flow based on the speed and
kinematic viscosity of the fluid and the length of the body:
l·v momentum
= Re = (4)
ν viscous
This quantity influences in different ways the forces acting on the foils. Studies made by Jones
and Whicker and Fehlner show that Re influences in a significant way the minimum drag
coefficient but very little the lift curve slope.

15
Another important geometry feature that affects in many ways the performance of the foil(s)
is its cross section. To maximise the performance of the foil its necessary to make an appropriate
decision regarding its section. This feature has a strong influence both on maximum lift and
stall, which occurs when at a further increase of AoA there is a loss of lift. Other than section
shape the stall angle is determined by the Re and from the operability conditions. The camber
of the foil section will also increase the lift. For symmetrical sections no lift is produced at 0◦
as opposed to cambered foils which produce lift.
Also the thickness of the section is an important dimension that can influence greatly a low
drag, high/lift ratio foil and delay cavitation (Molland and Turnock 2007).

2.1.4 Hydrofoil Design and Performance


By setting the computation to analyse consecutively different geometry variations (i.e. angle
of attack) it is possible to determine many performance properties of the foil as explained in
section 2.1.4. Lift to drag ratio is one of the key performance characteristics. The optimal angle
of attack where maximum trade off between lift respect to drag can be found. This will be the
most suitable condition of operation of the hydrofoil for a given profile section. Using lifting
line theory its possible to obtain the sectional drag coefficient Cd0 as well as the induced drag
coefficient Cdi which make up the total drag of the hydrofoil. In practical design drag needs to
minimised at very high speeds and for this reason this analysis is required.
Bal’s (2007) approach to analysing surface piercing hydrofoils was adopting a numerical
method based on Green’s theorem allowing the separation of the hydrofoil and the free surface
problem accounting for the interaction of the two in an iterative manner. Lift coefficient values
are calculated by the integration of pressure over the surface of the body. His results show how
the free surface affects substantially the characteristics and performance of foil(s) at or under
the free surface and thus declares the study of flow past hydrofoils of practical importance.
Two negative aspects that occur in viscous flows associated with high speed hydrofoil craft
are cavitation and ventilation. These phenomenon are similar and they influence negatively the
performance of hydrofoil reducing lift and increasing drag. Ventilation appears when the sub-
atmospheric pressure causes air to be drawn from above the free surface and forms a semi stable
cavity attached to the foil. This process has to be reduced to a minimum as it leads to loss in lift
and directional stability. The existence of a region of separated flow is necessary for ventilation
to occur for this gradual process to happen. Cavitation, which has similar physics involved, is
the formation of void or cavities in the liquid when the suction (i.e. negative pressure) increases
on a lifting device. For this reason during the design stage process its important to be able to
quantify the performance of surface piercing hydrofoils in every range of flow regime (Young
and Brizzolara 2013).
It is important to reiterate the importance of considering these two physical implications in
hydrofoil design but due to the nature of the project they won’t be analysed further because of
the increased complexity of simulation that would arise which is out of the scope of this project.
Experimental testing with complementary CFD was performed by Metcalfe et al. (2005)
on a surface piercing NACA 0024 foil. Their primary concern was about evaluating the wave-
induced boundary layer separation at the free surface. They show how at increasing Froude
numbers the separation increases and the influence of the free surface is seen at further depths
after which the foil surface pressure distribution behaves like 2D flow. The foil surface pressure
is high at the bow wave, low at the wave trough with a gradual rise towards the trailing edge.
The combination of experimental testing with CFD proved to be a valid approach to describe
with precision the flow physics.
The geometric characteristic of hydrofoils can be summarised as follows:

16
• Aspect ratio: The bigger the value is, the more the flow will resemble 2D performance (no
induced drag)

• Taper: the chord at hydrofoil tip is smaller than at the root

• Thickness ratio: The foil thickness varies over the span of the hydrofoil

• Sectional shape of airfoil may vary along the span

• Sweep: The tip of the hydrofoil may lie in front or behind the foil

As a summary the following are hydrofoil characteristics that are important to consider
when developing the design:

• Maximum lift coefficient and maximum lift angle of attack

• Lift curve slope

• Zero-lift angle

• Minimum drag coefficient

To compare the results of Palisupan with experimental data a conversion takes place to
process 3D Palisupan data into 2D data. This because when a hydrofoil of finite aspect ratio
is analysed the fluid flow happens also in the spanwise direction. This creates a vortex system
that ultimately induces a downward flow termed downwash (Molland and Turnock 2007). For
this reason what effectively happens is that the angle of attack is reduced by an angle  which
adds a component of drag called induced drag Di to the hydrofoil. So it is possible to consider
the hydrofoil of finite span at an angle α equivalent to a hydrofoil section of infinite span at an
angle of α -  with corresponding 2D values of CL and CD .
From Marine Rudders and Control Surfaces page 44 the conversion from 2D data to 3D
data is explained. For this case the opposite had to be made. First the 3D slope ai was found
using:

CL
ai =   (5)
57.3CL
α− πAR

after which the lift coefficient is obtained by:

CL = ai × α(AoA) (6)

17
3 Governing Principles of CFD
Computational fluid dynamics is a tool developed to provide a cost and time effective method,
alternative to expensive experimental testing, for analysing fluid flows. It is useful in optimising
the design and becomes almost necessary if scaling effects become problematic.
In CFD fluid is regarded as a continuous substance i.e. a continuum and its flow is caused
by externally applied forces. Properties of fluid such as density and viscosity are important
parameters to select appropriately. Their effects are very important close to walls (e.g. hydrofoil
surface).
The flow is analysed across very small regions called control volume where different assump-
tions are made:
dm
1. Conservation equation: =0
dt
d(mv) X
2. Conservation of momentum (Newton’s 2nd law of motion): = F
dt
Due to the very complex nature of the equations that describe the fluid flow (Navier-Stokes)
a discretisation method is applied which approximates the differential equation with a system
of algebraic equations solved by a computer applied to each control volume.
The following are components that are considered for a numerical solution method and the
specific one that was chosen for this project:
• Mathematical model: set of PDE and Boundary Conditions (BC’s)

• Discretisation method: Finite Volume Method

• Coordinate and basis vector system → Cartesian

• Numerical grid: 2D block structured grid with matching interface

• Finite approximations

• Solution method: iteration scheme

• Convergence criteria and geometrical progression: inner and outer iterations


The Navier-Stokes (N-S) equation are the governing equations that express the flow prop-
erties over a control volume. The equations apply Newton’s second law together with viscous
and pressure terms. However solution of these equations is still not possible even with the most
advanced supercomputers available.
The equations require time-averaging, idea first thought by Reynolds, to become the Reynolds-
averaged Navier-Stokes equation (RANS). The equation is used to describe turbulent flows and
are at the basis of most computational fluid dynamics packages available. The RANS equation
(7) along with the continuity equation (8) form the most popular turbulence model for CFD:
 
∂(uj ui ) ∂p ¯0 ¯0
ρ = ρf¯i − + − p̄δij + 2µS̄ij − ρui uj (7)
∂xj ∂xi
 
1 ∂ui ∂uj
where S̄ij = + is the main rate strain tensor, ū is the mean (time-averaged)
2 ∂xj ∂xi
0
component and u the fluctuating component.

∂ui
=0 (8)
∂xi

18
The appearance of the Reynolds stress makes the addition of a turbulence model necessary
to properly define and approximate the turbulent flow. In this project the k - ω SST model
developed by Menter (1994) was adopted to define the turbulence closure, recommended both by
Artur Lidtke and by documentation read about it (Menter 1994 and Mulvany et. al 2004). This
model embraces the best of two turbulent models: the k - ω at the close-wall in the boundary
layer and the k -  in the free-stream region.
The other numerical approach employed in this project to investigate lifting foils is with
a panel code. The principle of the panel method is based on the linear superposition of
source/sinks, vortices and/or doublet elements over a lifting surface, such that the boundary
conditions are satisfied on the body, across the wake and in the far field (Molland and Turnock
2007). The specific theory of the panel code used in the project, Palisupan, is outlined in section
3.1.
The idea of employing both panel and RANS codes to investigate forces on lifting surfaces
has been done before. Pittaluga and Becchi (2005) compared the calculations of flow past a
marine propeller between RANSE calculations and a panel method. By adopting viscous cor-
rection techniques to panel methods it is shown that the data is in good agreement between the
two codes. Describing panel methods “a fast and accurate tool for preliminary design” com-
pared to the more time consuming RANSE calculations they show how the viscous properties
can be evaluated either with boundary layer approximation techniques or with solution of the
RANS equations. In their RANS simulation the Morino Kutta Condition was preferred to the
Iterative Kutta Condition because the latter solver tries to minimise the numerical problem by
a correction factor affecting the solution with a non physical elaboration.

19
4 Potential Flow Theory Analysis
This study is based on a composite approach which takes into consideration a panel code and
a wave surface code. The two separate analysis will be combined together to make an overall
prediction of the performance of the hydrofoils investigated.

4.1 Palisupan
Palisupan is a code written by Professor S.R. Turnock that employs surface panel theory to
solve the equations for forces acting on a lifting body.

4.1.1 Theoretical Background


From Turnock (2000) and Molland and Turnock (2007) a summary of the theoretical back-
ground of the code is given:
The advantage of this panel code approach is that it can be used to model actual geometries.
The computational effort of this approach is much less compared to RANS method. It requires
only 1% of computational power respect to a RANS code. This because the computations are
carried out only on the body, wake and far field boundary surfaces and not through all the
fluid domain. The downside is that this method does not account directly for frictional drag,
separation and stall effects but employs numerical methods to do so.
In a lifting surface panel formulation the approximation of the full Navier-Stokes equations as-
sumes that the flow is inviscid, incompressible and irrotational, and satisfies Laplace’s potential
equation:

∇2 φ = 0 (9)
A later formulation was made by Lamb (1932) showing that a quantity satisfying Laplace’s
equation (Equation 9) can be written as an integral over the bounding surface S of a source
distribution per unit are σ and a normal dipole distribution per unit area µ distributed over S.
If v represents the disturbance velocity field due to the bounding surface and is defined as the
difference between the local velocity point at a point and that due to the free- stream velocity
then:

v = ∇φ (10)
where φ is defined as the disturbance potential. This can be expressed in terms of a surface
integral as:
ZZ " ! # ZZ !
1 ∂ 1 ∂ 1
φ= σ+ µ dS + µdS (11)
r ∂n r ∂n r
SB SW

Where SB is the surface body and SW a trailing wake sheet. r is the distance from the point
which the potential is being determined to the integration point on the surface and ∂/∂n is a
partial derivative in the outward normal direction to the local surface.
The conditions imposed on the disturbance potential are that (from Hess[6.8]):

1. the velocity potential satisfies Laplace’s equation everywhere outside of the body and
wake;

2. the disturbance potential due to the body vanishes at infinity;

20
3. the normal component of velocity is zero on the body surface;

4. the Kutta-Joukowsky condition of a finite velocity at the body trailing edge is satisfied.

5. the trailing wake sheet is a stream surface with equal pressure either side.

For a steady-state solution the wake dipole strength distribution is uniquely determined by
the application of the Kutta condition at the body trailing edge. As conditions (1) and (2) are
satisfied as functions of µ and σ, conditions (3) and (4) are used to determine µ and σ on the
body. The Kutta condition only applies at the trailing edge and some other relationship have
to be used to uniquely determine the distribution of µ and σ over the body. The numerical
resolution of this non- uniqueness is referred to as the singularity mix of the lifting-surface
method. One of the benefits of the lifting surface panel method is that it manages to compute
very efficiently the potential (or velocity) influence coefficients thanks to expressions derived by
Newman. The approach used is different from the original formulated by Hess and Smith but
is algebraically similar. It uses approximate expressions as the distance from the panel centroid
increases thus reducing the computational effort but nevertheless maintaining the accuracy of
the result.

4.1.2 Pre-processing
Palisupan requires the creation of a command file (.cmd) and an input geometry file (.pan) which
generates output result files. The geometry file consists of various parameters and information
regarding the conditions of the simulation such as flow velocity, panel size and scaling factors
and the geometry which is formed by sectional data points. The command file instead has
stored the rest of the information for a case to be run correctly. Information such as input and
output file location, number of cases to be run (i.e. AoA range) and convergence criteria is
found in this file.
Once the sectional points are obtained they have to be ordered correctly according to these
guidelines:

• The sections have to be ordered from the highest section going towards the lowest (i.e.
decreasing z value for each section)

• The order of the section has to go from the trailing edge towards the leading edge and
end at the same start point. Plotting the section on an x-y graph the proper order is
starting with the lower chord, pass through the nose and then the upper part finishing at
the trailing edge (Figure 4.1).

Figure 4.1: Order of the points in which each section has to be defines in Palisupan.

Once the .pan geometry file is defined and the .cmd file properly setup it is possible to run
the selected cases.
Palisupan outputs a number of files which are the following:

21
• .LOG file that contains each geometric variation where numerical data is stored. All
force and moment values are written into this file along with other information like Cp
distribution and panel size.

• .INP file that is possible to open with a visualisation software (Paraview is recommended).
Surface plots of velocity and CP distribution can be visualised as well as panel distribution
and check the geometry of the object.

• .xml file which has the non dimensional forces and moments stored for each case variation
including total number of panels and values of the case variation

The above diagram shows the setup of the case folder. In the analysis, the folder was located
conveniently in the C:\local directory of the university workstations.

Case Folder
.CMD file
Panel
Geometry
.PAN file
Results
.INP file
.LOG file
.xml file

To run the code another step is required which is to map the location of the Palisupan
application. The location is at \\soton.ac.uk\ude\courses\sesfsi\panels. Version 060519 of
Palisupan was selected. The command, from the terminal, to run the case will then be:

z:\exe\palisupan 060519.exe c:\local\<case name>.PAN

Different kind of errors may arise due to the improper set up of the .PAN file or the .CMD
file. If everything is properly defined then the solution should converge and the files mentioned
above should appear in the respective directories.
This concludes the pre-processing part for Palisupan.

4.1.3 NACRA F20


For the 3D analysis the NACRA F20 c-shaped hydrofoil was chosen. This because it reflects
an industry practical application as well as a project that the University is currently involved
with. It will be interesting to compare the results obtained from the analysis with Palisupan
and Wavelstl and compare them with the results obtained from the wind tunnel tests carried
out by PhD students of the Fluid Structure Interactions group.
The F20 NACRA foils are given in .STL file format. Manipulation is needed to obtain the
section data points to input into Palisupan because the code does not read or have the function
to import geometry files. To do this the geometry was first converted using Rhinoceros in DXF
format. Then this new file was opened with Solidworks, another CAD software available within
the University computers. After keeping only the section curves (horizontal) and deleting all
the vertical curves along the foil it was then possible through a macro to export the points
created on each section curve. A total of 28 sections were created, making the process very
time consuming. This approach was chosen to have the points in the correct order. A total

22
of 205 points were created for each section: 100 for each half curve, 1 for the leading edge
intersection between the two half curves and 4 for the truncated trailing edge. This particular
trailing edge is designed in this way mainly for manufacturing and safety purposes. The trailing
edge points were modified to create a sharp end which will give better results in the panel code
calculations. The wake therefore is initialised further downstream, the surface area is increased
because of more panels in that region increasing drag results. Overall the drag difference of
this modification can be negligible because the separation at the truncation trailing edge with
relative increase in drag is avoided.
The particulars of the geometry and the cases run for the NACRA F20 foils can be found
in Table 4.1.
The F20 foils will be analysed fully submerged at a range of yaw, heel and pitch angles.
Forces and moments will be investigated and presented.

F20 Nacra C-foils


Dimension Symbol Unit Value
Height h m 1.63
length l m 0.293

Speed V m s−1 10
Yaw angle φY AW degree -3,-2,-1,0,1,2,3
Pitch angle φP IT CH degree -3,-2,-1,0,1,2,3
Heel angle φHEEL degree 0,-2,-4,-6,-8,-10,-12

Panels spanwise Ns 50
Panels chordwise Nt 45
Wake panels I No Wake 50
Reflection plane m (0 0 1.63)
Reflection plane normal (0 0 -1)
Complete geometry file is found in Appendix

Table 4.1: Particulars of the F20 NACRA geometry and run conditions

4.1.4 NACA 0012 Hydrofoil


The NACA 0012 airfoil section was chosen due to its widespread application and research.
Many papers and sources were found with experimental and numerical data many of which were
gathered in the Turbulence Modelling Resource website page of the NASA Langley Research
Centre (website: http://turbmodels.larc.nasa.gov/naca0012 val.html). Duncan in 1981 carried
out towing tank experiments to analyse the free surface deformation and forces acting on a
2-dimensional hydrofoil moving horizontally. Pascarelli et al. in 2002 investigated the fluid flow
past a NACA0012 section positioned below the free surface solving the viscous, incompressible
N-S equations supplemented by the linearised dynamic and kinematic boundary conditions for
the free surface.
To model the geometry, data points were taken from airfoiltools.com, an internet website
that has a vast database of foil data. The file from the website is in Selig format which works
appropriately in Palisupan. The sections at y=0 and y=1 are created maintaining the same
data for the x and z axis. In this way a horizontal hydrofoil is reproduced with AR=8 and
results from this analysis will be compared with experimental data.

23
NACA 0012
Dimension Symbol Unit Value
Aspect Ratio AR 8
t/c t m 0.12
Chord c m 1

Speed V m/s 6
Angle of attack AoA degree 0-15 (1 deg increase)

Panels spanwise Ns 40
Panels chordwise Nt 30
Complete geometry file is found in Appendix

Table 4.2: Particulars of NACA 0012 horizontal foil geometry and run conditions

4.2 Wavelstl
This code is a Matlab app developed within the University of Southampton. The theoretical
background to this code is a numerical method to calculate the wave pattern and resistance of
body moving about the free surface. The idea is to compare the deformation of the free surface
and see how appropriate thinship is when analysing hydrofoilsl.

4.2.1 Theoretical Background


The potential, non lifting numerical method uses the linearised free surface condition for the
flow description. This method was originally developed for slender hulls, mainly of catamarans,
but given the geometry of the foil sections utilised (high Length:Breadth ratio) it is possible to
adopt it. The bodies are represented by planar arrays of the sources on the local hull centrelines.
According to Couser et al 1998 “the method presented offers orders of magnitude saving in
computational effort over the higher-order methods without sacrificing accuracy.”
The original method was formulated by Insel (1990).The hull is discretised into a large
number of quadrilateral panels; the source singularities are then placed adjacent to each panel
centre. The source strength, on the panel of the hull is proportional to the waterline slope. The
wave resistance of the sources is calculated from an expression derived by Insel which describes
the resistance in terms of far field Eggers coefficients for a source of finite channel. From this
original methodology several modification have been made which are listed below

1. Trim and sinkage

2. Transom stern effects

3. Hydrostatic correction

As a conclusion this method has been shown to provide excellent predictions of wave pattern
resistance for catamarans and monohulls with a variety of geometries (Courser et al 1998) and
is thus suitable for the application made of it in this project.

24
4.2.2 Pre-processing
The Matlab code is used to assess two different geometries. The study carried out is about
wave drag and associated free surface deformation both for the NACRA F20 foils and a NACA
0012 horizontal foil and how this is associated to immersion depth. The results will be used in
conjunction with Palisupan to compare and analyse the foil performance.

Wavelstl Analysis
NACRA F20 foil NACA 0012
1. Wave resistance at submergence depth 1. Wave resistance at varying depth
2. Yawed wave resistance

Wavelstl takes STL files in binary format as input geometry data. After importing the
geometry, from the menu there are several input parameters that can be modified which are:

- Fluid viscosity and form factor (1+k)

- Speed range (m/s) with number of interpolation speeds between minimum and maximum
value

- Tank dimensions (depth and width) in metres

The geometry tab gives the option to scale, rotate or translate the geometry of a given value.
To run the case the Calculate command needs to be selected. The output will be several
variables such as Wave resistance and relative coefficient, total resistance and coefficient and
viscous and friction resistance coefficients.
The other output is the free surface deformation. Before the calculation takes place input
speed, free surface area and mesh size of the free surface are defined in the dedicated menu.

4.2.3 NACRA F20 foils


The geometry used in the program is the same as Palisupan. Wavelstl reads binary .STL files
so a conversion from the ASCII file is required. The .STL file was opened with Rhinoceros and
exported as a BINARY .STL file. Following this quick procedure, the geometry tab is very
efficient in providing a quick way of changing angles or submergence depth.
As a first analysis the change in wave-resistance with respect to submergence depth was
evaluated for a range of speeds. This is a valid analysis to understand how the foils operate
at various submergence depths. The depth was incremented of 0.2 m each run ranging from a
depth of 0.0899 m to 1.2899m.
A second analysis was carried out with the foil varying angle of yaw between ± 5◦ . During
this analysis as well as the calculations of the code the free surface deformation was extracted.

4.2.4 NACA 0012 horizontal foil


To create the STL geometry for the code Solidworks was used. A curve through xyz points
was inserted into a sketch. The points for a NACA 0012 section were found on airfoiltools.com.
From the sketch the Extrude feature was used to define the horizontal foil. Solidworks has the
option to export a geometry in STL binary format so no further step was required to prepare
the geometry for Wavelstl.
The free surface deformation at different immersion depths was investigated for this foil
considering maximum wave elevation and aspect ratio of the foil.

25
5 RANS Analysis
OpenFOAM (Open Source Field Operation and Manipulation) is an open source CFD package.
Its the code chosen to carry out the 2D RANS analysis on the NACA 0012 sections. OpenFOAM
is composed of many solvers and utilities that are combined in applications to solve specific cases.
It has a wide range of application and can solve many different physical flows.
In this project a selection of these, specifically those to solve steady state turbulent flow
problems, will be employed. Like many other CFD packages, OpenFOAM comes with pre- and
post-processing environments that will help to set the case properly.
The mesh generation (pre-processing) will use the blockMesh utility in OpenFOAM com-
bined with a python script that will create a blockMeshDict.
For monitoring and checking convergence gnuplot and pyFoam were used. With these tools
its possible to constantly check the several variables that are calculated and if necessary to
change on-the-run the entries of the algorithm and solution control.
Paraview is the post-processing tool where graphical results are checked; like pressure and
velocity fields visualised within the domain.
The scope of using this code is to demonstrate the possibility to investigate at an increased
level of accuracy, respect to panel code, the forces acting on hydrofoils. The study will take
a geometry (NACA 0012 section) which performance is well known and prove the quality of
results from this code.

5.1 OpenFOAM Pre-processing & Guidelines


The pre-processing in OpenFOAM will ensure that the results obtained are accurate and reflect
the physical environment investigated. It is paramount to express the importance of a good
mesh. It is very important that time and attention is put into the mesh generation. Failure to
do so will inevitably make the solution either diverge or have no physical meaning.
The blockMeshDict file is required to create a mesh with the use of the blockMesh utility.
A precise definition of this file, written using keywords, is needed to make sure that the mesh
is created appropriately. Specific definition of this can be found in the OpenFOAM-2.3.0 User
Guide. In synthesis the file creates a series of hexahedral blocks, defined by 8 vertices, connected
to each other to form the whole mesh domain. For this reason strictly speaking OpenFOAM
does not deal with 2D geometries but it’s possible to do this by defining the third dimension
with a unit distance and defining as empty the patches in the third dimension.
The essential composition of a blockMeshDict, located in the \constant\polyMesh folder,
will include the following keywords:

1. convertToMeters: Scaling factor

2. vertices: vertex coordinates

3. edges: defines type of edges (i.e. arc or spline)

4. block : definition of the block including number of cells (grid refinement) in each direction
and expansion ratio

5. patches (boundary): choose the appropriate boundary property for each face

After creating the blockMeshDict, in the \0 folder the boundary conditions are set for each
face of the domain. Flow velocity, pressure and turbulent parameters have to be included in
this folder (Detail of this is found in tables D.3 and D.2).

26
.
The last folder is the \system folder. In here all the files that control the algorithm, discret-
isation schemes, time-step and convergence parameters are set.
Below is a tree diagram of the case folder and the most relevant files included in each
directory:

Case Folder
0
Include
frontAndBack
runConditions
constant
meshGen
polyMesh
blockMeshDict
RASProperties
transportProperties
system
controlDict
forceCoeffs
fvSchemes
fvSolution
Residuals
forceCoeffs

Another useful OpenFOAM utility which helps in the pre-processing environment is checkMesh
utility. This is a mesh tool that checks the validity of the mesh before the case is run. Several
geometry features (topology, skewness, smoothness and geometrical progression) are checked
and if any errors occur these will be written in files and can be visualised in Paraview. y+
(dimensionless wall distance) is another indicator of mesh quality. When using turbulent flow
solvers and y+ 6 1 it is possible to set a fixedValue boundary condition on the wing for nut
(kinematic turbulent viscosity) and k (turbulent kinetic energy).
The specific turbulence model selected is defined in constant\RASProperties. Another im-
portant parameter to set properly is the control solution of the analysis. This is done through
setting appropriately the discretisation schemes, residual control and under relaxation factors
which will improve both the accuracy of the solution and its convergence rate.

5.2 Flow around foil section


The first analysis was on a flow past a NACA 0012 foil section at varying angles of attack at
Re = 6 million. This foil section was chosen because of the extensive data available, making
comparison between the analysis and experimental testing possible.
An initial mesh with the Python code was created and run for a flow past the given airfoil
section. Table 5.1 contains the dimensions for the foil section used in the script to generate the
mesh.

27
NACA 0012 Sectional Geometry
Chord c m 1
Thickness t m 0.12
Angle of attack AoA deg 0 - 5, 10, 15

Table 5.1: Principal dimension of the NACA0012 used in OpenFOAM

The simpleFoam algorithm was used as solver for this flow. This code is used for steady-
state, incompressible and turbulent flows which is the best method to investigate the forces on
the hydrofoil section. Initially the one-equation Spalart-Allmaras turbulent model was used for
the solution of the transport equation. After further reading and documentation the k - ω SST
was chosen for its robustness and accuracy for this type of analysis.
It was a huge benefit to have given a mesh generation code, from postgraduate student Artur
Lidtke, that saved a large amount of time. The python code consists on a number of files that
write the blockMeshDict that is located in the constant\polyMesh folder of each OpenFOAM
case.
For the 2D flow past a foil section figure C.8 shows the principal geometry and dimensions
of the domain. Refinement regions (figure C.10) are found close to the foil boundary layer and
in the wake of the foil to determine with increased accuracy the forces. Around the foil there
are two regions of increasing mesh size. The arc inlet boundary (diameter = 15c) is a common
practice for this kind of flows as it facilitates the refinement region at the nose of the foil. The
domain extends 15c downstream and 7.5c for the top and bottom walls.
It goes without saying that starting the analysis with a well defined mesh saves a copious
amount of time. Avoiding to deal with boundary positions, cell size refinement regions and
improving the geometry to a satisfying level were very beneficial for the development of the
project.
The time saved could then be available to define precisely the boundary conditions to increase
the computation accuracy. Time was spent to find the best solution for the boundary definition
mainly at the foil wall and for inlet and outlet properties of the variables. For example slip for
top and bottom patches was chosen as it models the walls as inviscid.

5.2.1 Mesh grid independence analysis


Once a mesh is created, it is general practice to carry out a mesh independence study. This
requires to increase the mesh resolution of a factor until the converged solution becomes steady.
The mesh refinement defined under the meshDensity paragraph in the .py script was considered
to have a factor of 1 for the values given. The mesh independence analysis was carried out from
a mesh density factor of 0.5 up to 2.5 which corresponds to 1’391’250 cells.
This analysis was carried out at a foil angle of 5◦ and at Re = 6 million. The angle was
chosen so that there would be significant drag and lift values to better show the variation
compared to 0◦ where lift tends to zero and there is minimal drag. Higher angles could have
been investigated but stall for this section is seen around ∼15◦ so the chosen angle seemed a
better option.

28
0.56

0.54
CL

0.52

0.5
0 200’000 400’000 600’000 800’000 1’000’000 1’200’000 1’400’000
Number of cells of computational grid

Figure 5.1: Lift coefficient variation at increasing mesh refinement

0.018
0.016
0.014
CD

0.012
0.01
0.008
0 200’000 400’000 600’000 800’000 1’000’000 1’200’000 1’400’000
Number of cells of computational grid

Figure 5.2: Drag coefficient variation at increasing mesh refinement

Figures 5.1 and 5.2 show the results of the analysis. It is possible to notice the asymptotic
trend of lift and drag coefficient values tending to a constant value of ∼0.555 and ∼0.00885 for
lift and drag coefficients respectively.
This shows that a solution independent of the grid size will be achieved with a meshFactor
around 2.5 which corresponds to just under 1’400’000 cells. Due to the nature of the project,
the following analysis are all carried out with lower meshFactor values because of the limited
computing power available for the project. It was not possible to run OpenFOAM on the
workstation computers but it was installed on a personal computer through a partition of the
hard disk.
The results show that the increased mesh resolution converges towards experimental values
that are between 0.55 and 0.56 for lift coefficient and 0.0085 and 0.009 for drag coefficient.
The validation data was taken from http://turbmodels.larc.nasa.gov/naca0012 val sst.html, a
trusted source for experimental and theoretical data for NACA 0012 section.

29
meshFactor Average y + Cl Cd clError cdError
0.5 4.42351 0.5137 0.01672 8.293% 48.146%
0.6 3.65538 0.51146 0.01409 8.767% 38.467%
0.7 3.11417 0.51885 0.01275 7.218% 32.000%
0.8 2.71236 0.52729 0.01188 5.502% 27.020%
0.9 2.40227 0.53385 0.01125 4.205% 22.933%
1 2.15578 0.53875 0.01072 3.258% 19.123%
1.1 1.95511 0.54248 0.01031 2.548% 15.907%
1.2 1.7886 0.54532 0.01001 2.013% 13.387%
1.3 1.64822 0.5475 0.00975 1.607% 11.077%
1.4 1.52825 0.5492 0.00956 1.293% 9.310%
1.5 1.41917 0.55077 0.00938 1.004% 7.569%
1.6 1.32929 0.55199 0.00925 0.781% 6.270%
1.8 1.1836 0.55346 0.00905 0.513% 4.199%
2 1.06366 0.5547 0.00889 0.288% 2.475%
2.5 0.848634 0.5563 0.00867 0.000% 0.000%

Table 5.2: Lift and drag coefficient difference at varying mesh size including percentage differ-
ence compared to the most refined mesh.

5.2.2 Varying angle of attack analysis


Following the mesh independence analysis the forces on the wing at varying angle of attack was
investigated. This analysis is done to investigate many characteristics of the foil. The optimal
AoA of the foil can be found from the Lift/Drag ratio vs AoA curve. Stall characteristics can
be investigated which is important to know when developing a hydrofoil. Effects of Reynolds
number and minimum drag coefficient are also important variables to investigate. From a
practical point of view the most important consideration is to understand at what speed for a
given hydrofoil geometry and section enough force is generated to lift the boat.
In the present study, to reduce the computation time, the analysis was carried out with a
meshFactor of 0.5 and then results, according to the mesh independence analysis, were corrected
to more accurate values. This technique was adopted due to the numerous analysis that have
to be carried out in the project and the limited amount of time available. The results proved
to be a success with good data agreement.
The same grid geometry was used and a range of angles from 0 to 5, 10 and 15 degrees was
tested. The maximum Lift/Drag ratio was obtained at an angle of ≈ 3.33 obtained by fitting a
polynomial trend line through the points.

30
40
35.6976
30
L/D

20

10

0
0 1 2 3.328 4 5 10 15
Angle of attack

Figure 5.3: Lift/Drag ratio curve over angle of attack showing the maximum value at the given
angle obtained.

5.3 Free surface flow


5.3.1 Varying submergence depth
After the analysis of the flow past a foil section the free surface interface was added to the domain
and a study of the forces at various submergence depths was made. The volume of fluid (VOF)
method was adopted to capture the free surface interface. The analysis was carried out using
the steady-state solver of interFoam: LTSInterFoam (local time stepping) for 2 incompressible,
immiscible fluids using the VOF phase fraction based interface capturing approach (OpenFOAM
2.3.0 User Guide 2014).
As benchmark case towing tank experiments of a NACA0012 hydrofoil at various immersion
depth (Duncan 1980) were used and reproduced with small differences to compare this study
with the one done previously.

Mesh generation The first step was to take the python script used for the previous analysis
and modify this to accommodate the free surface. The main steps to do this were to:

1. Define a flat inlet boundary face from the circular inlet of the previous analysis

2. Extend the top boundary and insert a refined mesh region for the air-water interface

A first mesh didn’t bring accurate results when checkMesh utility was run. At the outer
nose refinement region, zero and negative volume cells would occur. This error is quite common
and won’t make the solver run because a floatation point error occurs. In order to avoid this
problem, after many unsuccessful variations of the mesh geometry, the second nose inlet region
arc was substituted by a line and a second region was added further upstream (Figures C.12
and C.13). In this way at the approaching node of negative volume cells the space was increased
and the mesh grading would adapt in a better way. This modification won’t affect the results
because going towards the nose of the foil there are further regions with a refined mesh.

Solver and Boundary Conditions The solver in use for the steady state free surface
flow requires some changes to the \system files and to the initial conditions \0.

31
- The ddtSchemes keyword in the fvSchemes file (\system) is changed from steadyState to
localEuler rDeltaT

- alpha.water and p rgh initial boundary conditions are created in the (\0) folder.

- A setFieldsDict is created in the (\system) folder.

The setFieldsDict is necessary to run the setFields utility that will define a phase fraction for
alpha.water and define the air-water interface with a value of 1 for the water phase and 0 for the
air phase. The distinction between the two fluids is defined in \constant\transportProperties
file.
The above steps are all necessary for the solver to work. If something has been overlooked
this will cause the solver to either diverge or crash. For this reason particular care must be
taken at this point to save time at a later stage trying to debug the errors from the code.
Initially runs were made adopting both the transient (interFoam) and the steady state
(LTSInterFoam) solvers. Once the case was set correctly the complete analysis made use of the
latter solver in order to obtain converged results both for the forces and for the free surface
deformation.

32
6 Results
6.1 NACA 0012 horizontal foil

6.1.1 Investigation of Panel size and distribution


Lift coefficient CL

0.385

0.38

0.375
0.0005 0.0006 0.0007 0.0008 0.0009 0.001 0.0011 0.0012 0.0013 0.0014
Inverse panel number 1/N
Drag coefficient CD

0.0131

0.013

0.0129

0.0128
0.0005 0.0006 0.0007 0.0008 0.0009 0.001 0.0011 0.0012 0.0013 0.0014
Inverse panel number 1/N

Figure 6.1: Lift [top] and Drag [bottom] coefficient panel number convergence for Palisupan

At first an analysis of the influence of panel number was made with Palisupan. This determines
the best combination of spanwise (Ns) and chordwise (Nt) panel number respect to force values.
Figure 6.1 shows lift and drag coefficient results plotted against the inverse of the number of
panels (1/N) for the NACA0012 horizontal foil. The best combination was found to be I Ns =
40 and I Nt= 30 for both forces and this was selected for the computation. For the chordwise
direction the number 3 panel distribution was selected. This refines the trailing and leading
edges to better approximate the nose curvature and the trailing edge sharp edge. A constant
distribution (number 0) was selected for the spanwise direction. This panel size will resemble
the shape of a square which will give the best results in Palisupan. After that the best panel
aspect ratio was found, it was kept constant but an increase in number of panel was set to find
out if there is convergence in the solution. Figure 6.1 shows that values begin to converge at
around 1/N=0.9 for lift meanwhile drag has a steady increase. For lift this tells us that at a
further increase in panel number the solution will remain stationary. The least number of panel
wants to be selected to cut down computation time.
In Paraview it is possible also to view the panel size and distribution which affects highly
the solution.

33
6.1.2 Angle of attack variation

OpenFOAM
Lift coefficient CL

1.5
Palisupan 2D
1 Ladson
Palisupan 3D
0.5

0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Angle of attack AoA [ ]

Figure 6.2: Lift coefficient at varying AoA at 6 m/s for NACA0012 section comparing Palisupan
and OpenFOAM results with experimental data by Ladson.
Drag coefficient CD

OpenFOAM
0.1 Palisupan 2D
Ladson
0.05 Palisupan 3D

0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Angle of attack AoA [ ]

Figure 6.3: Drag coefficient at varying AoA at 6 m/s for NACA0012 section comparing Palisupan
and OpenFOAM results with experimental data by Ladson.

Figures 6.2 and 6.3 compare the non dimensional force coefficients between the two computa-
tional methods and experimental data from Ladson. Generally there is a good agreement of the
data although some differences are noticeable. The lift force from Palisupan is lower than the
rest because data was obtained for a 3D horizontal foil (AR=8) and forces are expected to be
more significant when analysed in 2D respect to 3D analysis. For this reason a conversion was
made to obtain the equivalent 2D lift and drag values for Palisupan data. After the conversion
the data for 2D lift values matches quite well between the three sources.
OpenFOAM drag curve has a good agreement for the lower angles, up to 5◦ , and then starts
to increase respect to experimental data. This is due again to the fact that the computation was
steady state, which gives less accurate data for increasing angles of attack due to the unsteady
and unstable flow and increased turbulence that occurs.
C2
The total drag can be written as the sum of Cdo and CDi as: CD = CDo + ki L . Figure
AR
6.4 shows the plot of CL2 versus CD and from the equation fitted though the dots the value of
profile drag can be found. The overestimate of drag is seen also here for OpenFOAM that will
affect the results.

34
Drag coefficient CD
0.1
y = 0.0353x + 0.0078

0.05 Palisupan
OpenFOAM
y = 0.04521x + 0.00508
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
2
CL

Figure 6.4: CL 2 versus CD plot is useful in data analysis to extract profile drag Cd0

From OpenFOAM Figures C.15 to C.25 show how pressure and velocity distribution of
the flow field vary at an increasing angle of attack. This is possible because, differently from
Palisupan, OpenFOAM solves the equations along all the fluid domain and not just on the lifting
surface. At 0◦ , because the NACA 0012 is a symmetrical section, the plots are symmetrical
about the foil centreline. The velocity plot shows how there is an increase in velocity around
the foil which corresponds to a lower pressure seen in the pressure plot below. This trend is
noticeable in all the following graphs that are shown. There is however an important difference
in every velocity plot. That is the flow separation from the surface at an increasing length as
the AoA increases. This is due to the total pressure energy that decreases because of the viscous
stresses that arise which make the fluid velocity tend to zero earlier.
Pressure coefficient CP

−0.5

0.5 OpenFOAM
Gregory & O’Reilly
Palisupan
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x/c

Figure 6.5: Pressure coefficient distribution of foil’s upper surface along the chord for 0◦ AoA
comparing experimental results with OpenFOAM and Palisupan calculations.

35
Pressure coefficient CP

−4 Gregory & O’Reilly


Palisupan
−2 OpenFOAM

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1


x/c

Figure 6.6: Pressure coefficient distribution of foil’s upper and lower surface along the chord for
10◦ AoA comparing experimental results with OpenFOAM and Palisupan calculations.
Pressure coefficient CP

−10
Gregory & O’Reilly
−5 Palisupan
OpenFOAM

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1


x/c

Figure 6.7: Pressure coefficient distribution of foil’s upper and lower surface along the chord for
15◦ AoA comparing experimental results with OpenFOAM and Palisupan calculations.

Figures 6.5, 6.6 and 6.7 show the pressure coefficient distribution CP for angles of attack
of 0, 10 and 15 degrees. Comparison between experimental data (Gregory & O’Reilly) and
computational methods is made (OpenFOAM and Palisupan). At 0◦ AoA all three sources
have a very similar distribution with a slight decrease for Palisupan at the trailing edge. At
10◦ (Fig. 6.6) Palisupan and experimental data match well. OpenFOAM has a good agreement
for the lower surface plot but has a higher CP for the upper surface. The difference between
OpenFOAM and the other two data plots increases more for the 15 degrees analysis (Figure 6.7)
where Palisupan and experimental data still have good matching. The increasing difference is
given by the fact that the analysis in OpenFOAM was steady-state. For this reason as the angle
of attack increases the flow becomes increasingly unsteady making the steady state computation
less accurate. This reiterates the fact that flow past a foil can be treated as steady state only
for small angles of attack.
Pressure coefficient distributions can also be graphically viewed through Paraview, the post
processing tool used which reads both OpenFOAM and Palisupan data. Figures C.32, C.33
and show the pressure coefficient distribution surface plot from Palisupan and OpenFOAM.
Paraview is very useful when analysing and post-processing results. The graphical visualisation
makes possible to check if the geometry is defined correctly and that surface plots match what
its expected to happen.

36
6.1.3 Free surface study

0.02
Drag coefficient CD

0.018

0.016

0.014

0.012

0.01
1 1.5 2 2.5 3 3.5 4 4.5 5
d/c

Figure 6.8: Influence of foil depth from free surface on Drag coefficient studied with OpenFOAM

0.6
Lift coefficient CL

0.55

0.5
1 1.5 2 2.5 3 3.5 4 4.5 5
d/c

Figure 6.9: Influence of foil depth from free surface on Lift coefficient studied with OpenFOAM

The free surface was then added to OpenFOAM’s computational domain to study the deform-
ation caused by an approaching foil towards the air-water interface. Force values from the
previous analysis of the flow past the airfoil section were compared to this study to quantify
the influence of the free surface (Figures 6.8 and 6.9). Generally the performance of the foil is
affected negatively by the free surface as lift decreases and drag increases approaching the free
surface. The most significant change is seen in the lift coefficient which decreases of almost 0.1
from a depth of 5·c to a chord length depth. Drag coefficient values on the other hand don’t
increase significantly for the approaching foil. A similar trend for the resistance can be seen
when investigating with Wavelstl. As the foil approaches the free surface the percentage of total
resistance due to waves is very small.
Figure 6.10 shows the variation of pressure coefficient respect to the submergence of the
hydrofoil. Because the foil is kept at a constant angle of attack the general trend is very similar
with a variation only in the magnitude which is highest at the deepest position of the foil.

37
Pressure coefficient CP depth = 0.5c
depth = c
−1 depth = 1c
depth = 2c
depth = 4c
0 depth = 5c

1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x/c

Figure 6.10: Pressure coefficient distribution along the section foil at various submergence depth
investigated with OpenFOAM.

0.2

0
ζ/c

−0.2 depth = 5c
depth = 4c
depth = 3c
depth = 2c
−0.4 depth = c
depth = 0.5c
−14 −12 −10 −8 −6 −4 −2 0 2 4
x/c

Figure 6.11: Free surface plots at the immersion depths tested with LTSInterFoam (OpenFOAM
solver). There is a constant increase in the wave profile as the foil approaches the free surface.
The foil longitudinal position is at x/c=0 and the flow direction is from right to left.

Figure 6.11 compares the free surface deformation at the immersion depths studied with
OpenFOAM. An increase in surface elevation ζ is seen as the foil approaches the free surface.
From Figure B.4 the maximum surface amplitude is shown in detail for the OpenFOAM analysis.
Free surface results can be compared to Wavelstl free surface plots (Figure ??) and also
experimental data by Duncan (Fig. 6.13). The maximum wave amplitude from Wavelstl is of
a magnitude of less than half of the one computed from OpenFOAM. Wavelength comparison
shows as well disparate results. The position where the wave profile elevation becomes negative
is much further downstream for Wavelstl results.

38
0.04

0.02
ζ/c

0
depth = c
−0.02 depth = 1c
depth = 2c
−0.04 depth = 4c
15 14 13 12 11 10 9 8 7 6 5
depth 4= 5c 3 2 1 0
x/c

Figure 6.12: Wavelstl free surface deformation at a range of submergence depth for a horizontal
NACA0012 foil. The foil position is a x/c=0, the flow direction is from right to left.

When comparing the data with Duncan’s experimental results there is a big difference due
to the fact that the wave profile has repeated wave crests at small wavelength. This is not
experienced neither in OpenFOAM or in Wavelstl. The only similarity between OpenFOAM
and Duncan’s experiment is that the value of maximum wave elevation is similar. The difference
experienced with Wavelstl may arise because this code was developed for slender hulls, mainly
catamarans, and as such the conditions to use this for a submerged horizontal hydrofoil aren’t
valid.

0.1 Duncan
OpenFOAM
Wavelstl

0
ζ/c

−0.1

−15 −14 −13 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 0


x/c

Figure 6.13: Comparison of wave elevation taken from OpenFOAM, Wavelstl and experimental
values from Duncan. The foil position is for both three case at x/c=0

39
6.2 NACRA F20 hydrofoils

6.2.1 Investigation of Panel size and distribution

85
Drag force FX

80

75

70

0.0001 0.00015 0.0002 0.00025 0.0003 0.00035 0.0004 0.00045


Inverse of panel number 1/N

Figure 6.14: Drag coefficient panel number convergence for Palisupan

37
Drag force FZ

36

35

34

0.0001 0.00015 0.0002 0.00025 0.0003 0.00035 0.0004 0.00045


Inverse of panel number 1/N

Figure 6.15: Lift coefficient panel number convergence for Palisupan

For the NACRA F20 a panel size investigation was made similarly to the one done for the
NACA0012 horizontal foil. A ratio of I Nt= 50 and I Ns = 45 was found to give the best
aspect ratio. Following this figures 6.15 and 6.14 show the increased number of panels to see
if convergence was met. Here the trend is less linear and convergence can be seen just for the
drag force. For lift the values after some refinement tend to decrease steadily at increasing mesh
refinement. The changes however aren’t of a big degree for lift as opposed to drag where the
increasing refinement makes a significant improvement in FX .

6.2.2 Angle variation


Figures 6.16, 6.17 show the forces acting on the foil at the angles studied. Looking at figure
6.16 maximum drag force both for yaw and pitch is the highest 0◦ with a decrease approaching
± 3◦ . Lift however (Figure 6.17) starting from 3◦ has a linear increase up to -3◦ for both yaw
and pitch. It is interesting to notice how around -1◦ the foils start to generate positive lift force.
This is a particularly interesting performance criteria as this will be the condition necessary for
a foil to actually commence to create “positive” force.

40
Drag force FX [N] 160

140

Yaw
120 Pitch

−3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5



Angle variation θpitch,yaw [ ]

Figure 6.16: Drag force of NACRA F20 foils at a range of angles for yaw and pitch analysed
with Palisupan

Yaw
Lift force FZ [N]

Pitch
0

−100

−3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5



Angle variation θpitch,yaw [ ]

Figure 6.17: Lift force of NACRA F20 foils at a range of angles for yaw and pitch analysed with
Palisupan

0
Lift force FZ [N]

−10

−20

−11 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 0 1

Angle of heel θheel [ ]

Figure 6.18: Lift force of NACRA F20 foils at a range of angles of heel analysed with Palisupan

Heel was studied as well to simulate the heeling of the catamaran. Figure 6.18 show the
variation of heel from 0◦ to 10◦ . Although lift variation is constantly decreasing at larger angles
of heel the changes are very small.

41
6.2.3 Free surface deformation
In Wavelstl its possible to analyse the area and wake pattern of the surface piercing hydrofoil
condition. Wave energy and free surface deformation can be extracted from this software. By
analysing the surface deformation it’s possible to understand how much drag due to wave is
associated with foil. Figure 6.19 shows wave resistance values at corresponding depth of the
foil. The wave resistance has a non-linear increase with submergence depth. The resistance due
to wave can be quantified almost as negligible compared to the total resistance of the foil which
is 90.6773 N at a submergence depth of 1.2899 m.
A different assumption respect to a horizontal foil can be made when computing the deformed
free surface for the surface piercing F20 foils in Wavelstl. Because around the free surface this
geometry resembles a slender hull, high length to breadth ratio, the free surface may be better
approximated. The free surface deformation is presented in figures C.30 and C.30 but it is not
possible to compare this with an OpenFOAM study as only 2D studies were made.
Wave resistance Rw [N]

0.1
0.08
0.06
0.04
0.02
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4
d/c [m]

Figure 6.19: The plot shows the variation of the wave resistance respect to submergence depth..

6.3 Discussion
Once the general configuration of the foil is chosen, surface-piercing or fully submerged, a first
analysis is made with panel code. The reason for this is that panel codes can quickly predict
with good accuracy lift and lift induced drag. In this way it’s possible to test numerous designs
and assess their performance and select the most suitable. With Wavelstl, its possible to find
out how the hydrofoil will behave both below the free surface and when it pierces the water.
Once a solid design is selected with these two codes a 3D analysis with OpenFOAM can be
carried out which will give robust results. Although during the project it was noticed that some
of the data did not match with good accuracy, this is a marginal problem when it comes to
compare different designs between each other. This is the reason that initially a quick estimate
with panel codes proves to be more efficient. It is possible just to concentrate on the comparison
between different geometries and then extract more accurate data for the chosen design with
powerful codes such as OpenFOAM. It is important to specify that each code used captures
different physical aspects of the same structure and for this reason there will always be some
differences in results. These differences can be quantified if understanding of what part of the
physics is left out from each code.
Figure 6.20 provides the explained methodology idea for the practical hydrofoil design iter-
ation.

42
Weight Foil configuration Speed

Estimate Geometry

Update geometry

Palisupan Wavelstl

Results

Do results provide
no satisfactory values?

yes

OpenFOAM

End of design process

Figure 6.20: Diagram of the proposed methodology for hydrofoil design

43
7 Conclusion
The project has proven the feasibility of utilising CFD for a composite approach to assess several
aspects of hydrofoil performance and use results for practical hydrofoil design.
The main objectives posed at the beginning of the project have been met. Both 3D and 2D
studies have been carried out with the selected codes. This has made possible to formulate a
methodology to explain a logical sequence of work to begin a design spiral.
The project gave also a very good initial understanding of CFD techniques and the technical
background that comes with it. The main functionalities of the three codes have been used
arriving to an acceptable level of understanding for all three of them.
Results have proven to be partially in good agreement with experimental data. Palisupan is
able to capture quickly with good accuracy lift values. In Wavelstl the free surface deformation
and wave drag was found to be almost negligible. OpenFOAM can capture very well the viscous
forces and turbulence of the flow.
Differences arise due to the limitations of time and computational power available to carry
out further studies with better solution methods. Overall, after making the appropriate as-
sumption, the data obtained can be considered satisfactory when compared with experimental
testing. When developing a design it is also more important to capture well the trends of results
and not absolute force values. In this way it is possible to still make design assessments based
on the comparison of the data.
Many challenges were encountered throughout the project which were raised mainly due to
the time restriction of the project (one semester) and because completely new computational
environments had to be learned in very little time.

8 Recommendation for future work


Even though satisfactory conclusion can be drawn from the project there are still many elements
that can be improved or approached in further studied.
In first place a 3D analysis could be carried out in OpenFOAM to better compare and
quantify the accuracy of panel method compared to RANS codes and also better show the
trade off between computational time and accuracy of data. To better estimate the unsteady
characteristics that arise due to surface piercing hydrofoils a transient solution can be made.
In Palisupan the free surface was effectively modelled with the use of a reflection plane.
Further works can be made to find a method to effectively the free surface in the panel method
considering the hydrofoil as a lifting surface and the free surface as a non lifting surface.
In order to complete an analysis for a practical hydrofoil design a structural analysis of the
hydrofoil through finite element analysis (FEA) could be carried in parallel to CFD studies.
To conclude, there are numerous possibilities to further develop more complete hydrofoil
analysis using CFD but it is general good practice to set specific limits to one’s project to avoid
finish in increasing complex studies.

44
References
[1] I. H. Abbott and A. E. V. Doenhoff. Theory of Wing Sections. Dover Books, 1959.

[2] A. Acosta. Hydrofoils and hydrofoil craft. 1973.

[3] A. Alexander, J. Grogono, and D. Nigg. Hydrofoil Sailing. Juanita Kalerghi, 1972.

[4] S. Bal. High-speed submerged and surface piercing cavitating hydrofoils, including tandem
case. Technical report, Istanbul Technical University, 2007.

[5] M. P. Couser, D. J. Wellicome, and D. A. Molland. An improved method for the theoretical
prediction of the wave resistance of transom-stern hulls using a slender body approach.
Technical report, Department of Ship Science, University of Southampton, 1998.

[6] J. H. Duncan. The breaking and non-breaking wave resistance of a two-dimensional hy-
drofoil. pages 507–520. Journal of Fluid Mechanics, 1982.

[7] S. F. Hoerner. Fluid-Dynamic Drag. Sighard F. Hoerner, 1965.

[8] M. Insel. An Investigation into the Resistance Components of High Speed Displacement
Catamarans. PhD thesis, University of Southampton, 1990.

[9] J. Katz and A. Plotkin. Low-Speed Aerodynamics. Cambridge University Press, 2001.

[10] R. W. Kermeen. Water tunnel tests of naca4412 and walchner profile 7 hydrofoils in
noncavitating and cavitating flows. 1956.

[11] K. Matveev and R. Duncan. Development of the tool for predicting hydrofoil system
performance and simulating motion of hydrofoil-assisted boats.

[12] F. Menter. Two Equation Eddy-Viscosity Turbulence Models for Engineering Applications,
volume 32. AIAA Journal, August 1994.

[13] B. Metcalf, J. Longo, S. Ghosh, and F. Stern. Unsteady free-surface wave-induced


boundary-layer separation for a surface-piercing naca 0024 foil: Towing tank experiments.
Technical report, IIHR - Hydroscience and Engineering, The University of Iowa, 2005.

[14] A. F. Molland and S. R. Turnock. Marine Rudders and Control Surfaces. Elsevier, 2007.

[15] N. Mulvany, J. Tu, L. Chen, and B. Anderson. Assessment of two-equation turbulence


modelling for high Reynolds number hydrofoil flows, volume 45. International Journal for
Numerical Methods in Fluids, 2004.

[16] B. Pascarelli, G. Iaccarino, and M. Fatica. Toward the les of flow past a submerged
hydrofoil. 2002.

[17] C. Pittulaga and P. Becchi. Comparison between ranse calculations and panel method
results for the hydrodynamics analysis of marine propellers. Marine CFD 2005: 4th Inter-
national Conference on Marine Hydrodynamics, pages 65–77, 2005.

[18] Y. Semenov, G. Wu, and B. Yoon. A surface piercing body moving along the free surface.
Technical report, Department of Meshanical Engineering University College London and
School of Naval Architecture and Ocean Engineering University of Ulsan, 2012.

45
[19] S. R. Turnock. Palisupan User Guide. University of Southampton, 2nd edition, November
2000.

[20] H. K. Versteeg and W. Malalasekera. An introduction to Computational Fluid Dyanmics.


Longman Scientific and Technical, 1995.

[21] N. Xie and D. Vassalos. Performance analysis of 3D hydrofoil under free surface. Ocean
Engineering, 2006.

[22] Y. L. Young and S. Brizzolara. Numerical and Physical Investigation of a Surface-Piercing


Hydrofoil. Third International Symposium on Marine Propulsors, Launceston, Tasmania,
Australia, May 2013.

[23] Z. Zhi-rong, L. Hui, Z. Song-ping, and Z. Feng. Application of cfd in ship engineering
design practice and ship hydrodynamics. Technical report, School of Maths and Applied
Stats, University of Wollongong, 2006.

46
Appendices

47
A Hydrofoils

Figure A.1: F20 NACRA C-foil geometry and coordinate system

48
(a)

(b)

Figure A.2: (a) Shows a NACRA F20 sailing at an angle of heel with a good
view of the C-foil (Photo from http://www.nacrasailing.com/nacra-f20-carbon/) (b) Shows
the configuration of the foils and how they are attached to the hull (Photo from
http://www.icarussportsusa.com/boats/nacra-f20-carbon/)

49
Figure A.3: Shows the different configurations of hydrofoils: either
surface piercing (a) or fully submerged (b). Photo courtesy of
http://www.daviddarling.info/encyclopedia/H/hydrofoil.html

50
B Plots

0.16

0.14

0.12
depth = 5c
depth = 4c
0.1
depth = 3c
depth = 2c
0.08 depth = c
depth = 0.5c
ζ/c

0.06

0.04

0.02

−0.02

−0.04
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 4
x/c

Figure B.4: Detail of Figure 6.11, free surface plots in the region of maximum wave elevation.

150
Drag and Lift force FX,Z [N]

100

50

0 Drag
Lift
−50

−100

−3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5



Angle of yaw θyaw [ ]

Figure B.5: Drag force of NACRA F20 foils at a range of angles of yaw analysed with Palisupan

51
150
Lift and Drag force FZ,X [N]

100

50
Lift
0 Drag

−50

−100
−3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5

Angle of pitch θpitch [ ]

Figure B.6: Lift force of NACRA F20 foils at a range of angles of pitch analysed with Palisupan

100

90
Total resistance Rt [N]

80

70

60

50
−6 −5 −4 −3 −2 −1 0 1 2 3 4 5 6

Angle of heel θY AW [ ]

Figure B.7: Wavelstl total resistance experienced by the foil at a range of angles of yaw invest-
igated with Wavelstl

52
C Figures

Figure C.8: Figure showing the initial mesh created using a python script. It is possible to
notice the different regions of mesh refinement. These are made in areas where the flow has
the most significant changes. In this way the solution will be more accurate. Generally these
are around the foil boundary layer (wall) and at the trailing edge. The foil is positioned at x,y
(0,0). The flow direction is defined in the negative x direction.

53
54

Figure C.9: The mesh of the second analysis with OpenFOAM adding the free surface. The red region defines the water interface and the
blue region the air interface. The refinement regions are present around the foil like the first analysis and around the free surface interface
to capture with added accuracy the wave profile deformation.
Figure C.10: Detail of the mesh around the hydrofoil. In this zone the mesh is refined to better
approximate the boundary layer region of the foil. At the nose it is possible to notice an inner
and outer layer of mesh refinement..

Figure C.11: Division domain of the air-water interface with use of the setFIelds utility. This
is for a hydrofoil depth of 5c.

55
Figure C.12: Detail of region where the negative error cells occur before improving the geometry

Figure C.13: Improved geometry at the nose of the foil to avoid the negative volume cell error.

56
Figure C.14: Velocity magnitude field plot for Re = 6 million at 0◦ angle of attack.

Figure C.15: Pressure field plot for Re = 6 million at 0◦ angle of attack.

57
Figure C.16: Velocity magnitude field plot for Re = 6 million at 1◦ angle of attack.

Figure C.17: Pressure field plot for Re = 6 million at 1◦ angle of attack.

58
Figure C.18: Velocity magnitude field plot for Re = 6 million at 2◦ angle of attack.

Figure C.19: Pressure field plot for Re = 6 million at 2◦ angle of attack.

59
Figure C.20: Velocity magnitude field plot for Re = 6 million at 3◦ angle of attack.

Figure C.21: Pressure field plot for Re = 6 million at 3◦ angle of attack.

60
Figure C.22: Velocity magnitude field plot for Re = 6 million at 4◦ angle of attack.

Figure C.23: Pressure field plot for Re = 6 million at 4◦ angle of attack.

61
Figure C.24: Velocity magnitude field plot for Re = 6 million at 5◦ angle of attack.

Figure C.25: Pressure field plot for Re = 6 million at 5◦ angle of attack.

62
Figure C.26: Palisupan pressure coefficient distribution for the F20 c-foils obtained through
Paraview at an angle of yaw of -3◦

Figure C.27: Palisupan pressure coefficient distribution back view of F20 c-foils obtained
through Paraview at an angle of yaw of 0◦

63
Figure C.28: Palisupan pressure coefficient distribution front view of F20 c-foils obtained
through Paraview at an angle of yaw of 0◦

Figure C.29: Free surface deformation computed in Wavelstl for the F20 NACRA foil at 0◦
angle of yaw at 10 m/s at a submerged depth of 1.5m. The hydrofoil is positioned at x, y = (0,
-0.5) at the black dot.

64
Figure C.30: Free surface deformation computed in Wavelstl for the F20 NACRA foil at 0◦
angle of yaw at 5 m/s at a submerged depth of 1.5m. The hydrofoil is positioned at x, y = (0,
-0.5) at the white dot.

Figure C.31: Palisupan pressure coefficient distribution of the nose of the NACA0012 horizontal
foil, obtained through Paraview at an angle of attack of 5◦

65
Figure C.32: Palisupan pressure coefficient distribution of the top surface of the NACA0012
horizontal foil, obtained through Paraview at an angle of attack of 5◦

Figure C.33: OpenFOAM pressure coefficient distribution of the lower surface of the NACA0012
horizontal foil, obtained through Paraview at an angle of attack of 5◦

66
Figure C.34: OpenFOAM pressure coefficient distribution of the top surface of the NACA0012
horizontal foil, obtained through Paraview at an angle of attack of 5◦

67
D Files

Depth Pressure Viscous Porous


m N N N
0.5 3988.346 -2.54E-017 -20.5
1 4486.85 -3.04E-017 -21
2 4797.399 -3.58E-017 -21.5
3 4926.017 -3.7E-017 -21.6
4 5042.231 -3.74E-017 -22.3
5 5134.883 -3.73E-017 -22.2

Table D.1: OpenFOAM results force breakdown computed for the free surface case using LTSIn-
terFoam for the range of depth analysed.

inlet outlet wing


kqRWallFunction;
k freeStream freeStream
uniform 0;
calculated; calculated; nutkWallFunction;
nut
uniform 0; uniform 0; uniform 0
omegaWallFunction;
omega freeStream freeStream
uniform
p freestreamPressure freestreamPressure zeroGradient
fixedValue;
U freeStream freeStream
uniform (0 0 0);

Table D.2: Summary of boundary conditions for flow past hydrofoil OpenFOAM analysis. Front
and back patches are of type empty to provide the two dimensionally to the domain.

68
Variable inlet outlet wing top bottom
variableHeightFlowRate; inletOutlet;
alpha. fixedValue;
lowerBound 0; upperBound zeroGradient inletValue uniform zeroGradient
water uniform 0
1; value uniform 0 0; value uniform 0
inletOutlet; inletValue inletOutlet; inletValue
kqRWallFunction;
k turbulentKE; value turbulentKE; value slip slip
value turbulentK
turbulentKE turbulentKE
calculated; nutkWallFunction; calculated; uniform nutkWallFunction;
nut calculated; uniform 0;
uniform 0 uniform 0 0; uniform 0
inletOutlet; inletValue
fixedValue; omegaWallFunction;
69

omega turbulentOmega; value slip slip


turbulentOmega uniform
turbulentOmega
fixedFluxPressure; fixedFluxPressure; fixedFluxPressure;
p rgh fixedValue; uniform 0 totalPressure
uniform 0 uniform 0 uniform 0
p zeroGradient fixedValue; uniform 0 zeroGradient slip slip
pressureInletOut-
fixedValue; fixedValue; fixedValue;
U zeroGradient letVelocity; uniform
uniform flowVelocity uniform (0 0 0) uniform flowVelocity
flowVelocity

Table D.3: Summary of boundary conditions for free surface OpenFOAM analysis. Front and back patches are of type empty to provide the
two dimensionally to the domain.
View publication stats

70

Figure D.35: Screenshot from Battle of Boats - how technology won the America’s Cup explianing a plot to show and quantify how significant
the adoption of hydrofoils is to the performance of the AC72.

You might also like