MODELING HYDRODYNAMIC FORCES OF FLAPPING FINS FOR EFFICIENT
WATERCRAFT PROPULSION
CHAPTER 1
Introduction 1.1 General Overview
The concept of biomimetic flapping fins is inspired by observations of fish and cetaceans
(marine mammals like whales) in motion. Current studies on biomimetic fins have provided
valuable insights into their kinematics, demonstrating how these animals use their flapping tails
and fins to generate propulsive and maneuvering forces. To fully understand the mechanics behind
this force generation, it is necessary to integrate fin theory with modern fluid mechanics. The tails
of fish, which have a high aspect ratio, closely resemble the structure of fins, as illustrated in Fig.
1. In summary, these fins have been extensively studied through the development of kinematic and
dynamic models and by analyzing their hydrodynamic responses.
Figure 1: MIT laboratory robot "Robotuna" [1]
In general, a basic difference between a biomimetic thruster and a conventional propeller
is that the former absorbs its energy by two independent motions: the heaving motion and the
1
pitching (fin) motion, while for the propeller there is only rotational power feeding. In realistic sea
conditions, the ship undergoes a moderate or higher-amplitude oscillatory motion due to waves,
and the vertical ship motion could be exploited for providing one of the modes of
combined/complex oscillatory motion of a biomimetic propulsion system. At the same time, due
to waves, wind and other reasons, ship propulsion energy demand in rough sea is usually increased
well above the corresponding value in calm water for the same speed, especially in the case of
bow/quartering seas.
Numerous experimental and theoretical studies have been conducted on flapping fins.
Recent research and development efforts have shown that flapping fin systems, when optimized,
can generate high levels of thrust [2,3,4]. Additionally, environmental requirements and
intergovernmental regulations for ocean vehicles are becoming increasingly stringent. The Kyoto
Protocol emphasizes the reduction of pollution and environmental impact, particularly from ocean
vehicles, as crucial in combating global warming and climate change. For instance, the
environmental pollution caused by cargo ships worldwide has been identified as a significant
contributor to environmental degradation, primarily due to the use of poor-quality fuel in vessels
[5,6]. Satellite data further supports this concern, indicating that major maritime transportation
routes are characterized by areas with high concentrations of pollutants, largely due to emissions
from vessel engines [7]..
As aforementioned, the study of generation of thrust by flapping fins has been increased
significantly in recent years due to quest for clean energy generation by regulations. The principles
of flapping fins performance analysis is based on fluid mechanics especially hydrodynamics of fin
theory. In order to analyze thrust performance of flapping fins, it is getting approximated as fish-
like motion which is combination of harmonic heave and pitch motion as shown in Fig.2.
2
Figure 2: Considered laws of unsteady motion: A: combined translational-rotational oscillations, B: purely
translational oscillations, C: purely rotational oscillations, and D: advancing wave-type deformations [1]
The idea is that fish have their own thrust mechanism and they propel themselves in water
very efficiently due to rhythmic motion of their tail [8]. Propulsion is generated by means of the
combined motion of the tail. Thrust is manageable due to the fact that the combined translational
and rotational motions of the fin.
3
Figure 3: Vortex structures, forming behind the fin, performing heaving oscillations: A: near free surface, B:Near
solid flat ground and C: In unbounded fluid [1]
Flapping fins of a thin plate is approached as fish tails in steady forward motion. They are
thought of oscillating with combinations of harmonic heave and pitch motion during calculation
of their thrust performance. However, the hydrodynamics of oscillating fins or thin plates has been
mostly studied in experimental way because there is difficulty of analyzing vortex shedding
induced by boundary layer separation and estimation of nonlinear dynamics of the large vortices
generated by free layers as shown in Fig.3 [8].
Some may raise the question why many scientists have been interested in this field so far.
The answer is the useful advantages of flapping fin system [1] below;
• Can be considered as environmentally friendly
• Relatively low-frequency systems
• Sufficiently high efficiency systems
• Multi-functional being capable of operating in different regimes of motion
4
• Can provide static thrust
• Possess more acceptable cavitation characteristics than conventional propellers
• Can provide high maneuverability
1.2 Purpose of the Study
In this study, the data analysis of flapping fin oscillation which is combination of harmonic
heave and pitch motion is implemented for post-processing data from CFD code in order to
perform system identification for heave and surge force. The data results from CFD include 6
degrees of freedom outputs which are forces and moments. These forces are surge, sway and heave
forces and the moments are roll, pitch and yaw moments. In this work, especially heave and surge
force will be analyzed to develop system identification between inputs which is the combined
motions of flapping fin due to the forward motion of ship in waves and outputs which are heave
and surge force. Note that surge force is playing very important role to acquire thrust from flapping
fin and heave force is very important to get lift action. First of all, CFD software is very slow to
make analysis. If one parameter is changed for calculation, time needs to be dedicated. In addition,
CFD results are not a closed form relationship. Therefore, the relationship between force and
frequency, velocity or amplitude is unknown. As known, force is a function of the frequency
motion, the forward velocity of the vehicle and the amplitude of motion of the flapping fins. For
instance, once one parameter is changed, the effect of change of this parameter is unknown on
force due to unknown relationship. Furthermore, this parameter could be optimized for better
operations. To sum up, the objective of this thesis is to define a system identification based
computer algorithm through a set of equations which come from control system theory and can be
used to provide a prediction of relationships based on inputs and outputs.
5
CHAPTER 2
Review of Literature
This chapter includes a review of literature describing previous studies performed and the
fundamental information of flapping fins.
2.1 Previous Studies
The first study was implemented by Leonardo da Vinci to explain and apply the mechanism
of thrust generation by a flapping fin in 1490. Based on studies of Leonardo da Vinci, at the end
of the 19th century and the beginning of the 20th century, many studies were emerged in this field
to develop flight vehicles using flapping fins [9,10]. That time studies about bionics were limited
due to having inadequate scientific and engineering background.
The first explanation of the physics of flapping fins is given by Knoller and Betz [1] in 1909 and
1912 respectively. Both of them reached the conclusion that longitudinal thrust force and vertical
lift force is induced by flapping fin oscillations. However, an extensive research was done about
flapping fin by the end of the 19th century [1]. The first experimental study was implemented by
Katzmayr [11] to verify Knoller – Betz’s work in 1922.
Prandtl studied development of unsteady motion of a fin in incompressible flow and he
concluded that vortices are shed from a sharp trailing edge [1] in 1922. Birnbaum contributed a
linearized solution for Prandtl’s formulation and provided the resultant thrust force generated by a
flapping fin [1]. Also further studies of this unsteady airfoil theory were conducted by Wagner,
Kussner and Glauert. Keldysh and Lavrentiev provided a formulation for thrust generated by a
harmonically oscillating flat plate. The solution was obtained with conformal mapping method [1]
in 1935.
6
Golubev developed a flapping fin theory, different from Prandtl’s model, based on the
‘‘discrete’’ Karman form of the wake arrangement. Using the momentum theorem, he obtained an
integral equation whose solution allowed him to obtain the aerodynamic characteristics, including
the thrust of the flapping fin in the 1940s [1].
Polonsky and Bratt analyzed flow visualization experiments to verify von Karman and
Burgers’ observations. They illustrated the existence of different types of vortex structures behind
the oscillating airfoil [1].
Since 1960s, the investigations in this field have been intensified. Many researchers have
been conducting investigations about the aero-hydrodynamics of flapping fins especially
development of more extensive mathematical model [1].
The first symposium regarding bionics was held in Dayton, Ohio in 1960 [12]. Since that
day, there has been considerable progress in flapping fins.
7
2.2 The Investigation of Aerohydrobionts
Indeed, the first explanation about thrust generation by flapping fins is released in 1912
[1]. After that, the first study of analysis has been done in and after 1924 [13,14]. During 1960s,
bionics had been emerged in science as a separate field which covers studies of the basics of
flapping-fin propulsion [15]. Many publications in biology and biomechanics aimed to develop air
or water vehicles with flapping fins as well. The major purpose was the explanation of the
background of the efficient propulsion of fish, cetacean, birds, and insects. To sum up, after
understanding of the phenomenon of fish efficient propulsion, the desire of those studies is to
implement and apply that background to air or ocean vehicles with the principles of
hydromechanics mathematical model.
As mentioned earlier, many researchers were interested in flapping fin propulsors coming
from early observation of fish, insects and cetaceans who all utilize oscillating fin mechanisms for
thrust generation. It is useful to introduce the coefficient K which represents aero-hydrodynamic
mechanisms;
NL0
K
mU0
Where; N is the power of the system, m is the mass of system, U0 is the swimming or flying
speed, L0 is the maximum length of the object. The table below shows K coefficient for various
creatures and engineering system [1]
8
Table 1: Aerohydrodynamic characteristics of various systems [1]
Biological or technical
system
Coefficient of aerohyrodynamic perfection
K (kWs/ton)
Relative speed
U0/L0 (s-1)
Swimming in nature
Cateceans
4.40
0.3-1.0
High-speed fish
3.16
1.0-8.0
Dolphins
0.18
2.0-8.0
“Swimming” in technology
Submersibles 22.0
0.06-0.2
Flying in nature
Birds and insects
15.8
60.0-170.0
“Flying” in technology
Jet airplanes
81.0
4.0-9.0
As can be seen from the table, dolphins have the best hydrodynamic perfection in
comparison with the other systems. However, birds and insects have the highest relative speeds.
Many studies appeared about properties of insect flight, fish and cetaceans in 1920s. The research
was about the structure of the fins, flight modes, kinematics and deformations of the fins,
measurements of different characteristics, and simulations of the flight of insects [1]. Maxworthy
[16], Ellington [17,18], Freymuth [19], Liu [20,21] made research about aerodynamics of insect
flight. They realized that unsteady flow properties are very important for aerodynamics of insects.
9
Another study conducted by Schmidt-Nielsen [22] showing that it is more economical for
an animal to fly rather than to move on land. Table 2 illustrates the comparison between relative
speeds of motion for living creatures and man-made objects [23].
Table 2: Speed characteristic of some biological and technical systems [23]
Biological or technical system
Maximum speed (body lengths/s)
Human being
4
Cheetah
18
Supersonic aircraft
75
Starling (a type of bird)
120
Pershin [24] and Kozlov [25] have studied about fish and cetaceans to look into the
propulsive characteristics of these creatures. They systematized and analyzed bio-hydrodynamic
of properties of swimming. According to many studies about bio-hydrodynamics, there is a
correlation between the method of thrust generation, speed of displacement and typical regimes of
swimming for the creature. According to observations, low-speed marine creatures utilize wave-
type propulsors, like a whole part of body contributing for thrust generation. Therefore, body
flexibility is proportional to speed, i.e as increasing speed causes rise in flexibility of body.
High-speed sea creatures can be considered by a “system” that has three main functional
parts: a hull, a stem, and a fin. In follows, that “stem-fin” subsystem has 2 degrees of freedom.
The displacement of fin can be considered as heaving-pitching oscillations [1].
According to Pershin [24] and Kozlov [25], the analysis of statistical data of hydrobionts
10
(marine creatures) show that there is relationship between their geometric and kinematic
parameters. The most important parameters are the dimensional speed, frequency, and amplitude
of oscillation.
Many papers addressed kinematics of hydrobionts [24,26,27,28]. They showed that the
body, stem and fin have a relation with the swimming velocity and average period of oscillation.
Furthermore, the type of motion made by marine creature such as dolphins has been investigated
in literature. According to some papers, the motion of hydrobionts is sinusoidal type. However,
analysis of swimming motion based on underwater films of Cousteau, has shown that during
uniform translatory motion the oscillatory motion of the stem and the fin is close to a harmonic
motion.
The oscillations of fins are optimized in dolphins to avoid flow separation. Thus, kinematic
parameters of motion which are heave and pitch oscillations having have phase shift between them.
The pitch oscillations are following the heave oscillation by an angle close to ψ = π/2 [24].
To sum up, the main parameters of hydrobionts that can be as similar to human-made devices
can be [1];
• Dimensions, shape and platform of the fin
• Amplitude of the fin oscillations
• Frequency of oscillations
• Mass-rigidity characteristics of the fin (variable chordwise or spanwise)
• Large-fin deformations
11
2.3 The investigation of the flapping fin theory
The mathematical model of flapping fin theory is based on linearized inviscid
twodimensional (2-D) flow. Along with rise of computer technology, more complex 3-D fin
nonlinear inviscid and viscous approaches were introduced.
Keldysh and Lavrentiev implemented the linear theory to obtain thrust force acting on thin
oscillating plate. They used conformal mapping to solve the problem [1]. Furthermore,
Garrick provided a solution about same problem by using Theodorsen’s linearized solution for
incompressible unsteady flow past a flat plate.
Meticulous investigations of flapping fins propulsors were done by Khaskind, Sedov and
Nekrasov to make performance analysis for heaving, pitching and combined oscillations. The
forces acting on the oscillating thin plate were found by solving a singular integral equation [1].
Another analysis of combined heaving-pitching oscillating fin has been done by Gorelov
to look into the maximum thrust force achievement. He reached the conclusion that the maximum
thrust force is managed when there is a phase angle ψ=π/2 between heave and pitch
motions.
Gorelov was the person to develop the propulsive characteristics of a flapping fin in
nonlinear formulation. He obtained a solution by the method of discrete vortices and compared his
result with suction force obtained by linear theory.
Instead of the method of discrete vortices, panel method which is a nonlinear formulation
is used to investigate performance of the propulsive characteristics of flapping fins due to its
efficiency in comparison with the Navier - Stokes solver [29,30,31,32]. It was used to study the
12
effect of fin thickness, oscillation amplitude, phase shift, pitch axis location, and Strouhal number
[1].
In addition to these studies, many investigations have been conducted about mathematical
modeling of flow past flapping fin systems, two-dimensional flow models, three-dimensional flow
models and experimental investigations of flapping fins. Taking account of all these studies, the
application of flapping fins have been used in ship and offshore industry. Recently, many projects
are being conducted to improve water vehicles which utilize flapping fins.
Indeed, one of the first tests of full-scale vehicles with flapping fin propulsors was
implemented by Grebeshov in 1970s. He conducted the experiment a full-size cutter on hydro fins.
This cutter used hydro fins not only as lifting but also as propulsive components as shown Fig.4.
Figure 4: Cutter on underwater flapping fins [1]
The major work is recently going on in American, Canadian, British and German scientific
institutes. This work is published by Muller [33]. In addition, we have Russian and
Japanese scientific work about flapping fins development [34,35].
13
Figure 5: Russian ship with a fin wave energy extraction system [1]
Figure 6: Japanese ship with a fin wave energy extraction system [1]
Last but not least, other work about flapping fins development can be found in Refs.
[36,37,38,39,40,41,42,43,44,45]. Among these works, some of them are about biomechanics
aspects related to the structure of fins of birds and insects and the mechanics and dynamics of their
flight. The others are suggested aero-hydro-dynamics of flapping fins and to the practical
application of flapping fins vehicles as shown Fig.5 and Fig.6.
2.4 Kinematics of biomimetic fin thrusters
14
Figure 7: Description of motion of the flapping fin during CFD experiment [46]
In this section, it is indicated that the classical mechanisms which describes the motion of
flapping fins without consideration of the causes of motion applied by NTUA (National Technical
University of Athens) team (Dr. Kostas A. Belibassakis and PHD student Vasileios Tsarsitalidis).
As can be seen in Fig.7, the reference coordinate system is considered for flapping fins. Therefore,
the motion of body with flapping fins can be described in accordance to this coordinate system
such as translational (surge, heave and sway motions) and oscillatory (roll, pitch and yaw motions).
As shown in Fig.7 under constant parallel velocity, the flapping fin coordinate system is
experiencing a combination of both heave motion and pitch motion due to the waves. In this case,
two different basic frequencies have been considered for the motion of flapping fin, one of them
is relative frequency [47];
1 2f1
And the other frequency is flapping fin pitching frequency;
2 2f2 .
15
In the simple harmonic thrust producing case, relative frequency and flapping fin pitching
frequency is equal to each other [47].
1 2 And f1 f2 f
The surge motion of the flapping fin is shown as below
x(t) Ut
And t indicates time here. The vertical motion becomes;
h(t) h0 sin(2 ft)
Where; h0 is amplitude of vertical oscillation of the flapping fin. Simultaneously, the fin
undergoes a pitch oscillatory motion at a possibly different frequency yet as it is mentioned before
in the simple harmonic thrust producing case motions of the frequencies are equal. Therefore, pitch
oscillatory motion becomes [47];
(t) m 0 sin(2 ft )
Where θm is mean angle of attack, θ0 is amplitude of pitch oscillation of the flapping fin
and ψ is phase angle between the two movements.
2.5 Dynamics of biomimetic fin thrusters
The phase difference ψ between the two oscillatory motions is very important as far as the
efficiency of the thrust development by the flapping system is concerned. As it is discussed before
in the simple harmonic thrust producing case where: 1 2 , it usually takes value ψ = 900.
With the pivot point for the angular motion of the fin located around the 1/3 chord length from the
leading edge, a minimization of the required torque for pitching is achieved as shown Fig.8 [47].
16
Figure 8: Flapping fins simulation model in CFD environment
For flapping systems steadily advancing in unbounded liquid the main flow parameter
controlling the unsteady lift production mechanism is the Strouhal number [46];
St 2 fh0 /U ,
Note that the Reynolds number has a secondary role affecting viscous drag corrections. As
a result of the simultaneous heaving and pitching motions of the biomimetic fin the instantaneous
angle of attack is given by [47]:
(t) H (t) (t) tan 1(U1dh/dt) (t)
For relatively low amplitudes of purely harmonic motion and optimum phase difference
ψ=90o, the angle of attack becomes;
(t) (U1h0 0)cos( t)
17
which is equivalently achieved by setting the pitch angle θ(t) proportional to θH(t)
(t) wtan 1(U1dh/dt)
and thus
0 wh0/U
Where; w is termed the ‘pitch control parameter’ after [47], usually taking values in 0<w<1,
which is amenable to optimization. Decreasing the value of w, the maximum angle of attack is
reduced and the fin operates at lighter loads. On the contrary, by increasing the above parameter
the fin loading becomes higher and so is the chance of leading edge separation that would lead to
significant dynamic stall effects. Following [48], we exploit the above relation, as an active pitch
control rule of the flapping-fin thruster in the general polychromatic case, based on the time history
of vertical motion. In this case, the instantaneous angle of attack is [48]
(t) (1 w)tan 1(U1dh/dt)
2.6 Free Surface Effects
In the case of the biomimetic system under the calm or wavy free surface, additional
parameters enriches the above set, as the Froude number;
F U /(gL)1/2
where L denotes the characteristic (ship) length and g is gravitational acceleration, as well
as various frequency parameter(s) associated with the incoming wave, like μ = ω2L/g and τ = ωU/g,
τ has distinguish subcritical (τ <1 / 4) from supercritical (τ >1 / 4) condition [49].
18
2.7 The outline of CFD code by NTUA
The hydrodynamic analysis of ship and flapping fins are employed by NTUA team (Dr.
Kostas A. Belibassakis and PHD student Vasileios Tsarsitalidis). First, the hydrodynamic analysis
of ship is computed by developing algorithm in computational fluid dynamics environment. The
flapping fin is considered as appendage of a specific ship and by following the hydrodynamic
analysis of the ship, the hydrodynamic analysis of the flapping fin is employed as well. The
combined translational and pitch motion of the flapping fin is induced by the ship (49). Both
vertical and horizontal of arrangement of flapping fin is analyzed in this study (46,48,49) as shown
in Fig.9. The pitching motion of the fin about its pivot axis is selected properly in order to produce
thrust, with significant reduction of reduction of responses and generation of anti-rolling moment
by the vertical fin, useful for ship stabilization (49).
Figure 9: (a) Ship hull equipped with a horizontal flapping wing located below the keel, forward the mid-ship
section. (b) Same hull with a vertical flapping wing located below the keel, at mid-ship [49]
19
Ship hydrodynamic analysis has been applied by using linear theory using a Rankine
source-sink formulation and ship motions are calculated considering the additional forces and
moments because of unsteady propulsion systems.
Standard linearized sea-keeping analysis is used to achieve the motions and responses of
ship and flapping fin. The motion equations of ship (the coupled equation of heave and pitch
motion of the ship) can be calculated regarding the mass of ship, added mass and damping
coefficients (49).
A simplified lifting line model is applied for hydrodynamic analysis of flapping fin to
obtain expressions of the flapping fin forces (49). The generated flapping forces are considered as
two parts. One part is depending on the oscillatory ship amplitudes and the other part is dependent
on the incoming wave potential. The first part produces modifications of the hydrodynamic
coefficients of the system the other part adds on Froude-Krylov and diffraction forces in the right
hand side (49).
The horizontal arrangement of the flapping fin thruster and vertical flapping fin in
quartering and beam waves can be seen in detail (49). To sum up, the analysis of flapping fins
being appendage of ocean vehicles can be seen in detail (46,47,48,49).
20
2.8 Harvesting power and energy by flapping fins
Besides energy harvesting from application of propellers, flapping fins have a capability
of generating energy as well. We are getting energy from flapping fins induced by vortices,
freesurface waves and uniform currents. In the first case, at least two different methods are
considered to generate energy from the flapping fins; the constructive mode and the destructive
mode. In the constructive mode, the vortices created by the fin and the incoming vortex are in the
same phase and reinforce each other [50]. The destructive mode, on the other hand, is characterized
by a phase difference of approximately 1800 between the fin-generated vortices and the incoming
ones [50].
Related to these studies, generating flapping fin energy from free-surface waves has also
been experienced. According to some studies, submerged fins have the ability to propel itself
forward in sea conditions once it is right under a free surface.
According to some investigations, the coupling of different motion modes and external activation
is required in order to generate energy from flapping fins. For instance, one of modes of flapping
fin is acted as a periodic motion. Hydro-dynamically, this produces periodic variations in the lifting
and drag forces as well as pitching- rolling moments in an incoming flow. These timevarying
forces/moments in their turn can trigger other modes from which power extraction is achieved
through attached generators [50].
21
Figure 10: Schematics of a flow energy harvesting system based upon a flapping fin [50]
Two-dimensional fin is integrated with a damper c in incoming flow U as shown in
Fig.10. ρ is the fluid density. Here, chord length is 1m. The fin has a combined motion which is
heave and pitch as described before. Both motions are defined as harmonic motion. Heave motion
has already been described as
h(t) h0 sin(2 ft) And
also pitch motion has been defined as below;
(t) 0 cos(2 ft)
Inertia of the fin is ignored. The external moment is needed to trigger the pitching mode as
Me = -M. M is defined here as hydrodynamic pitching moment. The power input into the system
becomes;
Pi Me M
The power input is obtained by the damper c and it can be expressed as;
PO ch2
The mean power input becomes;
1
Pi t0T
Tt0 Pidt
And the mean power output becomes;
22
1
P0 t0T
Tt0 Podt
Where, T is defined period of the system. Finally, net power of the system becomes;
P Po Pi
In addition, the efficiency of system is analyzed with respect to net power of the system.
The power harvesting efficiency becomes
P
(1/2) U2YP
Where; YP is calculated as difference between the highest vertical position reached by the
leading and trailing edges and the lowest vertical position [50].
CHAPTER 3
Methodology
Flapping fins motion was previously describing the combination of heave and pitch motion
in chapter 2. All outputs which are forces and moments have been analyzed in CFD code based on
the input which is the motion of flapping fins. Therefore, we have discussed how to describe
physics behind flapping fin motion. In addition, in order to define system identification, the
correlation between input (motion of fin) and output (heave and surge force) needs to be
implemented. To do this, surge force (fx) and heave force (fy) depending on time needs to be
23
transformed in frequency domain. Hence, Fourier Transform, as well as the method used to
transform both forces in frequency domain will be discussed in this section. Furthermore, an
optimization method will be developed in MATLAB in order to establish approximating transfer
function model for heave force. Nonlinear model is implemented in order to find data points of
surge force signals.
3.1 Data Analysis for Heave Force
In this section, a brief presentation of Fourier analysis is given. What kind of data do we
get out of this model and why do we choose to work these data series. In addition to this, how do
we validate our Fourier methodology based on our data series?
In general, Fourier series can transform any periodic signal or function into harmonic
signals or sinusoidal functions. Therefore, all periodic functions can be analyzed easier by using
Fourier Transform. There are many reasons to utilize Fast Fourier Transform (FFT). However, the
fundamental idea is using FFT in many science fields to transform time-domain signals into
frequency-domain signals. This approach is very useful to define parameters of vibrating systems
[51].
The use of digital technology is on the rise in various applications because digital signal
processing has many advantages in comparison to analog signal processing. In this study, heave
and surge force data series are based upon discrete time signals. Therefore, FFT (Fast Fourier
Transform) has been applied to all heave force signals to transform into frequency domain by using
MATLAB. Unlike analog type of signals which are essentially continuous time signals, digital
technology encodes information using discrete-time signals as shown Fig.11 and Fig.12 [51].
24
Figure 11: Analog signal [51]
Figure 12: Discrete – time signal - low sampling rate [51]
If the variable to be represented requires fast transitions, it must be described by using a
higher sampling rate as shown in Fig.13 [51].
Figure 13: Discrete-time signal - high sampling rate [51]
In this study, the output which is heave and surge force data series is assumed to consist of
periodic functions. Therefore, any initial transient stage is ignored. For period –T≤ t ≤ T surge
force can be described as below.
25
y(t) c0 c2sin(2 ft )
Where; y(t) is the function (signal) in time domain and c0 and c2 are coefficients of series.
The output frequency is f and is phase angle.
With using Euler’s formula and Fourier integral, the transform of output function into
frequency domain becomes;
Y() y(t)ei2ftdt
Where ω is angular frequency and Y(ω) is a function of amplitude and phase spectrum of
output. This equation is called Fourier Transform of y(t). This analog Fourier Transform will be
needed to find data points of the system as well in the next section. However, due to the fact that
our data series of heave and surges force are composed of discrete signals (564 data), Fast Fourier
Transform methodology needs to be explained briefly. Based on the FFT algorithm, the expression
of frequency domain of heave and surge forces have been obtained in MATLAB for each data
series. First of all, the methodology of FFT will be explained and then how FFT algorithm is
implemented to surge force will be clarified.
The Fast Fourier Transform (FFT) is a very useful algorithm for Discrete Fourier Transform
(DFT). By using DFT computation time is decreased from N2 to Nlog2N where N is number of
samples of each surge force data series. Discretization of the time signal needed for
Discrete Fourier Transform is illustrated in Fig.14 [51].
26
Figure 14: Discretization of the time signal [51]
FFT algorithm is based on the fact that every discrete Fourier transform with N samples
can be divided into two Fourier transforms, each with N/2 samples (first with even samples and
second with odd samples) as shown in Fig.15 [51].
Figure 15: Dividing of the signal into the two new signals [51]
27
Fourier Transform becomes sum of two new smaller Fourier transforms:
N1 i2fr
Y f Yre N
r0
N N
2 1 i2f (2r) 2 1 i2f (2r1)
Y2re N Y(2r1)e N
r0 r0
N N
2 1 i 2fr i2f 2 1 i 2fr
Y2re N /2 e N Y(2r1)e N /2
r0 r0
Where r is sample number. We have two new Fourier transforms in equations above so that
they can be defined by real variables [51]
N
1 2fr
2 i
Af Y2re N /2
r0
N
1 2fr
2 i
Bf Y(2r1)e N /2
r0
And the complex variable becomes;
2
i
W e N
28
Because equations can be combined between each other, FFT equation which is used for
most of digital signal processing can be obtained as below [51];
Yf Af W f Bf
It is better to explain output data series structure before explaining how Fourier Transform
applied to surge force and heave force data series. For each selection of geometry, amplitude of
vertical oscillation of the flapping fin to chord ratio h/c, phase angle between heave motions and
pitch oscillation ψ, the mean angle of attack θm are created and all files
corresponding to this set are gathered. This typical set has simulations for five Strouhal numbers
ranging from 0.1 to 0.7 and the amplitude of pitch oscillation from 5 degree to maximum 75 degree.
For each run, the mean angle of attack θm is 0 degree and phase angle ψ is 900[46].
As shown in Appendix B, 564 output data series are obtained from CFD. Each column
includes iterations, time, surge force, heave force, sway force, roll moment, yaw moment and pitch
moment respectively [46]. Strouhal number ‘Str’ and pitch motion amplitude θ0 define each run,
‘its’ stands for the time steps used and time for the simulated duration. All forces and moments are
the mean values divided by water density ρ, pow stands for power (also divided by ρ), tra stands
for translational and rot stands for rotational., numbers 1,2,3 stands for the axes x,y,z and min,
max, dev stand for minimum, maximum and standard deviation values
Regarding geometrical parameter during CFD experiment, NACA 0012 standard fin is
used for all data series. As far as an individual flapping fin is concerned, the selection of platform
area, in conjunction with horizontal/vertical sweep and twist angles, and generating shapes ranging
from simple orthogonal or trapezoidal-like fins to fish-tail like forms, constitutes the set of the
most important geometrical parameters. Other important parameters are the fin aspect ratio, span-
29
wise distribution of chord, thickness and possibly camber of fin sections, as well as the specific
fin-sectional forms. Here, some of fin geometries regarding rectangular and fish type used for
experiment are presented in Fig.16, Fig.17 and Fig.18;
Figure 16: Rectangular Fin outline for s/c = 2, 4, 6 respectively [46]
Figure 17: Fish-like fin outline for s/c=4 and s = 150, 300, 450 respectively [46]
30
Figure 18: Cross section area for NACA 0012 [52]
As it was mentioned before two types of fins were used in data structure which is fish like
and rectangular fins. The geometric parameters are aspect ratio (AR) and skewback angle (such
fish15). For each selection of geometry, heave to chord ratio, phase angle, mean angle of attack
and position of pitching axis is created.
As known, 6 degrees of freedom 3 forces and 3 moments induced from the oscillating fins.
In this study, surge force and heave force are chosen to be analyzed. From the aspect of
hydrodynamics and control point of engineering, surge force is playing fundamental role to
generate thrust and heave force is playing important role for lift.
First of all, FFT (Fast Fourier Transform) needs to be calculated and illustrated in some
kind of framework in order to make data analysis for heave force data series. In MATLAB, an
algorithm is developed to transform heave force data into the Fourier domain. Due to the fact that
heave force time signals are discrete type of signal, Fast Fourier Transform, FFT is applied for all
heave force signal in order to compute dominant frequency. The excitation frequency is needed to
be computed in order to decide the system is whether linear or nonlinear and what kind of transfer
function model can be used to make data fitting for all heave force data series. To do this the
formula as below can be used to calculate the excitation frequency which is for heave and pitch
motion.
31
StU
fmotion
2h0
Where; fmotion becomes excitation frequency for each data, St is Strouhal number, U is flow
velocity and h0 is heave motion amplitude that can be calculated as below;
h
h0 c
h/c is heave to chord ratio which is defined for each data and c is chord length that is 1m
for all type of fins . Now, for all data FFT analysis results can be seen in Appendix D. The way of
obtaining heave force spectrums can be seen in Fig. 19 as well.
Figure 19: FFT flow chart
32
However, some of them that were identified as representative of important cases can be seen in
log-log plots for heave force data series in Fig.20, Fig.21, Fig.22, Fig.23, Fig.24 and Fig.25;
Figure 20: Heave force spectrum - frequency domain pattern for one type of data
33
Figure 21: Heave force spectrum - frequency domain pattern for one type of data
Figure 22: Heave force spectrum - frequency domain pattern for one type of data
Figure 23: Heave force spectrum - frequency domain pattern for one type of data
34
Figure 24: Heave force spectrum - frequency domain pattern for one type of data
Figure 25: Heave force spectrum - frequency domain pattern for one type of data
35
All plots generated in MATLAB have been observed to achieve data analysis and define
heave force system after Fast Fourier Transform for all data series. The first order frequency has
been become dominant frequency of all data spectrum which means the peak frequency is equal
to the excitation frequency which is the frequency of motions for almost each data series.
Therefore, heave force system is acting such as linear system. We can collect data points by
calculating the value of spectrum with respect to excitation frequency on the spectrum plots. All
data points have been calculated by developing codes in MATLAB. Now, the heave force system
is ready to develop describing function to build transfer function between inputs and outputs.
3.2 Describing Function for Heave Force
The data fitting concept achieved for heave force data series is explained in this section.
Transfer function will be generated between inputs and outputs for each amplitude to create
mathematical model for this system as shown in Fig.26.
Figure 26: Describing function heave force system
Since CFD data points have been gotten, second order transfer function is used for data
fitting with respect to different fish type geometries and different operational conditions depending
on heave and pitch motion amplitude. For suitable data fitting process, the sample second order
transfer function as below is used to compare to CFD data series for heave force system.
^ k(s z)
H(s)
(s 100)(s 200)
36
The purpose is to find the gain k and zero z along with suitable method for each specific
case depending on the geometry of fin and operational circumstances. Here, the poles in the
transfer function have been chosen arbitrarily because poles affect gain in this formula. If different
poles are defined, the gain will become different.
The optimization algorithm is developed in order to calculate gains and zeros effectively
in MATLAB. An optimization algorithm to make minimum error between approximation transfer
function and CFD data points has been developed in MATLAB. Thus, gains and zeros creating
small errors have been calculated by using suitable optimization method.
The objective function is defined as below which is calculating errors between transfer
function model and CFD data;
n ^ 2
ObjFunc log10(| Hi |) log10(| dpi |) 20
i1
Here gains and zeros are optimization parameters. In MATLAB, the optimization code is
created to calculate effectively these parameters for some specific fin geometry and operational
points. Unconstrained optimization minimization method is used to minimize the objective
function. In result section, the optimization result, gains and zeros can be seen for two geometries
and operational cases. Here the transfer function is compared with fish-like fin with skewback
angle 15 (fish15) with aspect ratio 4 (AR4) and skewback angle 15 (fish15) with aspect ratio 6
(AR6). The operational points 5 degrees, between 12-13.7 degrees, between 14.416.6 degrees,
between 20-23.7 degrees pitch amplitudes and 1 m, 1.5 m and 2 m heave amplitudes have been
investigated for data fitting.
37
All pass filter signal processing can be used to fix the phase of transfer function model.
The aim of all pass filter is to add phase shift (delay) to the response. An all pass filter is allowing
through all frequencies without changing the magnitude of transfer function model. Its magnitude
does not differentiate any frequency with all pass filter but it fixes our transfer function phase with
appropriate all pass filter design. The general transfer function of an all-pass filter can be seen as
below;
^^ c(n) c(n 1)s1 ... sn
H(s) 1 c(1) ... c(n)sn
Where; c are the coefficients here. They depend on the order of the all pass filter. c must
contain one, two, three, or four real elements.
For instance, in our case c with two elements generates a second order all pass filter.
Now, the fixed transfer function can be calculated by using an all-pass filter.
The all pass filter magnitude becomes;
^^
| H (i) | 1
The all pass filter phase becomes;
^^
H(i)
Now, the fixed transfer function model can be calculated based on the information. The
magnitude of transfer function model that found by using optimization algorithm does not alter
but the fixed phase of transfer function model becomes as below;
^ ^^
H(i) H(i) H(i)
38
3.3 Data Analysis for Surge Force
In MATLAB, code is developed in order to transform surge forces (fx) into the frequency
domain by implementing the FFT algorithm mentioned before. Due to the fact that surge force
time series of discrete type, Fast Fourier Transform is applied to all surge force signals in order to
determine the dominant frequency. Here, all data series is examined in order to define which the
dominant frequency that surge force signals is. Each data series has the same flow velocity which
is 2.3 m/s and the chord length of the fins is 1 m. In addition, Strouhal number, heave oscillating
amplitude and pitch oscillating amplitude are varying. First of all, the frequency of heave and pitch
motions f has been calculated by using the Strouhal formula as below;
StU
fmotion
2h0
After calculation of the frequency of motion, the FFT algorithm has been applied to all
surge force data series in order to obtain transform in frequency domain and semi-logarithmic
scale. Now, observations can be done for all surge forces. For instance, for some of the FFTs of
surge force data analysis, a pattern can be identified as shown in Fig.27, Fig.28, Fig.29, Fig.30,
Fig.31 and Fig. 32 semi-log plots;
Figure : Surge force spectrum - frequency domain pattern for one type of data
39
Figure 27: Surge force spectrum - frequency domain pattern for one type of data
Figure : Surge force spectrum - frequency domain pattern for one type of data
40
28
Figure : Surge force spectrum - frequency domain pattern for one type of data
41
Figure 29: Surge force spectrum - frequency domain pattern for one type of data
Figure : Surge force spectrum - frequency domain pattern for one type of data
42
30
Figure : Surge force spectrum - frequency domain pattern for one type of data
43
Figure 31: Surge force spectrum - frequency domain pattern for one type of data
Figure : Surge force spectrum - frequency domain pattern for one type of data
44
32
45
During the data analysis process, each surge force FFT set is observed in order to see
whether there is zero, second order harmonics. However, as can be seen in the analysis above, zero
and second harmonics are the dominant in comparison to other harmonics. Per the observation of
all output surge data series, motion frequency (combined heave and pitch oscillation) is half of the
dominant second order harmonic frequency of the FFT spectrum analysis for surge forces. Also,
the spectrum value which is a complex number can be calculated with respect to zero order and
second order frequency for all data series. Surge force system seems to be a nonlinear system based
on this information. Thus, the response’s Associated Transform will be used to make data analysis
for surge force.
In addition, in order to validate the FFT algorithm developed in MATLAB, comparison
was made between the period of some of surge data in time domain and the frequency of surge
force analyzed by FFT code in MATLAB.
1
f
T
Comparison for one data series is given in Fig.33 and Fig.34.
46
Figure 33: One of surge force in time domain for first data
Figure 34: FFT analysis spectrum - frequency for first data series
47
The dominant frequency in Fig.34 is 0.22 [1/s] which corresponds to a 4.54 [s] period and
the surge motion period in Fig.33 is close to 4.5 [s]. Using such comparisons, the FFT code can be
validated to work.
In summary, by using the FFT algorithm in MATLAB, surge force data series was
transformed from time domain to frequency domain. Thus, surge force can then be observed easily
and it is possible to observe whether surge force data series is linear or nonlinear system. As
mentioned earlier, most of surge force data outputs have double significant content of frequency
in comparison to the frequency of motion. Therefore, surge force will be analyzed as a nonlinear
system. Next step will be the explanation of the response’s Associated Transform and how to
implement for our surge force data series to obtain data points.
The inputs and the outputs of the flapping fin system have been defined before. It will now
be presented how to implement the response’s Associated Transform for input signals and output
signals. First of all, the nonlinear data analysis method will be explained with Associated
Transform. In addition to this, how the Associated Transform is used to compute data points of the
surge force system.
Here the surge force generation system can be considered as a black box along with
unknown nonlinear properties of flapping fin system. This black box is time invariant which means
that the properties of black box do not depend on time [53]. Signal y(t) is called the system’s
response and A(t) is called the virtual input signal that combines heave and pitch motion. The
nonlinear surge force system can be illustrated in Fig.35.
48
Figure 35: Nonlinear surge force system
The Associated Transform is used in order to calculate the data points of the system. The
Fourier Transform values of inputs and outputs have been used for computing transfer functions
of the system.
Before explaining how to obtain each data, it is better to provide some fundamental information about
response computation and the Associated Transform.
Associated Transform is needed in order to calculate data points. The Fourier Transform of
output signal Y(f) will be derived from Yn(f1,…,f2) by reducing all but one variables. Then a single
variable inverse Fourier Transform is necessary to compute y(t). However, in our process, Fourier
Transform is initially implemented for input and output signals. Then, the data points of the system
can be calculated given each point of Fourier Transform coming from CFD analysis.
The process is to calculate Y(f) is called the associated transform [54]. The
notation below is used to denote association of variables;
Y( f ) An[Yn ( f1,..., f2)]
In this study, n=2 since second order nonlinear system is considered. Now, the extension
to the general case can be made easily. For the Fourier Transform Ys(f1,f2), the associated transform
can be expressed as below;
49
Y( f ) A2[Y2( f1,, f2)]
Inverse Fourier Transform can be then expressed as below;
Y2( f 1, f2)ei2f1t1ei2f2t2df1df2 y2(t1,t2)
And along with setting t1=t2=t, the equation can be written as it is;
1 1 i2f1tdf1]ei2f2tdf2
y(t) [ 2 i Y2( f1, f2)e
2i
Now the variable of integration needs to be changed by writing
f f1 f2
1 1 i2( f f2)tdf ]ei2f2tdf2
y(t) 2i [ 2 i Y2( f f2, f2)e
All in all, once the order of integration is changed, the equation below is obtained;
Y2( f f2, f2)df2]e 2i
y(t) 1 [ 2 1i i2ftdf
It can be seen from the equation above, that term in the bracket should be equal to
50
Y( f ) F[y(t)]
And the proof is complete.
If the integral is rewritten, another formula for the association operation is obtained
1
Y ( f ) 2 i Y2( f1, f f1)df1
After providing basic information about the Associated Transform, the theory is
implemented for the CFD data series in MATLAB. First of all, the Fourier Transform is applied to
virtual input signals. This function is defined as “virtual fictitious” signal due to the fact that it
consists of heave motion and pitch motion. This combined input signal is created in order to
calculate the points of transfer functions of the system effectively. A(t) virtual input signal is
defined as below;
A(t) Asin(2 Ft)
Where; A is defined as virtual amplitude and defined as below;
A h02 02
Where; h0 is defined as the heave motion amplitude and θ0 as pitch motion amplitude. After
defining the virtual fictitious input signal, the Fourier Transform of the virtual input will be applied
for transform.
51
A( f ) A(t)ei2ftdt
Asin(2 Ft)ei2ftdt
Here, Euler’s theorem will be used for sine function;
cos(t) eit eit & sin(t) eit eit
2 2i
A( f ) A(ei2Ft
2ie i2Ft )ei2ftdt 2Ai
i2Ftei2ftdt ei2Ftei2ftdt]
[ e
By using Dirac’s delta function property, the virtual input signal in Fourier domain can be
obtained as below;
ei2Ftdt ( f F)
A
A( f ) [( f F)( f F)]
2i
52
This equation will be used for computing Fourier Transform of the output function with
Associated Transform.
Surge force data points can be calculated based on Associated Transform.
Surge force signals can be expressed before as below;
Y2( f1, f2) H( f1, f2)A( f1)A( f2)
Data points can be expressed as below;
Y2( f1, f2)
H( f1, f2)
A( f1)A( f2)
Along with the transfer function points relationship above, surge force can be expressed as
below;
1
Y( f ) 2 i Y2( f f2, f2)df2
Data points relationship can be substituted into surge force signal and the new expression between
input and output can be achieved along with Associated Transform.
f f1 f2
1
Y( f ) 2 i H( f f2, f2)A( f f2)A( f2)df2
This equation above shows the relationship between the combined motion of the flapping fin
and the surge force obtained by CFD code.
53
The input signal which is the motion of flapping fins can be found in Fourier domain to be;
A
A( f ) [( f F) ( f F)]
2i
The input signal can be expressed as below;
A
A( f1) [( f1 F) ( f1 F)]
2i
Along with Associated Transform, the equation can be expressed as below
f f1 f2
A
A( f f2) [( f f2 F) ( f f2 F)]
2i
Another input expression becomes;
A
A( f2) [( f2 F)( f2 F)]
2i
The correlation between the combined motion of the flapping fins and the surge force signal
has been found before;
1
Y( f ) 2 i H( f f2, f2)A( f f2)A( f2)df2
54
H represents data points of the surge force system. Now the input signals derived earlier will
be plugged into the relationship between input and output giving;
Y( f ) 21i H( f f2, f2) 2Ai [( f f2 F) ( f f2 F)] 2Ai [( f2 F) (
f2 F)] df2
After some algebra, the equation becomes;
( f f2 F)( f2 F)
8A2i H( f f2 Y ( f )
(( ff ff22
FF)) (( ff22 FF)) , f2)
( f f2 F)( f2 F)
And if more algebra is achieved, the equation becomes
A2
Y( f ) [
8i
H( f f2, f2)( f f2 F)( f2 F)df2
H( f f2, f2)( f f2 F)( f2 F)df2
H( f f2, f2)( f f2 F)( f2 F)df2
H( f f2, f2)( f f2 F)( f2 F)df2]
55
The response surge force equation has 4 terms; each term will be investigated one by one to
provide a simpler expression.
The first term in the equation becomes;
A2
Y( f ) 8 i H( f f2, f2)( f f2 F)( f2 F)df2...
This equation must be solved based on Dirac delta function property which is shown below;
( f2)df2 1
Therefore, once f2 = F is set for the first term of response surge force equation, the equation becomes;
A2
Y( f ) H( f F,F)( f 2F)...
8i
The second term in the equation becomes;
A2
Y( f ) 8 i [... H( f f2, f2)( f f2 F)( f2 F)df2...
Here, f2 should be set –F to solve the second factor, the equation becomes;
A2
Y( f ) [... H( f F,F)( f )...
8i
56
The third factor in the equation becomes;
A2
Y( f ) 8 i [... H( f f2, f2)( f f2 F)( f2 F)df2...
Here, f2 should be set F to solve the factor, the equation becomes;
A2
Y( f ) [... H( f F,F)( f )...
8i
The last factor in the equation becomes;
A2
Y( f ) 8 i [... H( f f2, f2 )( f f2 F)( f2 F)df2
Here, f2 should be set -F to solve the factor, the equation becomes;
A2
Y( f ) [...H( f F,F)( f 2F)]
8i
After these calculations, once all terms are arranged, surge force response spectrum becomes;
A2 H( f F,F)( f 2F) H( f F,F)( f )
Y( f ) 8i H( f F,F)( f ) H( f F,F)( f 2F)
Now, output signal which is surge force in Fourier domain can be calculated, once the
Dirac delta argument of each term is set to zero. Therefore, the equation becomes;
57
Y( f ) 8A2i HH(F(, FF), F()f( f2)F)H H( (FF,, FF) ) ((ff )
2F)
In order to find the data points of the system for surge force signals, surge force signals in
Fourier domain needs to be found by using CFD data. The output signal from CFD data has been
calculated in previous section as below;
y(t) c0 c2 sin(2 (2F)t )
Now this output signal will be calculated in order to find in Fourier domain as calculated before
for input signal in Fourier domain.
Y( f ) c0ei2ftdt c2 sin(2 (2F)t )ei2ftdt
Where; c0 and c2 coefficients in the equation.
By using Euler’s formula and Dirac delta properties, we are able to find output signal from CFD
data.
ei2Ftdt ( f F) c2
(ei(2 2Ft )
ei(2 2Ft ))ei2ftdt
Y( f ) c0( f )
2i
58
And then, the mathematical process can be applied to this equation and the equation becomes;
c22i ei(2 2Ft)eiei2ftdt e i(2 2Ft)eiei2ftdt
Y( f ) c0( f )
Y( f ) c0( f ) c2 eiei(2 2Ft)ei2ftdt c22i
eiei(2 2Ft)ei2ftdt 2i
Now by using Dirac delta function property, surge force signal from CFD data series in
Fourier domain can be computed as below;
Y( f ) c0( f ) c2 ei( f 2F) c2 ei( f 2F)
2i 2i
By using Associated Transform, surge force signal in Fourier domain has been calculated.
In addition to this, the surge force signal calculated from CFD data series has been used to compute
the surge force signal in Fourier domain. These equations represent the same mathematical
phenomenon. Now, the data points of the system can be computed.
Surge force signal in Fourier domain;
A2 A2 A2 A2
Y ( f ) 8i H(F, F)( f 2F) 8i H(F,F) 8i H(F, F)( f ) 8i
H(F,F)( f 2F)
The surge force signal in Fourier domain from CFD data
59
Y( f ) c0( f ) c2 ei( f 2F) c2 ei( f 2F)
2i 2i
Y( f ) b1( f ) b2ei( f 2F)b2ei( f 2F)
Here, b1 and b2 can be found from FFT analysis spectrum previously mentioned by developing code
in MATLAB.
As can be seen from the two equations, the coefficients of δ(f), δ(f-2F) and δ(f+2F) are the
same and data poimts can be calculated from the coefficients;
First of all, for δ(f-2F) coefficient, the equation can be written as below;
A2
b2 H(F,F)
8i
H(F,F) 8Ab22 i
In this equation, the coefficient b2 (for f=2F) is computed for each surge data in MATLAB.
In addition, A is defined as the virtual input signal amplitude developed for each data as well in
MATLAB.
A2 h02 Q02
h0 is defined as heave motion amplitude and Q0 as pitch motion amplitude.
For δ(f+2F) coefficient, the equation can be written as below;
60
A2
b2 H(F,F)
8i
H(F,F) 8bA22
i
The coefficient b2 (for f = 2F) and A are computed with the same method for each surge data
points in MATLAB.
For δ(f) coefficient, the equation can be written as below;
A2
H(F,F) H(F,F) b1
8i
Here, the coefficient b1 (for f = 0) is calculated from FFT data spectrum in MATLAB. It holds
that:
H(F,F) H(F,F)
If this property is plugged into the equation, it becomes;
H(F,F) H(F,F) b1
8A2
i
Therefore, the real part of H(-F,F) needs to be set to zero. The imaginary part can be found by
generating code in MATLAB as below;
8
Im(H(F,F)) b1 2A2
61
The imaginary part of transfer function H(F,-F) can be found as well.
Now, all data points of surge force have been computed. Each data point of the surge force
system can be seen in the chapter 4 in 3-D plots. The mathematical meaning of these results is
going to be discussed in the conclusion part.
CHAPTER 4
Results
This chapter shows and discusses all of the plots and results for flapping fins CFD data. The
methodology of the calculations using MATLAB was explained in the previous chapter.
First, the FFT spectrum analyses of all heave and surge force from CFD data is generated
in MATLAB. Then, all of these spectra were observed for suitable method in order to develop
mathematical model between inputs (the combined motion of the flapping fins) and outputs (the
surge or heave forces of the system). Transfer function model can be built for system identification
for heave force system because this system is pseudo-linear. Also, the codes generated in
MATLAB to compute data points of the surge force data series were based on nonlinearity. In this
chapter, the results are shown and discussed.
62
4.1 Comparison Transfer Function Model and CFD Data for Heave Force
As it was mentioned before, the comparison between transfer function model and CFD
data has been implemented for two types of fish-like fins which are fish15AR4 (fish type skewback
angle 15 and aspect ratio 4) and fishAR6 (fish type skewback angle 15 and aspect ratio 6). The
transfer function below has been used to generate mathematical model between input and output
for two types of fish-like fins.
Hˆ (s)
k(s z) (s 100)(s
200)
Where; k is the gain and z the zero for this type of system.
All gains and zeros for varying pitch motion amplitude are given in the following Table 3
- 8;
Table 3: fish15AR4, Heave Amplitude = 1m
Pitch Amplitude
Gain, k
Zero, z
θ0 = 50
1.5722x10+05
-9.8777x10-09
63
120≤ θ0 ≤ 13.70
1.1231x10+05
-3.4950x10-10
14.40≤ θ0 ≤ 16.60
7.7787x10+04
-1.0930x10-08
200≤ θ0 ≤ 23.70
9.9950x10+04
-5.5205x10-09
Table 4: fish15AR4, Heave Amplitude = 1.5 m
Pitch Amplitude
Gain, k
Zero, z
θ0 = 50
2.3998x10+05
4.8120x10-09
120≤ θ0 ≤ 13.70
1.7231x10+05
-5.1926x-09
14.40≤ θ0 ≤ 16.60
1.1178x+05
3.5734x10-09
200≤ θ0 ≤ 23.70
1.4936x10+05
-4.7219x10-09
Table 5: fish15AR4, Heave Amplitude = 2 m
Pitch Amplitude
Gain, k
Zero, z
θ0 = 50
3.2777x10+05
-1.1866x10-10
120≤ θ0 ≤ 13.70
2.3424ex10+05
-6.1904x10-09
14.40≤ θ0 ≤ 16.60
1.4729ex10+05
-6.4869x10-10
200≤ θ0 ≤ 23.70
2.0100ex10+05
4.6962x10-09
Table 6: fish15AR6, Heave Amplitude = 1 m
Pitch Amplitude
Gain, k
Zero, z
64
θ0 = 50
2.5481x10+05
-1.8217 x10-09
120≤ θ0 ≤ 13.70
1.8332 x10+05
7.1935 x10-09
14.40≤ θ0 ≤ 16.60
1.3114 x10+05
6.0017 x10-09
200≤ θ0 ≤ 23.70
1.6546 x10+05
-4.0756 x10-09
Table 7: fish15AR6, Heave Amplitude = 1.5 m
Pitch Amplitude
Gain, k
Zero, z
θ0 = 50
3.9681x10+05
1.4133 x10-09
120≤ θ0 ≤ 13.70
2.8419 x10+05
4.7066 x10-10
14.40≤ θ0 ≤ 16.60
1.8863 x10+05
-1.2517 x10-09
200≤ θ0 ≤ 23.70
2.4774 x10+05
7.3069 x10-10
Table 8: fish15AR6, Heave Amplitude = 2 m
Pitch Amplitude
Gain, k
Zero, z
θ0 = 50
5.4416x10+05
-2.2221 x10-09
120≤ θ0 ≤ 13.70
3.8906 x10+05
6.5648 x10-09
14.40≤ θ0 ≤ 16.60
2.4868 x10+05
5.1533 x10-09
200≤ θ0 ≤ 23.70
3.3470 x10+05
5.4166 x10-09
65
In the methodology section, it was mentioned that unconstrained optimization
minimization method is used to perform data fitting for different type of geometry and operational
points. The optimization results can be seen Table 9 - 14;
Table 9: Optimization result fish15AR4, Heave Amplitude = 1m
Pitch Amplitude
Obj. func. Initial Value
Obj. func. Initial Value
θ0 = 50
1.1600e+04
21.3995
120≤ θ0 ≤ 13.70
6.7871e+03
60.2657
14.40≤ θ0 ≤ 16.60
4.4287e+03
138.1106
200≤ θ0 ≤ 23.70
9.6731e+03
75.2356
Table 10: Optimization result fish15AR4, Heave Amplitude = 1.5 m
Pitch Amplitude
Obj. func. Initial Value
Obj. func. Initial Value
θ0 = 50
1.8148e+04
19.5714
120≤ θ0 ≤ 13.70
2.5081e+03
62.6374
14.40≤ θ0 ≤ 16.60
5.1996e+03
164.5401
200≤ θ0 ≤ 23.70
1.1431e+04
85.7753
Table 11: Optimization result fish15AR4, Heave Amplitude = 2 m
Pitch Amplitude
Obj. func. Initial Value
Obj. func. Initial Value
66
θ0 = 50
1.5205e+04
17.8378
120≤ θ0 ≤ 13.70
9.0496e+03
65.0498
14.40≤ θ0 ≤ 16.60
5.8203e+03
179.0440
200≤ θ0 ≤ 23.70
1.2824e+04
92.2964
Table 12: Optimization result fish15AR6, Heave Amplitude = 1 m
Pitch Amplitude
Obj. func. Initial Value
Obj. func. Initial Value
θ0 = 50
1.3915e+04
18.8520
120≤ θ0 ≤ 13.70
8.2494e+03
53.9990
14.40≤ θ0 ≤ 16.60
5.5020e+03
120.2066
200≤ θ0 ≤ 23.70
1.1879e+04
65.1946
Table 13: Optimization result fish15AR6, Heave Amplitude = 1.5 m
Pitch Amplitude
Obj. func. Initial Value
Obj. func. Initial Value
θ0 = 50
1.6222e+04
16.0768
120≤ θ0 ≤ 13.70
9.6898e+03
57.4684
14.40≤ θ0 ≤ 16.60
6.3642e+03
150.1683
200≤ θ0 ≤ 23.70
1.3832e+04
76.9731
Table 14: Optimization result fish15AR6, Heave Amplitude = 2 m
Pitch Amplitude
Obj. func. Initial Value
Obj. func. Initial Value
θ0 = 50
1.7978e+04
15.9085
120≤ θ0 ≤ 13.70
1.0794e+04
60.5454
67
14.40≤ θ0 ≤ 16.60
7.0549e+03
167.9887
200≤ θ0 ≤ 23.70
1.5382e+04
84.4968
Some results of the data fitting for heave force can be seen below in the plots. That follow
the process has been done for one type of flapping fin which is fish-like fin having 15 degree skew-
angle. Two fish like geometries have been investigated fish15AR4 and AR6 from the CFD data.
The transfer function model developed in MATLAB and CFD data can be seen in each plot for
every fish-like geometry and operational point.
For type fish15AR4 and operational point having 5 degrees pitch amplitude and 1 m heave
amplitude, it can be seen in Fig.36 how the transfer function model developed approaches the CFD
data. Each point represents one data point with respect to excitation frequency and the error can
be seen between the mathematical model and CFD data.
68
Figure 36: Comparison between CFD data and heave force transfer function model
Fish15AR4 has been investigated with specific operational points. However, since the
operational point is defined by the pitch motion amplitude, varying from 12 to 13.7 degrees
amplitude around 12 degrees pitch motion has been chosen to implement data fitting into CFD
data, along with 1.5 m heave motion amplitude. The data fitting result regarding this specific
operational point as shown in Fig.37.
69
Figure 37: Comparison between CFD data and heave force transfer function model
Fish15AR4 has been investigated with specific operational points. However, since the
operational point is defined by the pitch motion amplitude, varying from 14.4 to 16.6 degrees
amplitude around 15 degrees pitch motion has been chosen to implement data fitting into CFD
data, along with 2 m heave motion amplitude. The data fitting result regarding this specific
operational point as shown in Fig.38.
70
Figure 38: Comparison between CFD data and heave force transfer function model
Fish15AR6 has been investigated with specific operational points. However, since the
operational point is defined by the pitch motion amplitude, varying from 20 to 23.7 degrees
amplitude around 20 degrees pitch motion has been chosen to implement data fitting into CFD
data, along with 1 m heave motion amplitude. The data fitting result regarding this specific
operational point as shown in Fig.39.
71
Figure 39: Comparison between CFD data and heave force transfer function model
Fish15AR6 has been investigated with specific operational points. However, since the
operational point is defined by the pitch motion amplitude, varying from 14.4 to 16.6 degrees
amplitude around 15 degrees pitch motion has been chosen to implement data fitting into CFD
data, along with 1.5 m heave motion amplitude. The data fitting result regarding this specific
operational point as shown in Fig.40.
72
Figure 40: Comparison between CFD data and heave force transfer function model
For type fish15AR6 and operational point having 5 degrees pitch amplitude and 2 m heave
amplitude, it can be seen in Fig.41 how the transfer function model developed approaches the CFD
data. Each point represents one data point with respect to excitation frequency and the error can
be seen between the mathematical model and CFD data.
73
Figure 41: Comparison between CFD data and heave force transfer function model
4.2 Data points of Surge force
3D plots have been created for data points of surge force system. x-y plane is defined as
excitation frequency [1/s] plane and z axis is defined as magnitude or phase [degrees]. All data
points (fish15 and fish30) can be seen in plots.
74
Figure 42: Magnitudes of H(F,F)
The magnitude points of H(F,F) have been plotted to develop approximation model between inputs
and outputs as shown in Fig.42.
75
Figure 43: Phases of H(F,F)
All phase points of H(F,F) in degrees can be seen for surge force with respect to excitation
frequency as shown in Fig.43.
76
Figure 44: Magnitudes of H(-F,-F)
All magnitude points of H(-F,-F) can be seen for surge force with respect to excitation frequency
as shown in Fig.44.
77
Figure 45: Phases of H(-F,-F)
All phase points of H(-F,-F) in degrees can be seen for surge force with respect to excitation
frequency as shown in Fig.45.
78
Figure 46: Magnitudes of H(-F,F)
All magnitude points of H(-F,F) can be seen for surge force with respect to excitation frequency
as shown in Fig.46.
79
Figure 47: Phases of H(-F,F)
All phase points of H(-F,F) in degrees can be seen for surge force with respect to excitation
frequency as shown in Fig.47.
80
Figure 48: Magnitudes of H(F,-F)
All magnitude points of H(F,-F) can be seen for surge force with respect to excitation frequency
as shown in Fig.48.
81
Figure 49: Phases of H(F,-F)
All phase points of H(F,-F) in degrees can be seen for surge force with respect to excitation
frequency as shown in Fig.49.
82
CHAPTER 5
Conclusions
This study provides a closed-form phenomenological model specifically for heave force
that can produce the useful output signal without using CFD, if the specific conditions are met. In
addition, it is an identification method of nonlinear system generating e.g. heave force. Since the
data points of surge force are readily available, the study is open to expand to define system for
nonlinear surge force generation as well.
5.1 The Contribution of the Study
CFD software is normally slow, and expensive to buy. Also, the results are not in closed
form. For example, if one parameter’s value is changed, CFD run needs to be repeated. As it was
pointed out, the result of CFD is not in closed form relationship which means that the correlation
between force and motion, velocity or amplitude is unknown. Therefore, if one parameter changes,
the effect of the changing parameter is unknown. Since the mathematical system built between
inputs and outputs especially for heave force has a closed form, all these drawbacks are eliminated.
For each specific geometry and operational point, a transfer function is known. This comes with
mathematical advantage as well. We can look into all parameters insightfully and much more
efficiently.
In addition, since all data points are known for surge force response, this study can be
expanded to investigate in order to build a mathematical system as nonlinear black box for this
variable as well.
83
5.2 Future Work
A potential future study based on this work would include system identification for surge
force. Due to knowledge of all necessary information of surge force system, it can be expanded to
develop nonlinear transfer function (describing function) model by using appropriate
methodology. For this study, a Volterra model could be a suitable method to obtain mathematical
model that properly predicts the hydrodynamic response of flapping fins. Then such model can be
used to make some analysis without costly re-execution of the CFD code [55], since surge force
can be predicted by a nonlinear black box system. Volterra theory of nonlinear systems provides a
mathematically rigorous approximation technique to describe these unsteady hydrodynamic
effects.
Another track of future work might aim to investigate different kind of geometries like
other fish or rectangular type of fins. In addition to this, the mathematical system can be defined
for other degrees of freedom such as sway force.
84
Appendix
Appendix A: CFD folder tree made spectrum analysis
ay
pana
are
‘TITLE
‘VARIABLE:
ZONE
T="
an
a2
13
14
45
16
7
ae
19
20
21
22
23
24
28
26
27
28
23
30
31
32
33
34
35
36
37
38
39
40
41
42
43
“4
45
46
47
48
49
50
en
ee
Cae
weep
rete
*
system
Forces
and Moments
divided
by
rho
-
Ictal”
S
"ite",
"cime’
systen
",T-
0.231884E+00
0.347826E+00
0.463768E+00
0.579710E+00,
0..695652E+00
0.811594E+00
0.927536E+00
0.104348E+02
0.115942E+01,
0.127536E+01
0.139130E+01
0.150725E+01
0.162319E+01
0.173913E+01
0.185507E+01,
0.197101E+01
0.208696E+01
0.220290E+01
0.231884E+01
0.243478E+01
0.255072E+01
0.266667E+01
0.278261E+01
0.
2898SSE+01
0.301449E+02
0.313043E+02
0.324638E+01
0.336232E+01
0.347826E+01
0.359420E+01
0.371014E+01
0.382609E+01
0.394203E+01
0.405797E+01
0.417391E+01
0.428986E+01
0..440580E+01,
0.452174E+01
0.463768E+01
0.475362E+02
0.486957E+01
0.
498SS1E+02
0.510145E+02
0.521739E+01
0.533333E+01
0.544928E+01
0.556522E+01
0.
568116E+01
0.
579710401
Ex",
156,
J:
~0.122912E+01
~0.124216E+01
-0.124710E+01,
-0.122012E+01,
-0.116429E+01
-0.108637E+01,
-0.989541E+00
-0.879067E+00
~0.759705E+00
~0.635930E+00
~0.512115E+00
-0.392391E+00,
~0.280607E+00
~0.180450E+00
~0.949333E-01,
-0.266957E-01
0.219873E-01
0.496563E-01,
0.554453E-01
0.389048E-01
0.55332SE-03
-0.58716SE-01
-0.137353E+00
-0.233006E+00
-0.343335E+00
~0.465067E+00
-0.594855E+00,
~0.729176+00
-0.864102+00
-0.995293E+00,
~0.112126+01,
~0.123602E+01,
-0.133735E+01
~0.142217E+01
-0.148836E+01,
-0.153390E+01,
-0.155758E+01
-0.155270E+01,
-0.153711E+01
-0.149342E+01
-0.142935E+01
~0.134601E+01
~0.124631+01,
~0.113273E+01
-0.100856E+01
-0.877259+00
~0.742617E+00
~0.608016E+00
~0.477677E+00
-0.789712E+01
~0.716912E+01
~0.724108E+01,
~0.717S31E+01,
~0.701803E+01
-0.678359E+01
-0.648150E+01
-0.611813E+01
~0.570023E+01,
~0.523425E+01,
~0.472558E+01,
~0.417909E+01,
~9.359998E+01,
~0.299374E+01
~0.236415E+01
~0.171690E+01,
~9.105651E+01,
~0.387532E+00
0.2e5399E+00
0.957262E+00
0.162391E+01
0.228082E+01
0.292345E+01
0.384690E+01
0.414825E+01
0.472213E+01
0.52645eE+01
0.577160E+01
0.623968E+01
0.66643E+01
0.704503E+01
0.737424E+01
0.765373E+01
0.767762E+01
0.204741E+01
0.81590SE+01
0.821215E+01
0.220726E+01
0.814242E+01
0.801932E+01
0.78402SE+01
0.760366E+01
0.731424E+01
0.697244E+01
0.65273E+01
0.614538E+01
0.566786E+01
0.S14e19E+01
0.459674E+01
See
ene
,
Kel,
ZONETYPE=Ordered,
DATAPACKING=POINT
0.8344388-05
0.688726E-05
0.6925418-05
0.176291E-05
0.444817E-05
0.917042E-05
0.181247E-04
0.118887E-04
0.1232078-04
0.858124E-05
0.708475E-05
0.5382708-05
0.5045538-05
0.3218938-05
0.2924238-05
0.2225768-05
0.123322E-05
0.812749E-06
0.479991E-06
-0.856960E-07
-0.680339E-07
0.124498E-06
-0.508859E-07
0.5679918-06
0.345818E-06
0.1321938-05
0.3077138-05
0.3189288-05
0.3820348-05
0.1734582-05
0.5942938-05
0.671246E-05
0.7597908-05
0.7583508-05
0.8829788-05
0.881747E-05
0.9238968-05
0.9428428-05
0.319082E-05
0.692495E-05
0.987784-05
0.108076E-04
0.8027498-05
0.173376E-04
0.1008838-04
0.1017038-04
0.7426638-05
0.6512428-05
0.6836032-05
0.174441E-03
0.109562E-03
0.148294E-03
-0.185707E-03,
0.187052E-03
0.417461E-03
0.606706E-03
0.513817E-04
0.481839E-03
0.666172E-04
0.176420E-03
0.127530E-03
0.139813E-03
0.115866E-03
0.109487E-03
0.906497E-04
0.766180E-04
0.466642E-04
0.471922E-04
0.356322E-04
0..613023E-05
0.729104E-05
-0.108678E-04
0.186053E-05
-0.541985E-04
~0.684801E-04
-0.179995E-03
~-0.293061E-04
-0.996724E-04
0.120095E-03,
-0.376337E-03,
-0.857828E-04
-0.208736E-03
~0.114266E-03
-0.138731E-03
-0.170138E-03
-0.124290E-03,
-0.197509E-03
0.173158E-03
-0.265301E-03,
-0.212357E-03
-0.201051E-03
-0.366476E-04
-0.877737E-03,
0.142909E-03,
-0.222445E-03,
-0.171063E-03,
~0.132522E-03
~0.189092E-03
-0.2785S7E-04
-0.186209E-04
~0.213175E-04
0.145527E-04
-0.226298E-04
~0.474473E-04
-0.678707E-04
~0.230772E-04
~0.501874E-04
~0.190789E-04
~0.220023E-04
~0.161626E-04
~0.141329E-04
-9.930756E-05
-0.726243E-05
-9.518470E-05
~0.210627E-05
~0.127327E-06
0.593175E-06
0.222086E-05
0.293591E-05
0.234673E-05
0.207141E-05
0.20818SE-05
~0.389204E-06
~0.270759E-05
~0.125721E-04
-0.341816E-05
~0.111633E-04
0.999963E-05
~0.377933E-04
~0.146176E-04
-0.283220E-04
~0.205349E-04
~0.265344E-04
-0.294831E-04
-0.258977E-04
~0.339795E-04
0.100975E-04
-0.350964E-04
-0.36578SE-04
~0.318939E-04
~0.182452E-04
~0.955971E-04
0.
606289E-05
~0.339280E-04
~0.231766E-04
~0.197798E-04
~0.218734E-04
0.402298E+00
0.528743E+00
0.52072SE+00
0.516658E+00
0.508370E+00
0.4981¢8E+00
0.477363E+00
0.455536E+00
0.429878E+00
0.400548E+00
0.3678608+00
0.3322538+00
0.293924E+00
0.2531886+00
0.210431E+00
o.165941E+00
0.120059E+00
0.7318228-01
0.25667SE-01
~0.221587E-01
~0.698752E-01,
~0.117100E+00
-0.163506E+00
~0.208780E+00
~0.252343E+00
-0.294163E+00
-0.333871E+00
-0.371170E+00
-0.405741E+00
-0.437665E+00
~0.466089E+00
~0.491396E+00
~0.513230E+00
-9.531626E+00
~0.546173E+00
~0.886995SE+00
~0.564038E+00
~0.567169E+00
-0.566509E+00
-0.561938E+00
-0.$534:6E+00
-0.841347E+00
-0.525470E+00
-0.505938E+00
~0.482787E+00
~0.456482E+00
-0.426764E+00
~0.394348E+00
-9.388665E+00
85
Appendix B: One of the CFD output data
86
Appendix C: MATLAB codes
87
MATLAB code to calculate FFT spectrums and spectrum values clc;clear
all;close all
%%%processing data
%%%symbol('##' shows that we need to change this parameters aspect of each
%%%data
%%%last update 03Mar2014
%##########################################################################
%DATA INPUT
%##########################################################################
clc;clear;close all load '12_Forces_System_total
0.10_ 5.0.DAT' %##load data =
X12_Forces_System_total_0_10___5_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.1; %strouhal number %##St Q0
= 5.0; %pitch motion amplitude %##Q0 hc =
1; %heave to chord ratio %##h/c
h0 = hc*c; %heave motion amlitude
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_1__5_0_mot_freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.10_ 8.5.DAT' %##load
data = X12_Forces_System_total_0_10___8_5; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.1; %strouhal number %##St
Q0 = 8.5; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_1__8_5_mot_freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
88
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.10_ 12.0.DAT' %##load data
= X12_Forces_System_total_0_10__12_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.1; %strouhal number %##St
Q0 = 12.0; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_1__12_0_mot_freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n =
length(data(:,2)); %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.10_ 15.5.DAT' %##load data
= X12_Forces_System_total_0_10__15_5; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.1; %strouhal number %##St
Q0 = 15.5; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_1__15_5_mot_freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.10_ 19.0.DAT' %##load data
= X12_Forces_System_total_0_10__19_0; %##data
89
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.1; %strouhal number %##St
Q0 = 19.0; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_1__19_0_mot_freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.10_ 22.4.DAT' %##load
data = X12_Forces_System_total_0_10__22_4; %##data
U = 2.3; %flow velocity [m/s]; c = 1;
%chord length [m]
St = 0.1; %strouhal number %##St
Q0 = 22.4; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_1__22_4_mot_freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
clc;clear all;close all
%%%processing data
%%%symbol('##' shows that we need to change this parameters aspect of each
%%%data
%%%last update 04mar2014
%##########################################################################
%DATA INPUT
%##########################################################################
load '12_Forces_System_total 0.22_ 5.0.DAT' %##load
data = X12_Forces_System_total_0_22___5_0; %##data
90
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.22; %strouhal number %##St
Q0 = 5.0; %pitch motion amplitude %##Q0 hc
= 1; %heave to chord ratio %##h/c
h0 = hc*c; %heave motion amlitude f
= (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_22__5_0_motion freq [1/s]');disp(f)
T = 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode pause;
load '12_Forces_System_total 0.22_ 10.8.DAT' %##load data
= X12_Forces_System_total_0_22__10_8; %##data
U = 2.3; %flow velocity [m/s];
c = 1; %chord length [m]
St = 0.22; %strouhal number %##St
Q0 = 10.8; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_22__10_8 _motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.22_ 16.6.DAT' %##load data
= X12_Forces_System_total_0_22__16_6; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.22; %strouhal number %##St
Q0 = 16.6; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_22__16_6_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
91
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode pause;
load '12_Forces_System_total 0.22_ 22.3.DAT' %##load data
= X12_Forces_System_total_0_22__22_3; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.22; %strouhal number %##St
Q0 = 22.3; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_22__23_3_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.22_ 28.1.DAT' %##load data
= X12_Forces_System_total_0_22__28_1; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.22; %strouhal number %##St
Q0 = 28.1; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_22__28_1_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.22_ 33.9.DAT' %##load data
= X12_Forces_System_total_0_22__33_9; %##data
92
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.22; %strouhal number %##St
Q0 = 33.9; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_22__33_9_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.22_ 39.7.DAT' %##load data
= X12_Forces_System_total_0_22__39_7; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.22; %strouhal number %##St
Q0 = 39.7; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_22__39_7_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
clc;clear all;close all
%%%processing data
%%%symbol('##' shows that we need to change this parameters aspect of each
%%%data
%%%last update 01mar2014
%##########################################################################
%DATA INPUT
%##########################################################################
load '12_Forces_System_total 0.34_ 5.0.DAT' %##load
data = X12_Forces_System_total_0_34___5_0; %##data
93
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.34; %strouhal number %##St
Q0 = 5.0; %pitch motion amplitude %##Q0 hc
= 1; %heave to chord ratio %##h/c
h0 = hc*c; %heave motion amlitude
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_34__5_0_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.34_ 12.8.DAT' %##load
data = X12_Forces_System_total_0_34__12_8; %##data
U = 2.3; %flow velocity [m/s];
c = 1; %chord length [m]
St = 0.34; %strouhal number %##St
Q0 = 12.8; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_34__12_8_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.34_ 20.6.DAT' %##load data
= X12_Forces_System_total_0_34__20_6; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.34; %strouhal number %##St
Q0 = 20.6; %pitch motion amplitude %##Q0
94
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_34__20_6_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.34_ 28.4.DAT' %##load data
= X12_Forces_System_total_0_34__28_4; %##data U = 2.3;
%flow velocity [m/s]; c = 1; %chord length [m]
St = 0.34; %strouhal number %##St
Q0 = 28.4; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_34__28_4_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.34_ 36.3.DAT' %##load data
= X12_Forces_System_total_0_34__36_3; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.34; %strouhal number %##St
Q0 = 36.3; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_34__36_3_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load
95
'12_Forc
es_Syste
m_total
0.34_
44.1.DAT
'
%##load
data =
X12_Forc
es_Syste
m_total_
0_34__44
_1;
%##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.34; %strouhal number %##St
Q0 = 44.1; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_34__44_1_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.34_ 51.9.DAT' %##load data
= X12_Forces_System_total_0_34__51_9; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.34; %strouhal number %##St
Q0 = 51.9; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_34__51_9_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
96
basecode
clc;clea
r
all;clos
e all
%%%processing data
%%%symbol('##' shows that we need to change this parameters aspect of each
%%%data
%%%last update 01mar2014
%##########################################################################
%DATA INPUT
%##########################################################################
load '12_Forces_System_total 0.40_ 5.0.DAT' %##load
data = X12_Forces_System_total_0_40___5_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.40; %strouhal number %##St
Q0 = 5.0; %pitch motion amplitude %##Q0 hc
= 1; %heave to chord ratio %##h/c
h0 = hc*c; %heave motion amlitude f
= (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_40__5_0_motion freq [1/s]');disp(f)
T = 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.40_ 30.7.DAT' %##load
data = X12_Forces_System_total_0_40__30_7; %##data
U = 2.3; %flow velocity [m/s];
c = 1; %chord length [m]
St = 0.40; %strouhal number %##St
Q0 = 30.7; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_40__30_7_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
97
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.40_ 56.5.DAT' %##load
data = X12_Forces_System_total_0_40__56_5; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.40; %strouhal number %##St
Q0 = 56.5; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_40__56_5 _motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
clc;clear all;close
all
%%%processing data
%%%symbol('##' shows that we need to change this parameters aspect of each
%%%data
%%%last update 04mar2014
%##########################################################################
%DATA INPUT
%##########################################################################
load '12_Forces_System_total 0.46_ 5.0.DAT' %##load data =
X12_Forces_System_total_0_46___5_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
98
St = 0.46; %strouhal number %##St
Q0 = 5.0; %pitch motion amplitude %##Q0 hc
= 1; %heave to chord ratio %##h/c
h0 = hc*c; %heave motion amlitude
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_46__5_0_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.46_ 12.9.DAT' %##load data
= X12_Forces_System_total_0_46__12_9; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.46; %strouhal number %##St
Q0 = 12.9; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_46__12_9_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.46_ 20.8.DAT' %##load data
= X12_Forces_System_total_0_46__20_8; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.46; %strouhal number %##St
Q0 = 20.8; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_46__20_8_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
99
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.46_ 28.7.DAT' %##load data
= X12_Forces_System_total_0_46__28_7; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.46; %strouhal number %##St
Q0 = 28.7; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_46__28_7_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.46_ 36.6.DAT' %##load data
= X12_Forces_System_total_0_46__36_6; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.46; %strouhal number %##St
Q0 = 36.6; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_46__36_6_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.46_ 44.5.DAT' %##load data
= X12_Forces_System_total_0_46__44_5; %##data
100
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.46; %strouhal number %##St
Q0 = 44.5; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_46__44_5_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode pause;
load '12_Forces_System_total 0.46_ 52.4.DAT' %##load data
= X12_Forces_System_total_0_46__52_4; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.46; %strouhal number %##St
Q0 = 52.4; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_46__52_4_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.46_ 60.3.DAT' %##load data
= X12_Forces_System_total_0_46__60_3; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.46; %strouhal number %##St
Q0 = 60.3; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_46__60_3_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
101
psi = 90; %phase angle
t = data(:,2); %[s] cfd running time
n = length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
clc;clear all;close all
%%%processing data
%%%symbol('##' shows that we need to change this parameters aspect of each
%%%data
%%%last update 04mar2014
%##########################################################################
%DATA INPUT
%##########################################################################
load '12_Forces_System_total 0.58_ 5.0.DAT' %##load
data = X12_Forces_System_total_0_58___5_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.58; %strouhal number %##St
Q0 = 5.0; %pitch motion amplitude %##Q0 hc
= 1; %heave to chord ratio %##h/c
h0 = hc*c; %heave motion amlitude
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_58__5_0_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.58_ 13.7.DAT' %##load
data = X12_Forces_System_total_0_58__13_7; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.58; %strouhal number %##St
Q0 = 13.7; %pitch motion amplitude %##Q0
102
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_58__13_7_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.58_ 22.5.DAT' %##load data
= X12_Forces_System_total_0_58__22_5; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.58; %strouhal number %##St
Q0 = 22.5; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_58__22_5_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.58_ 31.2.DAT' %##load data
= X12_Forces_System_total_0_58__31_2; %##data U = 2.3;
%flow velocity [m/s]; c = 1; %chord length [m]
St = 0.58; %strouhal number %##St
Q0 = 31.2; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_58__31_2_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
103
basecode
pause;
load '12_Forces_System_total 0.58_ 40.0.DAT' %##load data
= X12_Forces_System_total_0_58__40_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.58; %strouhal number %##St
Q0 = 40.0; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_58__40_0_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load
'12_Forc
es_Syste
m_total
0.58_
48.7.DAT
'
%##load
data =
X12_Forc
es_Syste
m_total_
0_58__48
_7;
%##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.58; %strouhal number %##St
Q0 = 48.7; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_58__48_7_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
104
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.58_ 57.5.DAT' %##load data
= X12_Forces_System_total_0_58__57_5; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.58; %strouhal number %##St
Q0 = 57.5; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_58__57_5motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.58_ 66.2.DAT' %##load
data = X12_Forces_System_total_0_58__66_2; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.58; %strouhal number %##St
Q0 = 66.2; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_58__66_2_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
clc;clear all;close
all
%%%processing data
105
%%%symbol('##' shows that we need to change this parameters aspect of each
%%%data
%%%last update 01mar2014
%##########################################################################
%DATA INPUT
%##########################################################################
load '12_Forces_System_total 0.60_ 5.0.DAT' %##load
data = X12_Forces_System_total_0_60___5_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.60; %strouhal number %##St
Q0 = 5.0; %pitch motion amplitude %##Q0 hc
= 1; %heave to chord ratio %##h/c
h0 = hc*c; %heave motion amlitude
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_60__5_0_motion freq
[1/s]');disp(f) T = 1/f;
%motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.60_ 36.0.DAT' %##load data
= X12_Forces_System_total_0_60__36_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.60; %strouhal number %##St
Q0 = 36.0; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_60__36_0_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
106
basecode
pause;
load '12_Forces_System_total 0.60_ 67.1.DAT' %##load data
= X12_Forces_System_total_0_60__67_1; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.60; %strouhal number %##St
Q0 = 67.1; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_60__67_1_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
clc;clear all;close
all
%%%processing data
%%%symbol('##' shows that we need to change this parameters aspect of each
%%%data
%%%last update 01mar2014
%##########################################################################
%DATA INPUT
%##########################################################################
load '12_Forces_System_total 0.70_ 5.0.DAT' %##load
data = X12_Forces_System_total_0_70___5_0; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.70; %strouhal number %##St
Q0 = 5.0; %pitch motion amplitude %##Q0 hc
= 1; %heave to chord ratio %##h/c
h0 = hc*c; %heave motion amlitude
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_70__5_0motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
107
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.70_ 14.4.DAT' %##load data
= X12_Forces_System_total_0_70__14_4; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.70; %strouhal number %##St
Q0 = 14.4; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_70__14_4_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.70_ 23.7.DAT' %##load data
= X12_Forces_System_total_0_70__23_7; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.70; %strouhal number %##St
Q0 = 23.7; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_70__23_7_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode pause;
load '12_Forces_System_total 0.70_ 33.1.DAT' %##load data
= X12_Forces_System_total_0_70__33_1; %##data
108
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.70; %strouhal number %##St
Q0 = 33.1; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_70__33_1_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.70_ 42.5.DAT' %##load data
= X12_Forces_System_total_0_70__42_5; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.70; %strouhal number %##St
Q0 = 42.5; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_70__42_5_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running timee
basecode pause;
load '12_Forces_System_total 0.70_ 51.8.DAT' %##load data
= X12_Forces_System_total_0_70__51_8; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.70; %strouhal number %##St
Q0 = 51.8; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_70__51_8_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
109
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.70_ 61.2.DAT' %##load data
= X12_Forces_System_total_0_70__61_2; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.70; %strouhal number %##St
Q0 = 61.2; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_70__61_2_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle
t = data(:,2); %[s] cfd running time
n = length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
pause;
load '12_Forces_System_total 0.70_ 70.5.DAT' %##load
data = X12_Forces_System_total_0_70__70_5; %##data
U = 2.3; %flow velocity [m/s]; c
= 1; %chord length [m]
St = 0.70; %strouhal number %##St
Q0 = 70.5; %pitch motion amplitude %##Q0
f = (U*St)/(2*h0); %motion freq.[1/s]
fprintf('0_70__70_5_motion freq [1/s]');disp(f) T
= 1/f; %motion period [s]
Qm = 0; %mean angle of attact [deg]
psi = 90; %phase angle t =
data(:,2); %[s] cfd running time n
= length(data(:,2)) %time vector length
tim_run = t(n,1); %cfd running time
basecode
110
Basecode
%##########################################################################
%DATA OUTPUT
%##########################################################################
% omega shows motion of freq. [rd/s]
omega = 2*pi*f; %the motions of
wings %heave motion ht =
h0*sin(2*pi*f*t);
%pitch motion
Qt = Qm + Q0*sin(2*pi*f*t+psi);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plotting heave force, fy
% h=figure('visible','on');
% plot(data(:,2),data(:,4)),grid on
% title('heave force')
% xlabel('time [s]')
% ylabel('fy')
% finding magnitude of motion_freq in spectrum fy
[minVal, position] = min(abs(t-T)); hf
= ht(position,1);
%fft(fast foruier transform) for fy signal
h=figure('visible','on'); dt = tim_run/n;
df = 1/((n-1)*dt);
DSideSpectrumfy = fft(data(:,4));
Spectrumfy = DSideSpectrumfy(1:n)/n;
freq_axis = df*(0:1:floor(n/2-1));
Spectrumhalffy = abs(Spectrumfy(1:floor(n/2)));
semilogy(freq_axis,Spectrumhalffy),grid on;
title('fft heave force') xlabel('frequence
[1/s]') ylabel('Spectrum')
% finding x1magnitude of motion_freq in spectrum fy
[minVal,position]=min(abs(freq_axis - f)); fyf =
Spectrumfy(position,1) magfyf = abs(fyf) angradfyf
= angle(fyf); angfyf = radtodeg(angradfyf)
% %finding x3magnitude of motion_freq in spectrum fy
% [minVal,position]=min(abs(freq_axis - 3*f));
% fy3f = Spectrumfy(position,1); %
magfy3f = abs(fy3f);
% angradfy3f = angle(fy3f);
% angfy3f = radtodeg(angradfy3f);
% %finding x5magnitude of motion_freq in spectrum fy
% [minVal,position]=min(abs(freq_axis - 5*f));
% fy5f = Spectrumfy(position,1);
% magfy5f = abs(fy5f);
% angradfy5f = angle(fy5f);
% angfy5f = radtodeg(angradfy5f);
%fft angle for heave force
h=figure('visible','off'); angrad
= angle(DSideSpectrumfy); angdeg =
radtodeg(angrad); angdeghalf =
111
angdeg(1:floor(n/2));
plot(freq_axis,angdeghalf),grid on;
title('fft heave force angle')
xlabel('frequence [1/s]')
ylabel('degree')
MATLAB code to calculate FFT spectrums [Developed by Dr. Nikolas Xiros] clc;
clear all; close all;
FishTable = {'fish 15\','fish 30\'};
ARTable = {'AR4\','AR6\'};
HCTable = {'HC 1\','HC 1_5\','HC 2\'};
FilenameTable = {'12_Forces_System_total 0.10_ 5.0.DAT',...
'12_Forces_System_total 0.10_ 8.5.DAT','12_Forces_System_total 0.10_
12.0.DAT',...
'12_Forces_System_total 0.10_ 15.5.DAT','12_Forces_System_total 0.10_
19.0.DAT',...
'12_Forces_System_total 0.10_ 22.4.DAT','12_Forces_System_total 0.22_
5.0.DAT',...
'12_Forces_System_total 0.22_ 10.8.DAT','12_Forces_System_total 0.22_
16.6.DAT',...
'12_Forces_System_total 0.22_ 22.3.DAT','12_Forces_System_total 0.22_
28.1.DAT',...
'12_Forces_System_total 0.22_ 33.9.DAT','12_Forces_System_total 0.22_
39.7.DAT',...
'12_Forces_System_total 0.34_ 5.0.DAT','12_Forces_System_total 0.34_
12.8.DAT',...
'12_Forces_System_total 0.34_ 20.6.DAT','12_Forces_System_total 0.34_
28.4.DAT',...
'12_Forces_System_total 0.34_ 36.3.DAT','12_Forces_System_total 0.34_
44.1.DAT',...
'12_Forces_System_total 0.34_ 51.9.DAT','12_Forces_System_total 0.40_
5.0.DAT',...
'12_Forces_System_total 0.40_ 30.7.DAT','12_Forces_System_total 0.40_
56.5.DAT',...
'12_Forces_System_total 0.46_ 5.0.DAT','12_Forces_System_total 0.46_
12.9.DAT',...
'12_Forces_System_total 0.46_ 20.8.DAT','12_Forces_System_total 0.46_
28.7.DAT',...
'12_Forces_System_total 0.46_ 36.6.DAT','12_Forces_System_total 0.46_
44.5.DAT',...
'12_Forces_System_total 0.46_ 52.4.DAT','12_Forces_System_total 0.46_
60.3.DAT',...
'12_Forces_System_total 0.58_ 5.0.DAT','12_Forces_System_total 0.58_
13.7.DAT',...
'12_Forces_System_total 0.58_ 22.5.DAT','12_Forces_System_total 0.58_
31.2.DAT',...
'12_Forces_System_total 0.58_ 40.0.DAT','12_Forces_System_total 0.58_
48.7.DAT',...
'12_Forces_System_total 0.58_ 57.5.DAT','12_Forces_System_total 0.58_
66.2.DAT',...
'12_Forces_System_total 0.60_ 5.0.DAT','12_Forces_System_total 0.60_
36.0.DAT',...
'12_Forces_System_total 0.60_ 67.1.DAT','12_Forces_System_total 0.70_
5.0.DAT',...
112
'12_Forces_System_total 0.70_ 14.4.DAT','12_Forces_System_total 0.70_
23.7.DAT',...
'12_Forces_System_total 0.70_ 33.1.DAT','12_Forces_System_total 0.70_
42.5.DAT',...
'12_Forces_System_total 0.70_ 51.8.DAT','12_Forces_System_total 0.70_
61.2.DAT',...
'12_Forces_System_total 0.70_ 70.5.DAT'};
StrouhalTable = [0.10,...
0.10,0.10,...
0.10,0.10,... 0.10,0.22,...
0.22,0.22,...
0.22,0.22,...
0.22,0.22,...
0.34,0.34,... 0.34,0.34,...
0.34,0.34,... 0.34,0.40,...
0.40,0.40,... 0.46,0.46,...
0.46,0.46,... 0.46,0.46,...
0.46,0.46,... 0.58,0.58,...
0.58,0.58,... 0.58,0.58,...
0.58,0.58,... 0.60,0.60,...
0.60,0.70,... 0.70,0.70,...
0.70,0.70,...
0.70,0.70,...
0.70];
Q0Table = [5.0,...
8.5, 12.0,...
15.5,19.0,...
22.4,5.0,...
10.8,16.6,...
22.3,28.1,...
33.9,39.7,...
5.0,12.8,...
20.6,28.4,...
36.3,44.1,...
51.9,5.0,...
30.7,56.5,...
5.0,12.9,...
20.8,28.7,...
36.6,44.5,...
52.4,60.3,...
5.0,13.7,...
22.5,31.2,...
40.0,48.7,...
57.5,66.2,...
5.0,36.0,...
67.1,5.0,...
14.4,23.7,...
33.1,42.5,...
51.8,61.2,...
70.5];
U = 2.3; c = 1;
Qm = 0; %mean angle of attact [deg] psi
= 90; %phase angle
113
Str0 = 'C:\ftp NTUA data_procnew\'; % ATTENTION !! Define this before
running the program !!! Str00 = strcat(Str0,'fishFDP\');
Indx =
0;
for i1=1:length(FishTable) for i2 =
1:length(ARTable) for i3 =
1:length(HCTable) hc = .5 +
i3*.5; h0 = hc*c;
for i4 = 1:length(FilenameTable)
Indx = Indx+1;
Strouhal = StrouhalTable(i4);
Q0 = Q0Table(i4); f = (U*Strouhal)/(2*h0);
%motion freq.[1/s] %T = 1/f; %motion
period [s] ts1 =
strcat(FishTable{i1},ARTable{i2},HCTable{i3},FilenameTable{i4});
IndexVector(Indx).filename = ts1;
ts2 = strcat(Str00,'fdp_',num2str(Indx),'.mat');
IndexVector(Indx).dest = ts2;
PathStr = strcat(Str0,ts1);
DataSeries = load(PathStr);
%frequency analysis code here
if length(DataSeries) > 0
n = DataSeries(length(DataSeries),1);
tim = DataSeries(:,2); dt = tim(2)-
tim(1); df = 1/((n-1)*dt);
fx = DataSeries(:,3); fy =
DataSeries(:,4); fz =
DataSeries(:,5); mx =
DataSeries(:,6); my =
DataSeries(:,7); mz =
DataSeries(:,8); freqAxis =
df*(0:1:floor(n/2-1)); n2Freq =
length(freqAxis); FFT_ = fft(fx)/n;
fxFFT = FFT_(1:n2Freq); FFT_ =
fft(fy)/n; fyFFT = FFT_(1:n2Freq);
FFT_ = fft(fz)/n; fzFFT =
FFT_(1:n2Freq); FFT_ = fft(mx)/n;
mxFFT = FFT_(1:n2Freq); FFT_ =
fft(my)/n; myFFT = FFT_(1:n2Freq);
FFT_ = fft(mz)/n; mzFFT =
FFT_(1:n2Freq);
FFT_Tabl = [freqAxis', fxFFT, fyFFT, fzFFT, mxFFT, myFFT,
mzFFT]; else df = 0;
n2Freq = 0; FFT_Tabl = []; end;
DataTable.codename = PathStr;
DataTable.features = [Indx, length(DataSeries), Strouhal, h0, Q0,
f, df, n2Freq];
IndexVector(Indx).feats = [Indx, length(DataSeries),
Strouhal, h0, Q0, f, df, n2Freq];
DataTable.data = DataSeries;
DataTable.FFT = FFT_Tabl;
save(ts2,'DataTable'); end; end;
114
end; end; save
(strcat(Str00,'IndexVector.mat'),'IndexVector');
Indx
% DO NOT FORGET TO SAVE DataTable !!!!
%%-- END nx
MATLAB code to generate FFT Spectrum plots [Developed by Dr. Nikolas Xiros] clear
all; close all;
PLOTfx = figure;
PLOTfy = figure;
PLOTfz = figure;
PLOTmx = figure;
PLOTmy = figure;
PLOTmz = figure; load 'C:\ftp NTUA
data_procnew\fishFDP\IndexVector.mat';
for i =
1:length(IndexVector)
kati = IndexVector(i).feats; titlos = strcat('Run',
num2str(kati(1)),': Strouhal=', num2str(kati(3)),' h0=', num2str(kati(4)),'
Q0=', num2str(kati(5)),' F1=', num2str(kati(6)));
if
(kati(2)>0)
load(IndexVector(i).dest);
fft_ = DataTable.FFT;
figure(PLOTfx);
semilogy(fft_(:,1),abs(fft_(:,2)),'k-','Linewidth',2);
grid; xlabel('frequency F in Hz');
ylabel('|fx(F)|'); title(titlos); testr =
strcat('C:\ftp NTUA
data_procnew\FishSpecPlots\fx\',num2str(kati(1)),'_fx.jpg');
print(PLOTfx,'-djpeg',testr);
figure(PLOTfy);
semilogy(fft_(:,1),abs(fft_(:,3)),'b-','Linewidth',2); grid;
xlabel('frequency F in Hz');
ylabel('|fy(F)|');
title(titlos); testr =
strcat('C:\ftp NTUA
data_procnew\FishSpecPlots\fy\',num2str(kati(1)),'_fy.jpg');
print(PLOTfy,'-djpeg',testr);
figure(PLOTfz);
semilogy(fft_(:,1),abs(fft_(:,4)),'r-','Linewidth',2);
grid; xlabel('frequency F in Hz');
ylabel('|fz(F)|'); title(titlos); testr =
strcat('C:\ftp NTUA
data_procnew\FishSpecPlots\fz\',num2str(kati(1)),'_fz.jpg');
print(PLOTfz,'-djpeg',testr);
115
figure(PLOTmx);
semilogy(fft_(:,1),abs(fft_(:,5)),'k--','Linewidth',3);
grid; xlabel('frequency F in Hz');
ylabel('|mx(F)|'); title(titlos); testr =
strcat('C:\ftp NTUA
data_procnew\FishSpecPlots\mx\',num2str(kati(1)),'_mx.jpg');
print(PLOTmx,'-djpeg',testr);
figure(PLOTmy);
semilogy(fft_(:,1),abs(fft_(:,6)),'b--','Linewidth',3);
grid; xlabel('frequency F in Hz');
ylabel('|my(F)|'); title(titlos); testr =
strcat('C:\ftp NTUA
data_procnew\FishSpecPlots\my\',num2str(kati(1)),'_my.jpg');
print(PLOTmy,'-djpeg',testr);
figure(PLOTmz);
semilogy(fft_(:,1),abs(fft_(:,7)),'r--','Linewidth',3);
grid; xlabel('frequency F in Hz');
ylabel('|mz(F)|'); title(titlos); testr =
strcat('C:\ftp NTUA
data_procnew\FishSpecPlots\mz\',num2str(kati(1)),'_mz.jpg');
print(PLOTmx,'-djpeg',testr); end; end;
%close all;
MATLAB codes to find Transfer function model
The Code to calculate selected operational points for each case
%%Table for fish 15\AR4\HC1 for Q0=5 %
Table generating for generating H = 1
clc;clear all
%INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('data1.mat')
% load '12_Forces_System_total 0.10_ 5.0.DAT' %##load
% data = X12_Forces_System_total_0_10___5_0; %##data
%%%%CHANGE compnum
= data1(:,4).*exp(data1(:,5).*(pi/180)*sqrt(-1));
%%%%CHANGE
Table_fish_15_AR4_HC1 = [data1(:,1), data1(:,2),data1(:,3),compnum]; %%%%CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
116
squence_pitmot_H1_1 = sortrows(Table_fish_15_AR4_HC1,3);
%%% we have chosen to inversitgate fish_15_AR4 type of fin
Table_fish_15_AR4_pit5 = [squence_pitmot_H1_1(1,:);...
;squence_pitmot_H1_1(2,:)...
;squence_pitmot_H1_1(3,:)...
;squence_pitmot_H1_1(5,:)...
;squence_pitmot_H1_1(6,:)...
;squence_pitmot_H1_1(8,:)]
omega_5deg = (2*pi)*Table_fish_15_AR4_pit5(:,1); %2*pi*F radyal
freq.
datapoints_5deg = Table_fish_15_AR4_pit5(:,4); Table_omeg_datapoints
= [omega_5deg,datapoints_5deg];
%%%%%degistirmene gerek yok birseyi
squ_Table_omeg_datapoints = sortrows(Table_omeg_datapoints,1) %make squence with
respec to freqs.
%CFD result h=figure('visible','on');
y = 20*log10(abs(squ_Table_omeg_datapoints(:,2))); x
= squ_Table_omeg_datapoints(:,1);
semilogx(x,y,'r'),hold on,grid on
title('FinType=fish15AR4, PitchAmp=5Deg, HeaveAmp=1m')
xlabel('omega [rad/s]') ylabel('|CFD H(F)|')
save('squ_Table_omeg_datapoints'); omega
= squ_Table_omeg_datapoints(:,1) save
omega
datapoints = squ_Table_omeg_datapoints(:,2) save
datapoints
%saveas(h,[pwd '/fish15_AR4_pit5_heav_1.jpg'])
The Code to provide approx. Transfer function function
Hvalue = H_hat(u,omegvalue)
k =
u(1); z
= u(2);
s = j*omegvalue;
Hhat = k*(s+z)/((s+100)*(s+200));
Hvalue = Hhat;
end
The Code to provide the objective function
%%defining objective function
117
function J = Objfunc(u)
load('datapoints.mat') load('omega.mat')
Jaux=0; for count =
1:1:length(omega);
Jaux = Jaux + (20*(log10(abs(H_hat(u,omega(count)))) -
log10(abs(datapoints(count)))))^2;
end J =
Jaux;
The Code to provide Optimization
% Optimization clc;clear
all;close all
load('datapoints.mat') load('omega.mat')
u0 = [1000 .0001]; %random intial variables opts =
optimset('algorithm','interior-point','tolfun',10^-14000,'tolx'...
,10^-14000,'MaxFunEvals', 10^5,'MaxIter',10^5);
[u,fval] = fminunc(@Objfunc,u0, opts);
%design variables
k = u(1); z =
u(2); %optimum
results OptResuts
= [k z]; k =
OptResuts(:,1);
display(k)
z = OptResuts(:,2); display(z)
%Obj fun value with initial results
ObjFunc_InVal = Objfunc(u0);
display(ObjFunc_InVal)
%Obj fun value with optimum results
ObjFunc_OptVal = Objfunc(u);
display(ObjFunc_OptVal)
%plotting
H_aproxx = [H_hat(u,omega(1))... %*
;H_hat(u,omega(2))...
;H_hat(u,omega(3))...
;H_hat(u,omega(4))...
;H_hat(u,omega(5))...
;H_hat(u,omega(6))];
%CFD result
h=figure('visible','on'); %
h=figure('Units', 'pixels', ...
% 'Position', [100 100 500 375]);
%hold on; %if uncoomend real value is coming
%CFD data
118
y1 = 20*log10(abs(datapoints));
semilogx(omega,y1,'b-+','Linewidth',1.2);
hold on; grid on;
%Transfer function modal y2 =
20*log10(abs(H_aproxx));
semilogx(omega,y2,'k--*','Linewidth',2);
hold on; grid on;
title('CFD vs Tr. Func. Modal','FontWeight','Bold','FontSize',12)
xlabel('log \omega [rad/s]','FontWeight','Bold','FontSize',12) ylabel('log
|H(F)|','FontWeight','Bold','FontSize',12)
legend('CFD data','Tr. Func. Modal','Location','Best'); %text
(.77,32.5,'FinType:fish15AR4 (5deg,1m)') ylim=get(gca,'ylim');
xlim=get(gca,'xlim'); text(xlim(1)+0.03,ylim(2)-2.5,'FinType:fish15AR4
(5deg,1m)') %*
saveas(h,[pwd
'/comp_fish15_AR4_pit5_heav_1.jpg'])
MATLAB code to the points of transfer function values of the system and generate 3-D
plots
clc; clear all; close all;
load
'SurgeTable.mat'
X = surgeTable; %table of the surge force
%F[hertz]1 | h0(heave amp)2 | q0(pitch amp)deg3 | f(0)4 | posf(0)5 | f(2F)6 |
posf(2F) |
h0 = X(:,2); %heave amp
q0deg = X(:,3); %pitch amp deg
q0 = X(:,3)*(pi/180); %pitch amp rad
A = sqrt(h0.^2 + q0.^2); %virtual fictitious amp - combined amp
b = X(:,6); %this is coefficient coming from f(2F) in fourier domain
%%Now we can calculate our transfer function
%Calc of H(F,F) ((Fx(2F) f = 2F))
array1 = -(b.*8*pi*1i); array2 =
A.^2;
HFF = array1./array2; % check element-wise operations in matlab
%Calc of H(-F,-F) ((Fx(-2F) f = 2F))
%we can easily compte H(-F,-F) because of this property
%H(-F,-F) = - conj[H(F,F)] conf = Complex conjugate
H_F_F = - conj(HFF); %note that H_F_F means H(-F,-F)
%Calc of H(-F,F) ((Fx(0) f = 0))
%H(F,-F)+H(-F,F) = a %a is coefficient coming from f(0) in fourier domain
%H(F,-F) = -conj(H(-F,F))
a = X(:,4); array3 =
a.*8*pi*1i; array4 =
2*A.^2; Im_H_FF =
array3./array4;
%Calc of H(F,-F) ((Fx(0) f = 0))
%we can easily compte H(F,-F) because of this property
119
%H(F,-F) = - conj[H(-F,F)] conf = Complex conjugate
Im_HF_F = -conj(Im_H_FF);
%All Second order transfer functions is collected together
Volterra_Tr_Func = [HFF,H_F_F,Im_HF_F,Im_H_FF];
Volterra_Tr_Func_separate_real_and_imag = [real(HFF) imag(HFF),real(H_F_F)
imag(H_F_F),real(Im_HF_F) imag(Im_HF_F),real(Im_H_FF) imag(Im_H_FF)];
%Volterra_Tr_Func_2dec = floor(Volterra_Tr_Func*1e2) / 1e2; %2 decimals
%%Now we define Volterra_Tr_Func as magnitude and phase (another def of
%%complex number ) %for
H(F,F)
mag_HFF = abs(Volterra_Tr_Func(:,1)); %magnitude of HFF ph_HFF_radian
= angle(Volterra_Tr_Func(:,1)); %phase of HFF ph_HFF_deg =
ph_HFF_radian*(180/pi);
%for H(-F,-F)
mag_H_F_F = abs(Volterra_Tr_Func(:,2)); %magnitude of HFF
ph_H_F_F_radian = angle(Volterra_Tr_Func(:,2)); %phase of HFF ph_H_F_F_deg
= ph_H_F_F_radian*(180/pi);
%for H(-F,F) mag_H_FF = abs(Volterra_Tr_Func(:,3));
%magnitude of HFF ph_H_FF_radian =
angle(Volterra_Tr_Func(:,3)); %phase of HFF ph_H_FF_deg =
ph_H_FF_radian*(180/pi);
%for H(F,-F)
mag_HF_F = abs(Volterra_Tr_Func(:,4)); %magnitude of HFF ph_HF_F_radian
= angle(Volterra_Tr_Func(:,4)); %phase of HFF ph_HF_F_deg =
ph_HF_F_radian*(180/pi);
%% Making mag-ph table F H(F,F) H(-F,-F) H(-F,F) H(F,-F)
Volterra_Tr_Func_mag_ph = [X(:,1) mag_HFF ph_HFF_deg mag_H_F_F ph_H_F_F_deg
mag_H_FF ph_H_FF_deg mag_HF_F ph_HF_F_deg]; squ_Volterra_Tr_Func_mag_ph =
sortrows(Volterra_Tr_Func_mag_ph,1);
Y = squ_Volterra_Tr_Func_mag_ph;
F = squ_Volterra_Tr_Func_mag_ph(:,1);
%% 3d plot processing for H(F,F)
%magnitude
h=figure('visible','off'); squ_mag_HFF = Y(:,2);
plot3(F,F,squ_mag_HFF,'.'),grid on title('Volterra
Kernel Points') xlabel('F[1/s]') ylabel('F[1/s]')
zlabel('magnitude of H(F,F)') saveas(h,[pwd
'/3D_Tr_Func_plots/1_mag_H(F,F).jpg'])
%phase
h=figure('visible','off');
squ_ph_HFF = Y(:,3);
plot3(F,F,squ_ph_HFF,'.'),grid on
title('Volterra Kernel Points')
xlabel('F[1/s]') ylabel('F[1/s]')
zlabel('phase of H(F,F) [deg]') saveas(h,[pwd
'/3D_Tr_Func_plots/2_ph_H(F,F).jpg'])
%% 3d plot processing for H(-F,-F)
%magnitude
h=figure('visible','off'); squ_mag_H_F_F
= Y(:,4);
120
plot3(-F,-F,squ_mag_H_F_F,'.'),grid on
title('Volterra Kernel Points') xlabel('-F[1/s]')
ylabel('-F[1/s]') zlabel('magnitude of H(-F,-F)')
saveas(h,[pwd '/3D_Tr_Func_plots/3_mag_H(-F,-F).jpg'])
%phase
h=figure('visible','off'); squ_ph_H_F_F
= Y(:,5);
plot3(-F,-F,squ_ph_H_F_F,'.'),grid on
title('Volterra Kernel Points')
xlabel('-F[1/s]') ylabel('-F[1/s]')
zlabel('phase of H(-F,-F) [deg]')
saveas(h,[pwd '/3D_Tr_Func_plots/4_ph_H(-F,-F).jpg'])
%% 3d plot processing for H(-F,F)
%magnitude h=figure('visible','off'); squ_mag_H_FF
= Y(:,6); plot3(-F,F,squ_mag_H_FF,'.'),grid on
title('Volterra Kernel Points') xlabel('-F[1/s]')
ylabel('F[1/s]') zlabel('magnitude of H(-F,F)')
saveas(h,[pwd '/3D_Tr_Func_plots/5_mag_H(-F,F).jpg'])
%phase
h=figure('visible','off'); squ_ph_H_FF = Y(:,7);
plot3(-F,F,squ_ph_H_FF,'.'),grid on title('Volterra
Kernel Points') xlabel('-F[1/s]') ylabel('F[1/s]')
zlabel('phase of H(-F,F) [deg]') saveas(h,[pwd
'/3D_Tr_Func_plots/6_ph_H(-F,F).jpg'])
%% 3d plot processing for H(F,-F)
%magnitude h=figure('visible','off'); squ_mag_HF_F
= Y(:,8); plot3(F,-F,squ_mag_HF_F,'.'),grid on
title('Volterra Kernel Points') xlabel('F[1/s]')
ylabel('-F[1/s]') zlabel('magnitude of H(F,-F)')
saveas(h,[pwd '/3D_Tr_Func_plots/7_mag_H(F,-F).jpg'])
%phase
h=figure('visible','off');
squ_ph_HF_F = Y(:,9); plot3(F,-
F,squ_ph_HF_F,'.'),grid on
title('Volterra Kernel Points')
xlabel('F[1/s]') ylabel('-F[1/s]')
zlabel('phase of H(F,-F) [deg]')
saveas(h,[pwd '/3D_Tr_Func_plots/8_ph_H(F,-F).jpg'])
121
Appendix D: Spectrum plots for fish15AR4 and fish15AR6
Heave force
en
en
en
o
Ew
Ew
1"
1" 1"
ey
G
1
ey
s
0
ey
2 + 6 6
0
wz
1
omen
Fine
omen
Fine
renee
Finte
w
w
@
ge
2
ais
1"
Ww
ey
+66
0
wz
1
ey
2 + 6 6
0
wz
ee
2 + 6 6
0
wz
1
renee
Finte
renee
Finte
renee
Finte
w
o
E
Ew
1"
1"
Ww
ey
+66
0
wz
1
2 + 6 6
0
wz
ey
2
cy
6 6
0
wz
1
renee
Finte
renee
Finte
renee
Finte
o
Ew
Ew
1"
1" 1"
ey
+ 6
6
0
wz
1
ey
1
2
2
« s
6
7ea
ey
1
2
2
« s
6
7ea
1
renee
Finte
remeneFnt| remeneFnt|
w
w
E
Ew
1"
1"
w
ey
2
« s
6
7ea
1
1
2
2
« s
6
7ea
ey
1
2
2
« s
6
7ea
1
remeneFnt| remeneFnt|
renee
Fnt
w w
wh
W
= =
w
ww
wel
2
«s 7ea
0
1
2
2
«s 7ea
0
ee
1
2
2
« s
6
788
0
remeneFnt| remeneFnt| remeneFnt|
122
ie
ie
ie
aaa
.
78
+ + +
fromeney
Fine
fromeney
Fine
fremeney
Fine
ie
a
a
a
ie
ie
+
2 3 4 .
78
+e
Se
se
a
.
78
tremens
Fine
tremens
Fine
tremens
Fine
o
oo oo
ca
= = =
1
ra
10"
10" 10° 10°
o
+
234
sO
o
ost
4s
2
35
3
a5
4
as
5
o
ost
‘5
2
35
3
35
4
45
5
tremens
Fine
remnsy
Fine
remnsy
Fine
ir
4
a
a
ir
ir
Wo
st
sas
Wo
sts
as
Wo
st
sas
remency
Fine
remency
Fine
remency
Fine
ir
44
8
4&
ie
ir
ots
es
asa
ot
Ws
eS
°
oF
+
i
2
remency
Fine
remency
Fine
remensy
Fine
a
3 é
wi
a8
se
wrt
i
i
i
remensy
Fine
remensy
Fine
remeney
Fine
123
en
cea
en
0:
Sroind
7
Het
5006705
F1n0
9857
o o
0
2°
2
» »
«
a5
5
#3
3 2
5
#3
oa
roomy
Fine
roomy
Fine
ney
Fine
fats:
hei?
Scot
2
Fosse
.
ct
aa
7h
cost
FOS
.
Ica
ni
oO
F087
o
Ew
z
a a
»!
w=
o
2 +
6
e
10.
2“6
7%
~ ey
2 +
6
e
10.
2“6w
2 ee
2 +
6
e
10.
2“6w
2
acente acente acente
o
Ew Ew
——
onc
sata
1
»!
ais
:
ee ee
oa
oa
oa
acente acente acente
o
Ew Ew
—
,
ho
eeepers
:
A
»!
:
®
wo
oa
oi
oi
ney
Fine
sete
Fe
sete
Fe
o
z
Ew
a a
{
A
w
o
:
4
6 s 10
2
“
16
¥
s 10
2
“
ey
:
4
6 s 10
2
“
16
w
pbeeyrinte pbeeyrinte pbeeyrinte
Ew Ew
senses
.
con
ssasaae
4
wa
»!
2
wo wo
oi oi oi
pbeeyrinte pbeeyrinte pbeeyrinte
124
en
en
en
saa
@ *
fremensy
Fine
@ *
fremensy
Fine
[
fremensy
Fine
@ *
fremensy
Fine
Cn
remnsy
Fine
o
ca
cea
rca
cea
Cn
remnsy
Fine
5
remnsy
Fine
Cn
remnsy
Fine
wo
ca ca
= =
10" 10"
+
2
3
« 3
6
7@8
a
a
a
remency
Fine
remency
Fin
te
remency
Fin
te
co
§
gw
1
10"
7
23+5
‘
fremensy
Fin
te
3 ‘
fremensy
Fin
te
125
we
: E
‘0 ‘0
0}
wo we
a"
= 2
1.
ca ca
10"
+
73‘
> €
7
+
73‘
> €
7
we
+
73‘
>
seenaney
Fie
seenaney
Fie
seeps
Finite
oy
wo wo
wo
Eo Ew
Eo
1" 1"
ca
we
+
73‘
= €
7
oF
is
7
3
3
we
oF
+
is
seeps
Finite
sep
Fine
freer)
Fine
1
co
‘
= = =
7 7 7
we
oF
+
is
3
3
we
oF
+
is
7
3
3
we
oF
+
is
7
sep
Fine
sep
Fine
fees
Fine
.
eco
Seana
thoet
Som-5F-0Oe867
.
00
Seth
7
het
GOO
20
205
.
an
Sete
7
nt
GOW
250205
we
wo oT
2
2 2
1! 1!
w"
10 10
oF
+
is
3
3
°
5
7 a
®=%
°
5
7 a
~ =
freer)
Fine
freer
Fine
freer
Fine
0!
we
Ee
E
Ee
wo w'
“
ca
10° 10°
°
5
7 a
= %
5
7 a
®=%
°
5
7 a
fepany
Fine
fepany
Fine
fepany
Fine
Ee Ee Ee
wo
w'
w'
®e
5
7 a
® = % ®e
5
7
= % ®e
5
7
® =
Ey
frementy
Fine
frementy
Fine
frementy
Fine
126
ce
en
intron
£8
Wet
0-68.20
687
-
5
ie
o
ary
0
1% 20
wo
3
0 8
20
Es
wo
0 8
20
Es
—
maa maa
ro ro
0 8
20
wo
3
0 8
20
Es
wo
0 8
20
Es
maa maa maa
¢
ro ro
wo
0 8
20
wo
3
0 8
20
Es
wo
.
e
o24
oe
maa maa
amen
¢
ro
Ew
|
amen amen
aamean
¢
z
ro
s
Wo
ww
requ)
Fine
s
Wo
ww
requ)
Fine
remy
Fine
remy
Fine
127
wo
wo
5
Ew Ew
wo wo
10°
=
‘
wy
=
0
wy
=
0 6
fremeney
Fine
fremeney
Fine
fremency
Fine
oo
Ew Ew
5
1
wo wo
wy
=
0
‘
wy
=
wo
yt
3s
er
0
fremeney
Fine
fremeney
Fine
remency
Fine
.
oo
5
Ew Ew
1
10" 10"
10° 10°
o
+
@
3
«3
¢
788
x
+
ot
a
4
78
0
remensy
Fine
remensy
Fine
remensy
Fine
o
o
o
wo
ca
wo
= = =
10" 10" 10"
+
>t
se
3
x
rr
: 8
wo
yt
3s
er
0
remency
Fine
remency
Fine
remency
Fine
.
Ee Ee
Ew
10"
10"
10"
os
1
335
a5
4
ot
Sas
ase
eo
05
i
3
a5
3s
as
fremeney
Fine
fremeney
Fine
fremeney
Fine
wo wo
=
2”
=
10"
10" 10"
eo
05
1
335
a5 as
ot
Ss
eo
05
i
3
a5
3s
as
fremeney
Fine
fromeney
Fine
fremeney
Fine
128
w
2
sateete
1
e
5
6
e
5
6
wo
2
4
6 e 10
z
“
cs
w
ene
Fae
ene
Fae
omens
Fe
il
-
o
2
«
6 8
wi
il
i
r
remy
Fine
remy
Fine
remy
Fine
remy
Fine
remy
Fine
remy
Fine
i
freensy
Fine
o
2
«
6 8
wi
freensy
Fine
€ *
fremensy
Fine
w
=
sone
ww
1
1
a
ot or
renee
Finke
renee
Finke
renee
Finke
Ew
=
€ *
fremensy
Fine
129
co oo
a
= =
1 1
wo
wy
es
wy
ees
a
a
a
fremeney
Fine
remency
Fine
remency
Fine
“0
o
co oo
ca
=
2°
=
1
10
.|
10"
0
wy
yt
se
ree
wy
ees
a
remency
Fine
remency
Fine
remency
Fine
Ew
Ew
Ew
wo
wo
wo
wy
>
4
se
7 8
ex
wy
ees
a
a
a
remency
Fine
remency
Fine
remency
Fine
o
ca
5
Ew
wo
10"
x
2 ‘
$ °
0
wy
2 ‘
$
0
e
frome
Fine
frome
Fine
ca ca
= =
2°
cy
1
10"
0
10°
10°!
2 3
ese
o
+
2 3
4s
o
+
234
se
frees
Fine
frees
Fine
frees
Fine
o
ca
Ew
5
Ew
10"
wo
10°
12
3,
8 8
7
%
ec)
frees
Fine
130
ce
wo wo
2 =
10" 10"
1
2 3 ‘ 8
Wo
ost
15
2
25
3
35
4
45
5
Wo
ost
15
2
25
3
35
4
45
5
frequensy
Fine
requensy
Fine
requensy
Fine
wo
wo
Ew
=
= =
10"
10"
3
o
os
1
15
@
25
3
35
4
«5
5
Wo
ost
15
2
25
3
35
4
45
5
requensy
Fine
requensy
Fine
10!
wo
wo
=
Ew
= =
10"
10"
o
08
1
15
@
25
3
35
4
«5
5
Wo
ost
15
2
25
3
35
4
45
5
a
08
+
15
2
requensy
Fine
requensy
Fine
requeney
Fine
‘unt
Suhel
tho=2
Q0e19F00.0575
"
und
Stoxtaad
112
185
F120 0575
.
Runt0S
Suhel.
thO=2
00812
F
120.0575.
wo
= =
10"
=
10"
°
08
+
15
2 5
"0
08
+
15
2
25
a
08
+
15
2
requeney
Fine
requeney
Fine
requeney
Fine
wo
10!
Ew
gw
= 2
10"
10"
10"
08
+
15
2 5
08
+
15
2
25
o
2 #
6
6
® 2
wu
6 8
2»
requeney
Fine
requeney
Fine
"resueney
Fine
10!
wo
tw
=
= =
10"
10"
Cr
requ)
Fine
> 4
6
8
© @
ww
wm»
requ)
Fine
owe
requ)
Fine
131
en
en
en
tw Ew
Kove
|
oe
.
coc
ata ata ata
icc
..
o
e
10.
2
2
4
6 e 10
z
“
6
wo
2
4
6 e 10
z
“
6
w
ata
Bent Bent
tw Ew
i
=
|
her
Bent Bent Bent
5
Ew
Fy Fs
o
2
4
6 e 10
z
“
6
1"
2 6 e 10
z
“
6
e
2
4
6 e 10
z
“
6
w
Bent Bent
Bante
tw Ew
.
or
.
—
Bent Bent
masa
5
Ew
oo
.
a
—
1
o
2 4
0
wz
“
0
2 + 10
2
4
0
2 + 6 8 10
2
4
fremeney
Fine
132
o
«
«
«
E
Ew
5
|
«|
-
wl wl
a
CT
a
a a
a
So
eres
cereus
eres
«
«
Ew ge
z
sal
|
‘|
e
2 4 6 6
0
wz
e
2 4 6 6
0
wz
wo
6
0
2
eres eres
=
w
«
«
.
sess
:
=
.
cre
|
‘|
2
6 6
0
e
2
6
0
e
2
«
s
6
7ea
0
= =
rem
Fae
« F
«
z
ge gw
| |
wl wl
a
a
4
°
i
icra icra icra
«
«
Ew
ee
| |
«|
w
wl
wl
a
oa
°
ya
icra
einen
icra
«
z
ge gw
«| «|
w! wl
Fo
a
°
So
eeu eeu eeu
133
a
a
a
ai
a
aa
3 + 3 + 3 +
fremeney
Fine
fremeney
Fine
fremeney
Fine
Fins?
Shoue0
2h0at
5
00051016867
unt:
Sroua-a
et
$0024
F
1-0
oe887
Fant
Gueuhand
NOt
$COe19
10076867
ce
en
a
&
a
ce
3 +
is
is
fremeney
Fine
fromenny
Fine
fromenny
Fine
en
ce
en
5
08
¥
is
25
3
wy
os
1
is
z
25
3
wy
os
1
is
fromenny
Fine
fromenny
Fine
fremenny
Fine
“0 “0
Ew
gw gw
..
Ko
10"
wo wo
os
1
is
25
3
wy
5
0
= » wy
5
0
fromenny
Fine
frementy
Fine
frementy
Fine
o
ca
gw
=
gw
= = =
wo
wo
10"
10°
el
°
5
0
Ey
= »
5
0
Ey
= »
°
5
0
frementy
Fine
frementy
Fine
frementy
Fine
o
ca
=
gw gw
= = =
wo wo
10"
°
5
0
= » wy
5
0
= » wy
5
0
Ey Ey
frementy
Fine
frementy
Fine
frementy
Fine
134
en
a 4
w=
‘reaey
Fine
a
|
°
sateete
‘
uaatsees
:
*
;
«
2
4
6 e 10
z
“
6
1"
e 10
“
6
e
2
4
6 e 10
z
“
6
w
eats eats eats
#
*
«
”
:
077
os
: :
077
os
: :
077
os
eiicn eiicn
a
ae
:
Neues
:
—
#
*
” ”
:
777
os
: :
077
os
: :
077
os
eiicn eiicn
a
=aes5
:7
Recor
#
<a
*
” ”
:
077
os
y :
i
‘i
=
-
<a
_
eiicn eiicn
ata
.
°
*
«
e
10.
2“6
e
2 +
6
e
10.
2“6
1
2 wo
2 +
6
e
10.
2“6
2
eaten eaten eaten
#
*
j
«
z
135
8!
8
w
2 2
1 1
a
3 ‘
= :
i
y
3 ‘
= +
i
2
eae
Fe
teary
Fite
teary
Fite
we
2 2
1 1
3 ‘
=
°
+
y
3
y
3
8
teary
Fite
taney
Fite
taney
Fite
8!
oa
2
2
w'
wl
3
‘
y
3
8
taney
Fite
taney
Fite
taney
Fite
ry
aoe
—
|
fooae
w
w'
3
‘
y
3
wo
a
a
a
taney
Fite
taney
Fite
rears
we
Eo Eo
Ww Ww
3
ea
ee
a
wo
a
a
rears
Fite
rears
Fite
rears
Fite
1!
we
Ww
1
ta
a
se
ta
a
se
wo
a
rears
Fite
rears
Fite
rears
Fite
136
137
Surge force
“0
o
wo
ca
Ew
5
Ew
= = =
ry
wy
=
0
wy
=
6
=
6
fremeney
Fin
te
fremeney
Fin
te
fremeney
Fin
te
unztSuow0
7
hz
ona
4
0
4025
.
ur;
Sond
72
00
100
4028
urz60.
Sven
58
=2
0065
21003895
o
Ew
5
= =
wy
=
wy
=
0
15
fremeney
Fin
te
fremeney
Fin
te
ca
Ew
Ew
5
= = =
wy
z 7 @ ®
0
wy
z 7 @ ®
6 2
i“
wy
7@®
0
fremeney
Fin
te
fremeney
Fin
te
fremeney
Fin
te
«0 ae
ca oo
Ew Ew
Ew
= = =
ra ra
10
10
10°
°
z 7 @ ®
0 w
°
z 7 @ ®
0 w
i“
°
7@®
0
fremensy
Fine
fremensy
Fine
fremensy
Fine
oo
Ew
Ew
5
= = =
0
wy
z 7 @ ®
0
a
a
wy
a
a
fremeney
Fin
te
remency
Fin
te
remency
Fin
te
“0
o
wo
ca
ca ca
E =
Ew
= = =
0"
se
remnsy
Fine
138
i
i
:
cn Cn
remency
Fin
te
remency
Fin
te
remency
Fin
te
Oo
Oo
Oo
a
+ +
fromeney
Fine
fromeney
Fine
fromeney
Fine
:
/
|
7
8
|
7
8
re
fromeney
Fine
fromeney
Fine
tremens
Fine
er
Oo
er
o
+
3
34
.
78
3
a
fromeney
Fine
remency
Fine
remnsy
Fine
Oo
Oo
er
ot
Ws
eS
ot
Ws
eS
ots
es
35
remency
Fine
remency
Fine
remnsy
Fine
er
er
ie
23
+
remency
Fine
remency
Fine
remensy
Fine
139
at
Soiled
nod
ODe19F1=0.0875
@ «
#
i *
= €
10°
E E
e
t
|
——
|
|
eo
os
5 2 3
eo
os
+
15
2
cry
o
05
+
15
e
aasiene
ain
; ;
cana
cara
/
canarias
Sen
°
E10
& =
E E E
.
05
+
15
2 5
05
+
15
e
25
Wy
2 +
.
e 10 12
4 6
e
Ey
aasiene
ain aaa
.
°
Ew
&
Ew
E E E
#
- #
So Lo So
aaa aaa aaa
° ®
° °
= =
Ew
E E E
Lo Lo So
aaa aaa aaa
.
aacnereaneiaas
.
canna
pesaareat
/
aca
*
Ew
&
Ew
E E E
#
»
So Se
SS
eats
eerste eerste
#
°
°
°
&
Ew
&
E E E
rat
=
tat
=
rat
=
140
oo
OoererOo
10
z 7
6 w
Ew Ew
= =
oe
4 6
8
we
wo
a
a
ew
ee
eo
a 4
Se
frmeney
Fine
frmeney
Fine
frmeney
Fine
ca
Ew
5
= =
10°
°
z 7 @ ®
0
°
z 7 @ ®
0
i“
z 7 @ ®
0
fremensy
Fine
fremensy
Fine
fremensy
Fine
:
“0
oo
ca
=
Ew
= =
ra
°
z 7 @ ®
0
wy
z 7 @ ®
0
i“
wy
z 7 @ ®
0 w
fremeney
Fin
te
fremeney
Fin
te
fremeney
Fin
te
oo
Ew
5
= =
0
°
8 2
7
+
2 3
«3
¢78
8
fremensy
Fine
remnsy
Fine
o
+
@ 3
«
3@7@8
remency
Fin
te
o
ca
=
Ew
= =
10" 10°
o
+
2 3 :
7s
sw
o
+
2 3 * 8
w
o
+
2 9
a
a
os
remnsy
Fine
141
wo
=
Ew
=
= = =
0
0
10° 10"
°
7
23+
‘ °
7
23+5
‘
7
23+5
‘
7
fremeney
Fine
fremeney
Fine
fremeney
Fine
‘e wo wo
E E E
= = =
SS
ni
3 ‘
fremensy
Fin
te
3+5
‘
fremensy
Fin
te
arto
rou
{hoe
$ODe19
120076567
a a
#
«
Fe
fe
i
i i i
"
«
aie
ee
+
23+
6
ee
05
1
15
ery
ee
05
1
15
z
ery
3
Sasaivats
satan satan
w
2
w
2
set
"
‘
«|
‘
i
For
i
i
Po:
i
Ps ‘|
Se
05
1
15
ery
ee
08
1
15
ery
05
1
15
ery
3
satan satan satan
,
jacana
aaa
aaa
.
cai
orn
:
At
Sano
cso
e
«
: i
fe
i i i
‘
Pe
+
—5
#5
>
7
oi
s—>
satan
saateans
eon
.
fe
fe
i
i i i
ey
5
10
15
cy
ey
5
10
15
cy
ee
5
10
15
cy ~w
saateans saateans saateans
142
ry
ca
=
Ew
Ew
= = =
wy
°
5
= wy
°
5
wy
°
5
frementy
Fine
frementy
Fine
frementy
Fine
Ew Ew Ew
= = =
wy
050
wy
050
wy
050
fremeney
Fine
fremeney
Fine
fremeney
Fine
o
wo
Ew
5
Ew
= = =
wy
050
wy
050
wy
050
fremeney
Fine
fremeney
Fine
fremency
Fin
te
ca
=
Ew
Ew
= = =
wy
050
wy
050
wy
.
8
0 @
wwe
we
fremeney
Fine
fremeney
Fine
remenny
Fine
o
ca
Ew
5
Ew
= = =
wy
.
8
0
@
wwe
we
.
8
0 @
wwe
we
wy
.
8
0 2
www
remenny
Fine
remenny
Fine
remenny
Fine
1 o
cy
ca
wo
ca
E E E
= = =
wy
.
8
0
@
wwe
we
wy
sw
a
)
wy
sw
a
)
remenny
Fine
remenny
Fine
remenny
Fine
143
;
erer
-
fi
Oo
r
Ew
=
0
+o
frees
Fine
fremensy
Fin
te
fremensy
Fin
te
er
ca
E
=
fremensy
Fin
te
f
fremensy
Fin
te
fremensy
Fin
te
§
Ew
§
= = =
0
o
3
wo
3
0
18
a
reeney
Fine
reeney
Fine
reensy
Fine
o
wo
Ew
§ =
= = =
0
a
So
a
2 3
4
3678
8
o
1
2
ssw
reensy
Fine
reensy
Fine
wo
wo
0
§
Ew
Ew
= = =
0
0
+
2 3
«
3@788
+
2 3
«@
3¢788
Ww
+
2 3
«
3@788
freee
Fine
freee
Fine
freee
Fine
ants
Steud
ht
00224
F120.115
ani:
Suen
1
et
GOH
F80.18
ants
Stout
GOSS
F1e0.15
wo
o
wp
0
Ew
§
= =
0
|
I
-
.|
|
Wo
ost
a
oO
ostS
as
fremeney
Fine
fremeney
Fine
fremeney
Fine
144
o
wo
a
0
ee ow
= = =
= = =
co
a
sal
Wo
os
tS
Wo
os
tS
rr
fremeney
Fine
fremeney
Fine
fremeney
Fine
wo
ot
wo
§
Ew
§
= = =
5
0
ba
5
6
5
reeney
Fine
reeney
Fine
reeney
Fine
wo
wo
Ew
Ew Ew
= = =
0
o
3
wo
3
0
8
wo
3
0
reeney
Fine
reeney
Fine
fremensy
Fine
Funtee:
SuebaiO
7
hoe?
Oot.
4
04025
ants:
Sturte07
02
005
F1s0
4025
2
unt42
Sues
62
000671
FH0.345
wo
§
Ew
§
= = =
10"
40
0
°
3
0
°
5
6
o
2468
0 @
ww
freaeney
Fine
freaeney
Fine
freuen
Fine
wo
§
Ew
Ew
= = =
2468
0 @
ww
Wo
2 4 6 8 wo
@
Ww
w
wo
z a 6 @
6 7
fremensy
Fine
fremensy
Fine
reensy
Fine
wo
‘o!
wo
wo
wo
= §
Ew
= = =
0
10
0
40"
0
10
°
z a 6 @
6 7
1
°
2 + 6
0 7
1“
°
z + 6 ®
0 7
1"
reensy
Fine
reensy
Fine
reensy
Fine
145
erererOo
;
f
i
erOo
wo
=
=
°
z a 6 @
6 7
wo
z a 6 @
0 7
2 + 6 ®
0 7
1“
freensy
Fine
reensy
Fine
reensy
Fine
wo wo
= =
= =
0 0
o
2
#6
ea
a a
reensy
Fine
reensy
Fine
reensy
Fine
wo
10!
wo wo
Ew
§
= =
10)
ot
@ 3
«
$@7@8
re
a
a
reensy
Fine
reensy
Fine
reensy
Fine
‘o!
wo
=
Ew
= =
o
+
@ 3
«
3@7@8
o
+
@ 3
«
3@7@8
0
remency
Fin
te
remency
Fin
te
remency
Fin
te
wo
wo
E E
= =
0
°
2‘€
0
2‘€
0 s
2‘€
0
fremeney
Fine
fremeney
Fine
fremeney
Fine
Ew
5
= =
0
o
+
3
34
7
Se
1
ae
7
Se
1
ane
.
78
fromeney
Fine
fromeney
Fine
fromeney
Fine
146
er
Oo
er
i
:
er
a
3
tremens
Fine
re
tremens
Fine
Fant
12:
Seno
22
2
00339
F150
1265
er
Oo
a
3
tremens
Fine
a
remnsy
Fine
a
remnsy
Fine
Ew Ew
= =
0 0
a
wo
Wo
st
sas
Wo
st
sas
remency
Fine
remency
Fine
remency
Fine
un:
Sueno
220-2
008108
F1s0
1265
ant?
Seoul
2202.
GO-5F1=01285
.
ntO6
Sead
tho2
082240
0578
0
@
0
E
=
a
ws
a
a
is
7
ess
Fie
ess
Fie
ecpaey
Fae
we
i
aif
er
as
8
5
as
+
8
7
3s
a
is
7
fans)
Finke
fans)
Finke
fans)
Finke
s
tO
Sine
nO
00-5
100575
.
fn:
Seino
7h
50708
10
59687
0
eo
=
Ew
E
a
=
I
I
as
+ c
3s
"
as
+
8
7
3s
a
ecpaey
Fae
ecpaey
Fae
‘cones
Fine
8
147
0
wo
Ew
Ew
= =
0"
o
2 4
€
6
0 @
ww»
a
a
a
2 4
6
8
0 @
ww wm
remenny
Fine
remenny
Fine
remenny
Fine
wo
=
Ew
=
= = =
10°
o
2 4
©
6
0 @
ww»
o
2 4
©
8
© @
ww»
2 4
6
8
0 @
ww wm
remenny
Fine
remenny
Fine
remenny
Fine
inga
Stosai0
7
hot
500-5
F1=053667
.
RirS2
Sound
Oe
SOOW67.4
F048
‘rot
Sout
6x02
500-36
F150.48
o
ca
E E E
= = =
o
2 4
€
6
0 @
ww»
oa
4S
a
ar
a
a
remenny
Fine
frmeney
Fine
frmeney
Fine
Ew Ew Ew
= = =
10° 10°
oe
4 6
8
6
oe
« 6 8 wo
ww
oe
4 6
8
6
fremency
Fine
fremency
Fine
fremency
Fine
wo
wo oo
Ew
5
Ew
= = =
ry
10°
> 4 6 8 wo
ww
>468
ww
oe
468
0
6
fremency
Fine
fremency
Fine
fremency
Fine
1
aot
wo ca
E =
Ew
= = =
10° 10°
oe
«68
wo
ww
oe
«68
wo
ww
oe
4 6
8
6
fremency
Fine
fremency
Fine
fremency
Fine
148
er
7
er
co
ee
er
=
@ *
fremeney
Fin
te
@ ®
0
fremeney
Fin
te
@ ®
0
fremeney
Fin
te
i
ios
eee
@ *
fremeney
Fin
te
@ ®
0
fremeney
Fin
te
7
er
—
@ *
fremeney
Fin
te
z 7
0 w
@ *
fremeney
Fin
te
fremeney
Fine
if
en
Ke
2 2 2
ae
2 6
0
ey
2 6
0
2
ey
1
2 >
«
s67ea
0
a a
etn
Fo Fo Ew
« « .
a
oa
a
a
reeset
reeset reeset
«
2
=
staaseee
Are
ee
ee
1
2 >
«
7ea
0
ae
1
2 > e a
0
ee
1
2> 7ea
0
se
remency
Fin
te
149
un
Srae0.
22015
0633
9
F1n0.
18067
i
Oo
=
oo
=
3 ‘
fremensy
Fin
te
3 ‘
fremensy
Fin
te
7
2 5
‘
3 ‘
fremensy
Fin
te
if
Ew
=
0
7
23+5
‘
fremensy
Fin
te
3 ‘
fremensy
Fin
te
Fans:
Svante
a0e24
F160
06857
°
7
2 5
*
3 7
fremensy
Fin
te
Oo
ho
°
7
2 3
+
5
‘
fremensy
Fin
te
°
os
+
fremensy
Fine
25
fremensy
Fine
fremensy
Fine
°
oF
+
fremensy
Fine
25
oF
+
25
fremensy
Fine
=
rs
remy
Fine
oF
+
is
25
°
5
°
5
=
°
5
°
5
=
fremensy
Fine
remy
Fine
remy
Fine
Ew Ew
= =
°
5
°
= wy
5
°
Ey
= wy
5
°
Ey
=
remy
Fine
150
wo
Ew
§
Ew
= = =
ba
3
0
1
Es 7”
wo
w
6
=
EY
e
=
freensy
Fine
freensy
Fine
.
‘nt2
Srowahd
601
O0667
1
FIn063
.
Fant
Stounana
hos
0053615089
a
wo
wo
= §
Ew
= = =
0
re
o
2468
0 @
ww
2 6 8
0 2
ww
Ww
°
a
freuen
Fine
freuen
Fine
freuen
Fine
cy
wo
Ew
Ew
§
= = =
40"
40"
re
°
3
0 6
20
8
°
0 6
20
8
°
0 6
20
freuen
Fine
freuen
Fine
freuen
Fine
wo
§
Ew
Ew
= = =
10"
40"
re
°
3
0 6
20
8
°
0 6
20
8
°
0 6
20
freuen
Fine
freuen
Fine
freuen
Fine
wo
§
Ew
Ew
= = =
3
0 6
20
8 wo
06=
3 ba
6
8
0 @
wi
We
Treensy
Fine
Treensy
Fine
remensy
Fine
wo wo
§
Ew
§
= = =
0
0
oe
oa
4
6
6
0 4
wi
Ww
#
frees
Fine
8
frees
Fine
2
+o
frees
Fine
151
ererOoerer
wo
Ew
§
= =
re
o
2 4
6
8
0 @
ww
2
o
2 ¢
6
8
0 @
ww
2
reenny
Fine
freaeney
Fine
wo
Ew
§
= =
wo
z<6
0
i2
2‘6
0
frees
Fine
frees
Fine
wo
=
Ew
= =
0
0
°
2‘6
0
ba
5
6
Wo
5
i
frees
Fine
reeney
Fine
reeney
Fine
Ew Ew
= =
0
°
3
0
wo
3
0 6
5
i
6
reeney
Fine
reeney
Fine
reeney
Fine
Ew Ew
= =
0
5
0
5
0 6
+
2 3
«
3@788
reeney
Fine
reeney
Fine
freee
Fine
o
‘ee
= =
= =
ot
2 3
4
3
ssw
So
a
2 3
4
3678
8
So
a
2 3
4
3 7 8
8
reensy
Fine
reensy
Fine
reensy
Fine
152
153
Appendix E: Comparison plots for fish15AR4 and fish15AR6
CFD
vs
Tr.
Func.
Modal
35;
FinType:fish1SAR4
(ar12deg,
1m)
log
@
[rad/s]
CFD
vs
Tr.
Func.
Modal
35;
FinType:fish15AR4
(ar20deg,
1.5m)
10°
log
@
[rad/s]
CFD
vs
Tr.
Func.
Modal
35;
FinType-fish15AR4
(artSdeg,
1.5m)
10°
log
@
[rad/s]
log
|H(F)I
log
|H(F)I
35;
30)
25]
log
|H(F)I
35;
CFD
vs
Tr.
Func.
Modal
FinType:fish15AR4
(Sdeg,1m)
10°
log
[rad/s]
CFD
vs
Tr.
Func.
Modal
FinType-fish15ARé4
(ar
2deg,1.5m)
35;
10°
log
@
[rad/s]
CFD
vs
Tr.
Func.
Modal
30F
20}
FinType-fish15AR4
(Sdeg,
1.5m)
—+—
CFD
de
Tr.
Func.
Modal
10°
log
@
[rad/s]
154
35;
log
|H(F)I
log
|H(F)I
log
|H(F)I
CFD
vs
Tr.
Func.
Modal
FinType:fish15AR4
(ar20deg,2m)
10°
log
@
[rad/s]
CFD
vs
Tr.
Func.
Modal
35
FinType-fish15AR4
(art2deg,2m)
10°
log
@
[rad/s]
CFD
vs
Tr.
Func.
Modal
40,
FinType:fish15AR6
(ar20deg,1m)
10°
log
@
[rad/s]
CFD
vs
Tr.
Func.
Modal
35;
FinType:fish1SAR4
(ar15deg,2m)
10°
log
@
[rad/s]
CFD
vs
Tr.
Func.
Modal
35
FinType:fish1SAR4
(Sdeg,2m)
30)
25
log
|H(F)I
—+—
CFD
de
Tr.
Func.
Modal
10°
log
@
[rad/s]
CFD
vs
Tr.
Func.
Modal
40,
FinType:fish15AR6
(ar15deg,
1m)
log
@
[rad/s]
155
35;
y
R
log
|H(F)I
20
log
|H(F)I
log
|H(F)I
CFD
vs
Tr.
Func.
Modal
FinType:fish1SAR6
(ar12deg,1m)
CFD
vs
Tr.
Func.
Modal
log
@
[rad/s]
40,
Y
R
FinType:fish15ARG
(ar20deg,1.5m)
10°
CFD
vs
Tr.
Func.
Modal
log
@
[rad/s]
40,
35,
30-
ry
&
20+
FinType-fish15AR6
(ar12deg,1.5m)
10°
log
@
[rad/s]
40,
35]
log
|H(F)I
8
x
R
log
|H(F)I
log
|H(F)I
CFD
vs
Tr.
Func.
Modal
FinType-fish15AR6
(5deg,1m)
20)
10°
CFD
vs
Tr.
Func.
Modal
log
@
[rad/s]
40,
FinType-fish15AR6
(ar15deg,1.5m)
35]
30)
y
R
x
8
10°
CFD
log
@
[rad/s]
vs
Tr.
Func.
Modal
40,
FinType-fish15ARG
(Sdeg,
1.5m)
35]
8
8
x
R
10°
log
@
[rad/s]
156
157
158
Vita
Erdem Aktosun was born in Izmir, Turkey. After completing his schoolwork at Ahmet Yesevi High
School in Izmir, He entered Yildiz Technical University in Istanbul. He obtained his
Bachelor’s degree in Naval Architecture and Marine Engineering from Yildiz Technical University in 2010.
He joined the University of New Orleans Naval Architecture and Marine
Engineering graduate program to pursue a M.S in this program in 2013.
