Computational Methods For Hydrofoil Design - A Composite Analysis Using Panel Code and RANS
Computational Methods For Hydrofoil Design - A Composite Analysis Using Panel Code and RANS
Computational Methods For Hydrofoil Design - A Composite Analysis Using Panel Code and RANS
net/publication/271847208
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.
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
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
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
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
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.
(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.
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;
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.
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.
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
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 [ ]
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).
16
• Aspect ratio: The bigger the value is, the more the flow will resemble 2D performance (no
induced drag)
• Thickness ratio: The foil thickness varies over the span of the hydrofoil
• 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:
• Zero-lift angle
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
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)
• Finite approximations
∂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.
∇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;
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:
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.
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.
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
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.
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:
- Speed range (m/s) with number of interpolation speeds between minimum and maximum
value
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.
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.
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.
27
NACA 0012 Sectional Geometry
Chord c m 1
Thickness t m 0.12
Angle of attack AoA deg 0 - 5, 10, 15
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.
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
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
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.
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.
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.
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
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
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
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
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
85
Drag force FX
80
75
70
37
Drag force FZ
36
35
34
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 .
40
Drag force FX [N] 160
140
Yaw
120 Pitch
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
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
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.
44
References
[1] I. H. Abbott and A. E. V. Doenhoff. Theory of Wing Sections. Dover Books, 1959.
[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.
[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.
[14] A. F. Molland and S. R. Turnock. Marine Rudders and Control Surfaces. Elsevier, 2007.
[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.
[21] N. Xie and D. Vassalos. Performance analysis of 3D hydrofoil under free surface. Ocean
Engineering, 2006.
[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
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
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.
57
Figure C.16: Velocity magnitude 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.
59
Figure C.20: Velocity magnitude 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.
61
Figure C.24: Velocity magnitude 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
Table D.1: OpenFOAM results force breakdown computed for the free surface case using LTSIn-
terFoam for the range of depth analysed.
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
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.