• No results found

Modeling, Simulation and Dynamic control of a Wave Energy Converter

N/A
N/A
Protected

Academic year: 2021

Share "Modeling, Simulation and Dynamic control of a Wave Energy Converter"

Copied!
66
0
0

Loading.... (view fulltext now)

Full text

(1)

Modeling, Simulation and Dynamic control of a Wave Energy Converter

M A R I A B Å N K E S T A D

Master of Science Thesis Stockholm, Sweden 2013

(2)
(3)

Modeling, Simulation and Dynamic control of a Wave Energy Converter

M A R I A B Å N K E S T A D

Master’s Thesis in Numerical Analysis (30 ECTS credits) Degree Progr. in Engineering Physics 300 credits Royal Institute of Technology year 2013 Supervisor at CorPower Ocean was Patrik Möller Supervisor at KTH was Christina Carlsund Levin Examiner was Michael Hanke

TRITA-MAT-E 2013:55 ISRN-KTH/MAT/E--13/55--SE

Royal Institute of Technology School of Engineering Sciences

KTH SCI SE-100 44 Stockholm, Sweden URL: www.kth.se/sci

(4)
(5)

Abstract

The energy in ocean waves is a renewable energy resource not yet fully exploited. Research in converting ocean energy to useful electricity has been ongoing for about 40 years, but no one has so far succeed to do it at sufficiently low cost. CorPower Ocean has developed a method, which in theory can achieve this. It uses a light buoy and a control strategy called Phase Control.

The purpose of this thesis is to develop a mathematical model of this method—using Linear Wave Theory to derive the hydrodynamic forces—

and from the simulated results analyze the energy output of the method.

In the process we create a program that will help realizing and improv- ing the method further.

The model is implemented and simulated in the simulation program Simulink. On the basis of the simulated results, we can concludes that the CorPower Ocean method is promising. The outcome shows that the energy output increases—up to five times—compared to conven- tional methods.

(6)
(7)

Modellering, simulering och dynamisk kontroll av ett vågkraftverk

Sammanfattning

Vågenergi är en förnyelsebar energikälla som ännu inte utnyttjas fullt ut. Forskning inom konvertering av vågenergi till användbar elektricitet har pågått i cirka 40 år, men ingen har hittills lyckas att göra det till- räckligt kostnadseffektivt. CorPower Ocean har utvecklat en metod, som i teorin kan uppnå detta. De använder en lätt boj och en kon- trollstrategi kallad Phase Control.

Syftet med detta examensarbete är att utveckla en matematisk modell av metoden—genom att använda Linear Wave Theory för att härleda de hydrodynamiska krafterna—och från de simulerade resultaten analysera energiutbytet. Under arbetets gång skapades också ett simuleringspro- gram som hjälpmedel till att realisera och förbättra metoden.

Modellen implementeras och simuleras i programmet Simulink. Utifrån de simulerade resultaten kan vi dra slutsatsen att CorPower Oceans metod är lovande. Resultatet visar att energiutbytet ökar—upp till fem gånger—jämfört med konventionella metoder.

(8)
(9)

Acknowledgments

I would like to thank my colleges Patrik Möller and Gunnar Ásgeirsson on CorPower Ocean for their support during the project. I also want the thank Marco Alves from WacEC for helping me with his expertise.

My warmest gratitude goes to my supervisor Ninni Carlsund for her enthusiasm, interest end support during writing this thesis.

(10)
(11)

Contents

1 Introduction 1

2 Theoretical background 5

2.1 Basic concepts . . . 5

2.1.1 Buoy motion . . . 6

2.2 Wave theory . . . 7

2.2.1 Linear Wave Theory . . . 7

2.2.2 Ocean Waves . . . 11

2.3 Hydrodynamic Forces . . . 12

2.3.1 Radiation Force . . . 14

2.3.2 Excitation Force . . . 16

2.3.3 Hydrostatic Force . . . 17

2.3.4 WAMIT . . . 17

2.3.5 Drag Force . . . 18

2.3.6 The equations of motion . . . 18

2.4 The Power Take Off (PTO) . . . 18

2.5 The system of Partial Differential Equations for the WEC . . . 21

2.6 Phase control . . . 22

2.6.1 Resonant buoy oscillation . . . 22

2.6.2 Latching strategy . . . 23

3 Modeling 25 3.1 Modeling Approach . . . 25

3.2 Modeling of the buoy . . . 27

3.2.1 Simulink implementation . . . 27

3.2.2 Radiation Force . . . 28

3.2.3 Excitation Force . . . 29

3.2.4 Hydrostatic and Drag Force . . . 30

3.3 Implementing the PTO . . . 30

3.4 Linear damping . . . 32

4 Simulation results 35

5 Discussion 41

(12)

5.1 Summary . . . 43

Bibliography 45

Appendices 46

A Radiation Force Solution 47

B State-Space approximations 49

(13)

Nomenclature

h water depth

H wave height

Hs significant wave height

T wave period

Te energy wave period

λ wave length

ω angular wave frequency η surface elevation

r position vector

q generalized position vector g acceleration of gravity

ρ water density

p pressure

φ total velocity potential φ0 initial potential

φd diffraction potential φrad radiation potential

ˆφ(r) complex potential amplitude Sb buoy wetted surface

n normal vector to the surface Sb

˜n generalized normal vector to the surface Sb

Frad radiation force Fexc excitation force Fhyst hydrostatic force Mrad radiation torque Mexc excitation torque Mhyst hydrostatic torque

τrad generalized radiation force

µrad convolution part of the radiation force m infinite added mass

A(ω) radiation added mass matrix B(ω) radiation damping matrix

z(t) state vector in state-space approximation A˜ constant matrix in state-space approximation B˜ constant matrix in state-space approximation

(14)

CONTENTS

C˜ constant matrix in state-space approximation D˜ constant matrix in state-space approximation τexc generalized excitation force

τhyst generalized hydrostatic force τdrag generalized drag force CD drag coefficient

Sc cross section area of the buoy

mb buoy mass

mosc oscillator mass

xosc oscillator position vector

˙ξj velocity of flywheel j Jf w flywheel inertia

Pf rac,j energy absorption variable for energy device j Tgen,j generator torque for energy device j

u gear ratio

rpin pinion radius

Ftrans transmission force for energy devices in PTO p0u initial pressure in upper gas reservior

p0l initial pressure in lower gas reservior V0u initial volume in upper gas reservior V0l initial volume in lower gas reservior dcyl cylinder diameter in gas spring

lstoke length from oscillator bottom to oscillator equilibrium position Fgas force from gas spring in PTO

Fmoscg gravitational force from oscillatorin PTO Fµ force due to oscillator friction in PTO FP T O total PTO force

ωN natural period of the buoy in heave W T wire tension variable

2

(15)

Chapter 1

Introduction

In today’s energy demanding world and with increasing greenhouse gases from fossil fuel in the atmosphere, the need for renewable energy sources is greater than ever before. One that has not yet been fully utilized is the energy in ocean surface waves.

Wave energy originates from the wind, which in turn originates from the sun. When the wind blows over the ocean, friction gives rise to water movement and waves are generated. Even though the total energy of waves on earth is much lower than the total solar energy, it is much more dense, specifically wave energy is about five times more dense than solar energy; with power density up to 2–3 kW/m2 on the surface [8]. This makes wave energy a good renewable energy source, but only if we can find an efficient way to extract it. The study of converting the mechanical energy in the ocean waves to useful electricity began in the 1970’s [8], but still no one has found a way to make it profitable.

The device that converts ocean wave energy to electricity is called a Wave En- ergy Converter (WEC). A WEC have various shapes and sizes. The one we model in this thesis is a so called point absorber: it has a buoy that moves with the wave motion. The buoy is connected to an energy converter called Power Take Off (PTO), so that the movement of the buoy is transmitted to the PTO. Inside the PTO, a generator converts the mechanical energy to electrical energy and transmits it to the grid (see Figure 1.1).

A mathematical model of a product, in our case a WEC, assists during the de- veloping process. With the mathematical model we can test the behavior of the product before we have a physical product to test on. The objective of this thesis is to develop a simulation tool of a WEC developed by the company CorPower Ocean.

CorPower Oceanhas developed a WEC of the point absorbing type that aims to have higher energy output per capital cost, than those on the market today. The

(16)

CHAPTER 1. INTRODUCTION

Figure 1.1: A sketch of a Wave Energy Converter of point absorber type. A buoy is connected with a wire to an energy converter unit which converts the mechanical buoy movement to electricity.

special characteristics of the WEC is that the buoy is light compared with similar WEC and uses a control strategy called Phase Control. Figure 1.2 shows a photo of the CorPower buoy in scale one to thirty. The final buoy will have a height of around 9 m, a width of around 8 m and weigh around 25 tons.

Figure 1.2: A photo of the CorPower buoy. The buoy in the photo is thirty times smaller than the final buoy.

The highest power output is achieved when the buoy moves in resonance (phase) with the incident wave. The CorPower WEC uses Phase Control, which makes it possible to use it in an optimized way in almost all sea states (waves of various lengths and amplitudes).

The Phase Control is achieved by locking and unlocking the buoy motion at certain times. This is called latching since we are latching and unlatching the buoy.

2

(17)

The mathematical model of the WEC is developed in the numerical computing program MATLAB and the graphical simulation program Simulink. It consists of the mathematical model of the buoy and the mathematical model of the PTO merged into one system of partial differential equations. A controlling block is im- plemented in the Simulink model, where the latching algorithm for Phase Control is located. A graphical user interface is developed where the user can set the input parameters, start the simulation and view the simulation result.

The thesis consists of four different chapters following the introduction chapter.

In turn, they cover the theoretical background, the modeling and implementation, the results and finally a discussion.

Chapter 2 - Theoretical background

The second chapter provides the theoretical background of a WEC: the buoy motion the PTO and the control algorithm- phase control.

Chapter 3 - Modeling

The third chapter describes, using the theory from Chapter 2, the mathematical modeling of the WEC and the implementation in Simulink.

Chapter 4 - Simulation results

The forth chapter presents results and model verification.

Chapter 5 - Discussion

The final chapter discusses the result from the WEC model and its limitations. It also covers future areas of improvement.

(18)
(19)

Chapter 2

Theoretical background

This chapter goes through the theoretical background needed to understand the modeling of surface waves on the oceans.

The chapter covers four parts: the basic theory of ocean surface waves; the deriva- tion of main forces acting on the buoy, calculated from the theory of surface waves;

the Power Take Off (PTO), which converts mechanical energy to electricity; the theory of power absorption and the control strategy Phase Control.

Before we start with the theory of waves, we first define some basic concepts for the WEC modeling.

2.1 Basic concepts

In this thesis we only consider waves acting in two dimension, the x−z-plane. These waves are called plane waves, since their wave front moves in parallel planes. The coordinate system for the wave is defined such that the x-axis is in the direction of the wave propagation and the z-axis is vertical to the wave motion. The z-axis is zero where the water surface is located when the wave is at restl (see figure 2.1). The buoy is symmetric, so the (x, z)−axes can rotate to coincide with the plane move- ment of the waves. Hence, it is sufficient to model the buoy as two–dimensional.

We need some definitions to describe the waves: the wave length λ, the wave an- gular frequency ω and the wave period T . The wave height H defines as twice the wave amplitude (figure 2.1).

Waves can be classified as shallow water waves, intermediate waves or deep wa- ter waves, depending on the relationship between the wave length λ and the water depth h (see Figure 2.1). In this thesis we only consider deep water waves that occur when h/λ >> 1 [9].

(20)

CHAPTER 2. THEORETICAL BACKGROUND

z = 0

z= −h

λ z

η(x, t)

H h

Figure 2.1: Drawing of the main parameters of the wave: the wave length λ, the wave height H, the sea bottom depth h, and the surface elevation η(x, t).

To describe an ocean wave we need an expression for the free surface movement in the z-direction. The free surface is the surface between the air and the water and is by definition zero at z = 0. The free surface elevation, located at x at time t, is denoted η(x, t) (figure 2.1).

2.1.1 Buoy motion

The general buoy motion is defined by six degrees of freedom; the movement in x, y, z directions and the rotations around the x, y, z axes. By common notation we call these six degrees of freedom: surge, sway, heave, roll, pitch and yaw (see Figure 2.2).

x(surge)

y(swey) z(heave)

yaw

pitch roll

Figure 2.2: The coordinates used when simulation the movement of a buoy in a surface wave. The six degrees of freedom are (surge , sway , heave , roll , pitch , yaw).

In this thesis, only three degrees of freedom are considered, since the waves are 6

(21)

2.2. WAVE THEORY

moving in two dimensions and the buoy is symmetric. These are surge, heave and pitch, where surge is the x-movement, heave is the z-movement and pitch is the rotation around the y-axis. The derivations following this chapter could be done in the same way for six degrees of freedom, but we limit the derivation to three.

The position vector, (x, z) is denoted r. If we also include the motion in pitch we get the generalized position vector q, representing the position in three degrees of freedom.

2.2 Wave theory

We have to understand the ocean waves and their interaction with the buoy to model the buoy movement. This section covers the theory of surface waves and their interaction with an oscillating buoy.

Ocean waves arise when the wind is blowing over the ocean, forcing the water surface to move. The gravity is then working as a restoring force. Ocean waves can also originate from other sources such as earthquakes and Coriolis force. We only consider the wind generated waves since these are the waves with a frequency spectrum interesting for a Wave Energy Converter.

In this thesis, we consider both regular and irregular surface waves. A regular wave is sinusoidal and occurs if a steady wind blows during a long time period.

Ocean waves mostly consist of super-positioned regular waves, since the wind speed varies. A wave that consist of a range of regular waves with varying shapes is called an irregular wave.

The first part of this section goes through Linear Wave Theory. The second part describes how regular and irregular surface waves are modeled.

2.2.1 Linear Wave Theory

In this section we derive the equations that describes the motion and the pressure of the fluid, using Linear Wave Theory. The theory, also called Airy Wave Theory [5], seeks to describe the wave motion as a velocity potential with a few assumptions and approximations.

We assume the fluid to be incompressible (constant density) and inviscid (zero viscosity). This is an approximation of reality and to get a more correct model, a viscous force called drag force is later implemented. Furthermore, we assume the flow to be irrotational (no local rotation of fluid elements). We linearize the equa- tion by assuming higher order terms to be negligible.

We start by looking at two basic fluid dynamic equations: the Continuity equa-

(22)

CHAPTER 2. THEORETICAL BACKGROUND

tion and the Navier-Stokes equation [9]. For an incompressible flow the Continuity equation is

∇ ·˙r = 0, (2.1)

where r is the position vector (x, z). We use a simplified version of the Navier- Stokes equation—assuming an incompressible, inviscid and irrotational flow—called the Euler equation [9]

ρD˙r

Dt = ρ˙r

∂t + (˙r · ∇)˙r= −∇p + f. (2.2) The term D ˙rDt is the material derivative of the velocity, which is the rate of change of the velocity for a particle moving with the fluid. The term p is the total pressure of the fluid and f is the external force acting on the fluid, in our case the gravitational restoring force, f = −ρgez = −∇ρgz.

Velocity Potential

If a flow is both irrotational and incompressible; its velocity can be described with a velocity potential φ.

˙r = ∇φ ⇔ (∂φ

∂x = ˙x

∂φ

∂z = ˙z. (2.3)

When the velocity potential is used in the Continuity equation (2.1), we get the Laplace equation

∇ ·˙r = ∇2φ= 0. (2.4)

The Linearized Bernoulli equation

W replace the velocities with the velocity potential (2.3) in the Euler equation, (2.2)

∂φ

∂t +(∇φ)2 2 +p

ρ + gz

!

= 0, (2.5)

and then integrates the equation to get the Bernoulli equation [9],

∂φ

∂t + (∇φ)2 2 + p

ρ + gz = constant. (2.6)

We linearize the Bernoulli equation by eliminating the nonlinear term (∇φ)2

∂φ

∂t +p

ρ + gz = constant, (2.7)

and then rearrange it to

p= −ρ∂φ

∂t − ρgz+ constant. (2.8)

8

(23)

2.2. WAVE THEORY

We now look at the static case, which occurs when the water is still and ∂φ∂t = 0.

The pressure at the free surface η(0, t) = 0 is patm (atmospheric pressure). This gives us that the constant in (2.8) is patm.

p= −ρ∂φ

∂t − ρgz+ patm (2.9)

Here, −ρ∂φ∂t = dynamic pressure, −ρgz = hydrostatic pressure, and p0 = atmo- spheric pressure. We later use these pressures to calculate the forces acting on the buoy.

In conclusion, we have derived the Laplace equation of the velocity potential (2.4), and the linearized Bernoulli equation (2.9). Both these equations must be satisfied throughout the fluid, so they need to satisfy conditions on the boundaries.

Boundary conditions

We have four different boundaries to consider: the boundary at the sea bottom z= −h, the wetted surface of the oscillating buoy Sb (the buoy surface in contact with water), the sea bottom z = −h, and the free surface of the water η(x, t) (see figure 2.3).

z= −h η(x, t) Sb

n vn

∂φ∂z = ∂η(x,t)∂t on η(x, t)

∂φ∂n = vn on Sb (moving buoy).

∂φ∂n = 0 on Sb (fixed buoy).

∂φ∂z = 0 at z = −h.

Figure 2.3: Sketch of the boundaries for the velocity potential with its corresponding boundary conditions.

Let us first look at the free water surface η(x, t). There we have a so called kinematic boundary condition, which means that the fluid velocity component normal to the free surface, must be the same as the free surface velocity in the same direction [9].

It can be approximated to [5]

∂η

∂t = ˙z = ∂φ

∂z on η(x, t). (2.10)

(24)

CHAPTER 2. THEORETICAL BACKGROUND

The pressure of the fluid is constant at the free surface (atmospheric pressure, patm).

If we use the relation in equation 2.9 we get the boundary condition,

∂p

∂t = 0 ⇒ 2φ

∂t2 + g∂η

∂t = 0 at η(x, t). (2.11)

Now consider the wetted surface of the buoy Sb(see figure 2.3). No flow is permitted through the buoy surface. The fluid velocity perpendicular to the buoy surface is therefore equal to the normal component of the buoy velocity vn(see figure 2.3) [5].

This gives us the boundary conditions

∂φ

∂n = vn, on Sb, (2.12)

when the buoy is moving, and a special case when the buoy is fixed,

∂φ

∂n = 0, on Sb. (2.13)

The fluid velocity is zero at the sea bottom z = −h, so we have the boundary condition

˙z = ∂φ

∂z = 0, when z = −h. (2.14)

These boundary conditions will be extended with one at infinite distance from the buoy.

Summary of the linearized wave equations

To summarize, we have the linearized Bernoulli equation and the Laplace equation of the velocity potential,

p = − ρ∂φ

∂t − ρgz+ patm, (2.15)

2φ = 0, (2.16)

with boundary conditions,

2φ

∂t2 = − g∂φ

∂t, at the free surface η(x, t), (2.17)

∂φ

∂n = vn, on Sb, when the buoy is moving, (2.18)

∂φ

∂n = 0, on Sb, when the buoy is fixed, (2.19)

∂φ

∂z = 0, at the water depth z = −h. (2.20)

These equations are used when we derive the forces acting on the buoy. Before we derive the forces, the mathematical background for ocean waves is covered.

10

(25)

2.2. WAVE THEORY

2.2.2 Ocean Waves

Waves on an ocean are typically irregular and consists of multiple regular waves with varying frequencies and amplitudes, traveling on top on each other. Figure 2.4 illustrates the free surface elevation as a function of time for a regular and an irregular wave.

Figure 2.4: Example of a regular and an irregular wave.

First we describe the behavior of regular waves. The theory of regular waves is then used to describe irregular waves.

Regular Waves

Regular surface waves arise when the wind blows with constant speed for a long time over a large area. We let x be the direction of wave propagation and since we assume the waves to be two dimensional, x together with the time t are the only coordinates we need to describe the wave elevation. The surface elevation of a regular wave can be expressed with a sinusoidal function

η(x, t) = H

2 cos (ωt − kx) (2.21)

where H/2 is the wave amplitude, ω is the angular frequency, k = 2π/λ is the wave number and λ in turn is the wave length.

Irregular Waves

The period T and the wave height H are important properties of a wave. For irregular waves we use corresponding properties called significant wave height Hs— calculated by taking the average value of the highest one-third heights of the in- coming waves—and energy period, Te, which is the wave period of a regular wave, carrying the same amount of energy as the irregular wave (more information is found in [6]).

(26)

CHAPTER 2. THEORETICAL BACKGROUND

Figure 2.5: The surface elevation of two regular waves with different wave height an, angular frequency θn and phase shift ωn and the two regular waves super-positioned to an irregular wave.

The wave elevation for an irregular wave is described as a sum of waves of different shapes.

η(x, t) = XM

n=1

ancos(ωnt − knx+ θn) (2.22) Here, M is the number of waves added together, anis the wave amplitude, and θnis the phase shift of the nth wave. The phase shift takes a random value in the range [0, 2π]. Figure 2.5 shows an example of two super-positioned waves. The angular frequencies ωn, the energy period Teand the significant wave height Hs depends on the sea state at the geographic location we wish to simulate.

To describe a sea state mathematically, we use a so called energy density spectrum, which describes the amount of energy transported in the wave per wave angular frequency ωn. The spectrum depends on Hs and Te. An example of an energy spectrum is seen in figure 2.6, which shows the power density S(ω) as a function of the wave angular frequency ω. The surface elevation, calculated from the spectra is also shown, where the wave height an is obtained from the value S(ωn), taken from the spectrum.

In this thesis the JONSWAP (Joint North Sea Wave Project) spectrum [2], is used to model the irregular waves. It is derived from theory together with data from real measurements. The JONSWAP spectrum is commonly used for deep water waves (the interested reader can read more in [2]).

2.3 Hydrodynamic Forces

In this section we derive the different forces acting on the buoy, by using the equa- tions for the velocity potential and the pressure summarized in section 2.2.1.

12

(27)

2.3. HYDRODYNAMIC FORCES

Figure 2.6: To the left we see the energy density spectrum for a wave with significant wave height, Hs= 2.26 m and energy period, Te= 8.0 s. To the right we see the surface elevation of a wave, derived from the spectrum by equation (2.22), where an depends on the value Sn).

In section 2.2.1, we saw that the fluid motion is described by a velocity potential φ.

The velocity potential can be separated into three different potentials.

φ= φ0+ φd+ φrad (2.23)

The initial potential φ0 is the velocity potential in an incident wave if no buoy is present and the diffraction potential, φd, compensates for the presence of a fixed buoy in the incident wave. These two together represent the potential of a fixed buoy in an incident wave. The radiation potential, φrad, is due to the oscillation of a buoy in still water. The potential partition gives us a way to separate the forces acting on the buoy.

We assume when we calculate the different hydrodynamic forces that the oscil- lating buoy is initially at equilibrium position. We also assume that the sea bottom is sufficiently deep not to affect the free surface elevation. Irregular waves are super- positioned regular waves (2.22); thus, the derivations are limited to regular waves.

The velocity potential can be expressed as [5],

φ(r, t) = Rehˆφ(r)eiωti, (2.24) for a sinusoidal time varying wave. Here ˆφ(r) is the complex potential amplitude, r is the location vector, ω is the angular frequency of the incoming wave and t is the time. From now on we will use an accentˆ for all complex entities.

The forces on the buoy are derived from the fluid pressure on the buoy wetted surface Sb. The pressure is obtained from the Bernoulli equation (2.15) and by using the complex potential (2.24).

p= −ρ∂φ

∂t − ρgz= −ρRehˆφ0+ ˆφd+ ˆφrad

eiωti− ρgz. (2.25)

(28)

CHAPTER 2. THEORETICAL BACKGROUND

Here, the term patmis canceled out, since we assume that the buoy initially is in its equilibrium position. The forces are obtained by integrating the pressure over the surface Sb,

F=Z Z

Sb

pndS, (2.26)

where n is the normal vector to the surface Sb. The torque acting on the buoy is calculated with

M=Z Z

Sb

p(r × n)dS. (2.27)

We introduce a generalized force, which is a force vector in surge, heave and pitch.

τ = (F, M) (2.28)

By introducing the vector

˜n = (n, (r × n)), (2.29)

we can write the generalized force τ , as τ =Z Z

Sb

p ˜ndS. (2.30)

The three different potentials and the hydrostatic term −ρgz cause three different forces acting on the buoy: τrad, τexcand τhyst. We assume the fluid to be inviscid, which is a simplification of reality. Therefore, one additional force is added, to get a more accurate model. The force is called the drag force τdrag and is due to the friction in the water [9]. The forces are listed in table 2.1.

Table 2.1: The four hydrodynamic and hydrostatic forces.

Name Abbreviation Origin Description

Hydrostatic force τhyst −ρgz hydrostatic restoring force.

Radiation force τrad ˆφrad due to the radiating wave the oscil- lating buoy creates.

Excitation force τexc ˆφ0, ˆφd due to a fixed buoy in an incident wave.

Drag force τdrag due to the water friction.

The five following sections covers the calculation of these forces.

2.3.1 Radiation Force

Consider a still ocean. If we let a buoy oscillate in it, a wave will emerge; this wave gives rise to the radiation wave potential φrad.

14

(29)

2.3. HYDRODYNAMIC FORCES

The complex radiation potential is a superposition of the wave potentials in the three degrees of freedom [5]. This lets us write the radiation potential as

ˆφrad=

3

X

k=1

ˆukϕˆk= iω

3

X

k=1

ˆqkϕˆk, (2.31)

where ˆϕk is the so called unit-amplitude radiation potential, ˆuk is the generalized complex velocity vector and ˆqk is the generalized complex position in the k’th de- gree of freedom. The expression comes from the fact that the velocity in complex notation is written as

ˆuk= iωˆqkeiωt. (2.32)

The boundary condition ∂φ/∂n = vnon the boundary Sb for a moving buoy (2.18), can—using complex units— be written as

∂ ˆφrad

∂n = ∇ ˆφrad·˜n = iωˆq · ˜n = v · ˜n, (2.33) where ˜n is defined in (2.29). This gives us a boundary condition on Sb for ˆϕk:

∂ϕˆk

∂n = ˜nk (2.34)

The potential ˆϕ will also satisfy the Laplace equation (2.4) along with the free water surface boundary condition (2.17) and the sea bottom boundary condition (2.20).

The radiation potential needs one more boundary to be completed that says that a wave from the oscillating buoy goes to infinity.

The generalized radiation force can then be written as (2.30) τrad= (Frad, Mrad) = −ρZ Z

Sb

˜n∂φrad

∂t dS= ρZ Z

Sb

˜n Rehω2 ˆϕT ˆq eiωtidS.

(2.35) This equation is simplified into (see Appendix A for a detailed calculation),

τrad = −B(ω) ˙q − A(ω)¨q, (2.36)

where A(ω) and B(ω) ∈ R3×3, are two frequency depending symmetric matrices, and ˙q and ¨q are the generalized velocity and acceleration vectors. The matrix A(ω) is called the added mass matrix and B(ω) is called the damping matrix. If we let the frequency go to infinity in A(ω), we get the so call infinitely added mass m, which is a matrix depending on the infinity added mas in the different degrees of freedom. This equation represent the radiation force in frequency domain. If we instead look at the same force in time domain we get [12]

τrad= −m¨q −Z t

0

K(t − t0) ˙q(t0)dt0 = −m¨q + µrad. (2.37)

(30)

CHAPTER 2. THEORETICAL BACKGROUND

The matrix K(t) in the convolution term represents a memory effect due to the waves from the oscillating buoy.

The Fourier transform of the convolution integral µrad in equation (2.37) is

ˆµrad(ω) = −ˆu(ω) ˆK(ω). (2.38) where ˆu is the Fourier transform of the generalized velocity ˙q. This is an input- ouput system, where ˆK(ω) is the transfer function. The transfer function is [12]

Kˆ(ω) = B(ω) + iω [A(ω) − m] . (2.39) To get an approximation of (2.38) in time domain a state-space approximation is used [3].

˙z(t) = ˜Az(t) + ˜B˙q (2.40) µrad= ˜Cz(t) + ˜D˙q. (2.41) Here, z(t) is the state vector, and ˜A, ˜B, ˜C and ˜D are constant matrices (See Appendix B for a more detailed explanation). The accuracy of the state-space approximation depends on the number of states included in the state vector z(t).

The approximation tends to the true solution as the number of states goes to infinity.

2.3.2 Excitation Force

The excitation force is the force acting on a buoy fixed in an incident wave. To calculate the excitation force we consider the incident wave potential φ0 and the diffraction potential φd. The generalized excitation force is calculated by

τexc= −ρZ Z

Sb

˜nRehˆφ0+ ˆφd

eiωtidS. (2.42)

From the boundary condition (2.19) we know that the velocity potential has a Neumann boundary condition ∂φ∂n = 0 on the wetted surface Sb, when the buoy is fixed.

∂ ˆφ0

∂n = −∂ ˆφd

∂n, on the boundary, Sb. (2.43) We also have the boundary condition for the potential at the sea bottom (2.20) and the boundary condition at the water surface (2.17).

16

(31)

2.3. HYDRODYNAMIC FORCES

2.3.3 Hydrostatic Force

The hydrostatic force is also called the gravitational restoring force, since it always strives to reach the equilibrium of the buoy, at z = 0. The generalized hydrostatic force is calculated as

τhyst = −ρgZ Z

Sb

z˜ndS. (2.44)

Further calculations on the hydrostatic force will not be done in this thesis, the interested reader can find more information in [5].

The generalized hydrostatic force can be written as

τhyst = Sq, (2.45)

where S is the 3 × 3 hydrostatic stiffness matrix and q is the displacement from the equilibrium position.

We now have covered the three hydrodynamic forces acting on the buoy. The excitation force depends on the surface elevation and the shape of the buoy. The remaining hydrodynamic forces, the radiation force τrad and the hydrostatic force τhyst, depends on the shape and movement of the buoy.

The forces are calculated by solving the velocity potential. The velocity poten- tial is frequency dependent, as evident from the complex velocity potential (2.24).

There are standard tools used to calculate the hydrodynamic forces. One of them is WAMIT [1], which is the one we use in this thesis.

2.3.4 WAMIT

WAMIT (Wave Analysis At Massachusetts Institute of Technology) is a computer program used to analyze submerged bodies in presence of waves and is based on Linear Wave Theory [1].

The input to WAMIT is the shape and equilibrium point of the buoy. WAMIT calculates the velocity potentials in equation (2.23), with the associated boundary conditions (2.17-2.20). The frequency dependent forces are calculated using the velocity potentials as explained earlier in this chapter.

The excitation force is an array with the forces for the frequencies in a range of choice. The outputs for the radiation force are the matrices A(ω), B(ω) and m

(explained in section 2.3.1) and the output for the hydrostatic force is the stiffness matrix S (explained in section 2.3.3).

(32)

CHAPTER 2. THEORETICAL BACKGROUND

2.3.5 Drag Force

The drag force is due to the fluid resistance and its magnitude depends on the ve- locity of the buoy. The fluid resistance depends on the viscous friction in the fluid, which was neglected when the hydrodynamic forces were calculated.

The force is calculated by [9]

τdrag,k= 1

2ρ˙qk2Sc,kCD,k (2.46) where ρ is the water density, Sc,k is the cross section area of the buoy perpendicular to the motion in the k:th degree of freedom, and CD,k is the corresponding drag coefficient, which depends on the buoy shape and Reynolds number [9]. The drag coefficients are obtained from empirical measurements.

2.3.6 The equations of motion

The motion of the buoy is calculated by adding up the hydrodynamic and hydro- static forces.

mb¨q = τexc+ τhyst+ τrad+ τdrag, (2.47) Here, mb is the buoy mass. By using the relation in (2.37) it can be rewritten as

(Imb+ m) ¨q = τexc+ τhyst+ µrad+ τdrag, (2.48) where I is the identity matrix. This relation is used when we model the complete WEC.

The second part of the WEC is the Power Take Off (PTO) device, which is de- scribed in the following section.

2.4 The Power Take Off (PTO)

This section describes the different parts of the Power Take Off (PTO) device: the oscillator, two energy converting devices and the gas spring.

The motion of the buoy is transferred to the PTO through a wire, where an os- cillator moves with the buoy. The oscillator is connected to two energy converter devices with a gear (a cogwheel). One of the devices is used for the up motion and one is used for the down motion. The up motion energy converter is engaged during the oscillator up motion and free during the down motion, and the down motion energy converter is engaged during the down motion and free during the up motion (see figure 2.7).

18

(33)

2.4. THE POWER TAKE OFF (PTO)

Wire connected with buoy.

Oscillator.

Energy converter 1, engaged when oscillator is moving up.

Energy converter 2, engaged when oscillator is moving down.

Figure 2.7: Drawing of the oscillator with the energy conversion units. The wire is con- nected to the oscillator, transmitting the motion of the buoy to the PTO. The energy converter device 1 is engaged during oscillator up motion and energy converter device 2 is engaged during oscillator down motion.

The energy converter consists of a flywheel and a generator. The generator ex- tracts energy from the flywheel—which works as an energy storage—when best suited. The largest share of the energy is extracted from the flywheel to the gen- erator when the energy conversion unit is free. This decreases the damping on the buoy motion, since the flywheels can move almost freely during engagement. Figure 2.7 shows a sketch of the energy conversion devices and the oscillator.

The last parts of the PTO is an upper and a lower gas reservoir, which work as a gas spring for the oscillator (see figure 2.8). This gas spring protects the PTO from large forces, and works as a pre-tension force that puts the initial position for the buoy at its midpoint.

The free PTO (not connected with the buoy) is modeled by adding up the forces acting on the oscillator.

mosc¨xosc= k · Ftrans+ Fgas+ Fmoscg+ Fµ= FP T O (2.49) Here, moscis the oscillator mass, ¨xosc is the oscillator acceleration, Fgas is the force from the gas spring, Ftrans is a force from the energy conversion devices, Fmoscg is the gravitational force and Fµis a friction force. The variable k = 1 when an energy converting device is engaged, and k = 0 otherwise.

The energy conversion devices, when engaged to the oscillator, gives rise to the transmission force

Ftrans = −sign( ˙xosc) Jf w

¨ξju

rpin + Pf rac,jTgen,j· u rpin

!

, (2.50)

where the first term is the inertial force from the flywheel and the second term is the resistant force from the generator. Here, j = 1 when the energy conversion device one is engaged, and j = 2 when energy conversion device two is engaged. ˙xosc is the velocity of the oscillator, Jf w is the flywheel inertia, u is the gear ratio, rpin

(34)

CHAPTER 2. THEORETICAL BACKGROUND

xosc = 0

lstoke

Wire connected to buoy.

Wire connected to ground.

Vu pu

Vl pl dcyl

Figure 2.8: Drawing of the gas spring of the oscillator. The oscillator movement (up and down) changes the pressure and volume in the upper and lower reservoir and works as a gas spring.

is the pinion radius (the pinion connects the oscillator with the flywheel), Tgen,j is the torque from the generator j and Pf rac,j is a controller factor for generator j, deciding the amount of energy to be transferred from the energy conversion device to the grid. The energy conversion device, disengaged to the oscillator, is modeled with the equation

Jf w¨ξj = −Pf rac,jTgen,j, (2.51) where ¨ξj is the flywheel acceleration.

The gas spring gives rise to a gas force Fgas. It is modeled with the equation

Fgas =

p0l

V0l

V0l+πd42cylxosc+lstoke2 

γ

− p0u

V0u V0uπd

2 cyl

4

xosc+lstoke2 

γ

πd2cyl

4 . (2.52) The first term is the force from the lower reservoir acting upward on the oscillator and the second term is the force form the upper reservoir, acting downward on the oscillator (see figure 2.8). The terms p0u, p0l, V0uand V0l are the gas pressures and volumes when the oscillator is at the bottom position −lstoke; lstoke is the length from the bottom position to the equilibrium position for the oscillator. dcyl is the diameter of cylinder in which the piston moves. The term xoscis the position of the oscillator and γ is the specific heat ratio. The pressure of the lower gas reservoir is adjusted when the WEC is installed, so that the oscillator moves from −lstoke to its

20

(35)

2.5. THE SYSTEM OF PARTIAL DIFFERENTIAL EQUATIONS FOR THE WEC

equilibrium position at xosc= 0.

The weight of the oscillator gives rise to a gravitational force

Fmoscg = −mosc· g, (2.53)

where g is the acceleration due to gravity. The last PTO force is the friction force from the oscillator

Fµ= −µ · ˙xosc, (2.54)

where ˙xosc is the oscillator velocity and µ is a friction constant determined from empirical measurements.

We have now the equations for the buoy and the PTO. These are in the follow- ing section merged to one system of Partial Differential Equations.

2.5 The system of Partial Differential Equations for the WEC

We have derived the hydrodynamic forces acting on the buoy and the mechanical forces from the PTO. These forces together constitute a system of partial differen- tial equations that describes the motion of the whole WEC. The buoy is modeled in three degree of freedom, while the PTO is modeled in one degree of freedom.

The PTO and the buoy are connected by a wire. This gives us two different situa- tions: one when there is wire tension and the PTO is connected to the buoy, and one when there is no wire tension, so that the buoy and the PTO are two different systems modeled separately. The two equations (2.48) and (2.49) for the buoy and the PTO are when there is no tension:

(Imb+ m) ¨q = τexc+ τhyst+ µrad+ τdrag, (2.55)

mosc¨xosc= FP T O. (2.56)

The PTO is connected to the buoy when the wire is tense and the movement of the oscillator and the buoy is the same. Then the system is modeled as

(I(mb+ mosc) + m) ¨q = τexc+ τhyst+ µrad+ τdrag+ FP T O (2.57)

¨xosc= ¨xbuoy, (2.58)

where FP T O is the PTO force, split into the components of the coordinate system of the buoy, and ¨xbuoy is the acceleration of the buoy projected on the direction of the wire.

The last component we add to the model is a strategy to control the buoy mo- tion, with the aim to increase the power output. This strategy is called Phase Control.

(36)

CHAPTER 2. THEORETICAL BACKGROUND

2.6 Phase control

The goal for a WEC is to absorb as much of the energy in the waves as possible and convert it to electricity. This section describes the strategy we use to achieve optimal power absorption.

One way of achieving this is to use the control strategy Phase Control by latch- ing. This control strategy acts by locking and unlocking the motion of the buoy at certain moments, to force the buoy to move in phase with the incident wave. This maximizes the buoy velocity and thereby maximizes the energy output. Figure 2.9 illustrates the motion of the buoy when Phase Control with latching is used.

Figure 2.9: The surface elevation and the buoy motion when latching is used, when only looking at heave.

2.6.1 Resonant buoy oscillation

A buoy oscillating in an ocean has a natural frequency. The natural frequency is the frequency with which the buoy oscillates in still water. The natural angular frequency of a buoy in heave (2.2) is

ωN =

s ρgA

mb+ m∞,heave

, (2.59)

where g is the acceleration of gravity, A is the buoy plane area, mb is the buoy mass and mis the infinitely added mass for the heave coordinate (due to the movement of the surrounding fluid as discussed in section 2.3.1).

The optimal movement of a buoy (maximal velocity and position) occurs when the natural frequency of the buoy is the same as the frequency of the incoming wave [5], so that the buoy is at resonance with the wave. This typically requires a

22

(37)

2.6. PHASE CONTROL

large and heavy buoy to match the relatively low frequency of the wave. Since the frequency of the waves vary, a buoy must regulate its mass mb or the buoy plane area A to be in resonance.

There is an alternative way to achieve the resonance effect, without regulating the buoy mass and shape. The idea is to have a small buoy with natural frequency higher than most wave frequencies, and force the buoy to move in phase with the wave by locking the motion of the buoy at certain times. The buoy can then move with its natural frequency, even though it is much higher than the frequency of the incoming wave. The control strategy of locking the buoy motion is called Phase Control by latching. Figure 2.10 shows the movement of a buoy having the same natural frequency as the one of the incoming wave, together with a buoy with higher natural frequency than the incoming wave, but forced to move with its natural fre- quency using latching.

a) t

b) c)

a) Surface elevation of incident wave b) Movement of light buoy,

with phase control.

c) Movement of heavy buoy, with optimal mass and size.

Figure 2.10: Drawing adapted from [5] showing the movement of a heavy buoy, with the same natural period as for the incident wave, along with the movement of a small buoy, using the control strategy latching to achieve the same phase as the incident wave.

If we choose a buoy light and small enough we can force the buoy to move with its natural frequency for almost all incoming waves, since the frequency always is larger than the wave frequency.

2.6.2 Latching strategy

For a regular steady sinusoidal wave it is an easy task to predict when to lock and unlock the motion to get the optimal power absorption. For an unsteady or irreg- ular wave on the other hand, we need a way to predict when this should occur.

The locking of the motion occur when the buoy velocity is zero. The unlocking is more trickier to decide, so we need a way to predict when the optimal timing is.

One way to do this is to have a threshold for the surface elevation. When the wave passes the threshold, we unlock the buoy. The optimized threshold for the surface elevation was calculated from tests results in [10] to

ηthres= H 2 sinπ

2

1 −TN T



, (2.60)

(38)

CHAPTER 2. THEORETICAL BACKGROUND

where TN is the natural period of the buoy, H is the wave height and T the wave period. For irregular waves we can replace the wave period T with the average energy period Te and the wave height H with the significant wave height Hs. In the model the wave height is read from time series vectors of the surface elevation and in real life the wave height is estimated by sensors placed close to or on the buoy.

The method of wave elevation threshold works well for steady waves where the average energy period Te and the significant wave height Hs stays constant. For waves changing much with time, a way to predict the wave appearance—before it reaches the buoy—is needed. There are several ideas on how to estimate the wave shape before the actual wave reaches the buoy. One of them is to use sensors to sample the wave shape at a few location before the wave reaches the buoy. The wave shape at the buoy can then be estimated e.g, by a Kalman filter. Wave prediction is not included in this thesis, so we assume that Teand Hsare constant when irregular waves are modeled.

24

(39)

Chapter 3

Modeling

This chapter covers the numeric modeling of the WEC where the three sections covers: the modeling approach and some assumption that was made to simplify the model; the modeling of the buoy, with the hydrodynamic forces acting on the buoy;

the model extension to also include the PTO system, giving us the final model of the whole WEC.

3.1 Modeling Approach

The model is restricted to ocean waves moving with two degrees of freedom, x and z, where the surface elevation is denoted η(x, t). Because of the symmetry of the buoy it can be considered to move in two dimensions as well. When modeling the buoy we use the degrees of freedom: surge, heave and pitch (see figure 2.2), where surge and heave are the motions in x and z-direction and pitch is the rotation around the y-axis (see figure 3.1).

Figure 3.1: A sketch of the degrees of freedom used when simulating the buoy movement in a wave (see figure 2.2).

The hydrodynamic forces acting on the buoy are derived from Linear Wave The- ory (see section 2.2). We use generalized forces to represent the forces acting on the buoy, which include the forces in surge and heave, and the torque acting in pitch. The generalized hydrodynamic forces are obtained using the computer pro-

(40)

CHAPTER 3. MODELING

gram WAMIT (see section 2.3.4).

We use Simulink [14] as the simulation program to model the WEC. Simulink is preferred to regular MATLAB code since it better overviews the model, which in turn makes it easy for other users to understand the model and develop it further.

To build a GUI (graphical user interface), we use the MATLAB user interface tool GUIDE [11], which communicates with the Simulink model. The GUI is developed to make the model easy to use.

The numerical solver we use is ode45, since it is both robust and relatively fast.

ode45is provided by MATLAB and is based on the Dormand-Prince Runga-Kutta (4-5) formula [13].

Overview of the Simulink model

The mathematical model of the WEC is implemented in Simulink. The model includes a block calculating the hydrodynamic forces (see section 2.3), a block cal- culating the PTO forces (see section 2.4) and a controller block that handles the PTO settings and the Phase Control. Figure 3.2 overviews the complete Simulink implementation.

Switching states block

1 s

1 s

FP T O block

τhydro

block

Controller

¨x ˙x x

FP T O

τhydro

Figure 3.2: A sketch of the complete WEC model in SIMULINK. The position vector x includes the buoy q, the oscillator xosc and the flywheels ξi. The 1/s block represent an integration.

In the following section the implementation of the buoy with the hydrodynamic and hydrostatic forces is covered.

26

(41)

3.2. MODELING OF THE BUOY

3.2 Modeling of the buoy

The modeling of the buoy comes down to the generalized hydrodynamic and hydro- static forces acting on the buoy: τrad, τexc, τhyst and τdrag (see section 2.3 ).

In this section we describe the model implementation of these forces. Before we go into detail, we first overview the complete Simulink implementation of the buoy.

3.2.1 Simulink implementation

Simulink is a graphical modeling and simulation program, where blocks are used to create a block diagram of a model. Blocks included in Simulink are e.g. the integrator block, used for integration and the gain block, used to multiply an input with a constant or matrix. The equation of motion for a buoy oscillating in a wave (2.48),

(m+ Imb)¨q = τexc+ τhyst+ µrad+ τdrag, (3.1) is implemented in Simulink. An outline of the Simulink model is seen in figure 3.3 where the excitation force is the input and the driving force of the model.

Figure 3.3: An overview of the Simulink model of the buoy. The buoy acceleration ¨q is calculated by dividing the force vector with the buoy mass. The acceleration vector is then integrated one time to get the velocity ˙q and then one more time to get the position vector q. The position and velocity vector is then used to calculate the forces (see equation 3.1).

In the following sections the modeling of the forces are explained in more detail.

(42)

CHAPTER 3. MODELING

3.2.2 Radiation Force

The derivation of the generalized radiation force was done in section 2.3.1, where we saw that the radiation force can be written as (2.37)

τrad = −m¨q + µrad (3.2)

The infinitely added mass m is constant, and it can therefore be seen as a part of the buoy mass mb. The infinitely added mass is an output from the computer program WAMIT.

The computer program WAMIT gives the added mass matrix A(ω) and the damping matrix B(ω) explained in (2.36), as an output. The frequency dependent transfer function ˆK(ω) can then be calculated by (2.39).

The output-input system in (2.38), with the transfer function ˆK(ω), can in time domain be approximated with a state-space approximation [3] (see Appendix B for a detailed explanation).

˙z(t) = ˜Az(t) + ˜B˙q(t) (3.3) µrad= ˜Cz(t) + ˜D˙q(t) (3.4) The state-space approximation, for all three degrees of freedom, is calculated with two toolboxes: the SS Fitting toolbox [4] that provides a state-space model from the WAMIT data and the MSS FSI (Marine Systems Simulator) [7] that identifies the fluid memory effects of the buoy. The inputs to the toolbox are the WAMIT data: m, A(ω) and B(ω) .

When we have the state-space matrices in equation (3.3-3.4), we calculate the ra- diation force with the Simulink state-space block (see figure 3.4).

Figure 3.4: The Simulink State-Space block which takes the buoy velocity ˙q as an input and calculates the convolution part of the radiation force µradby a state-space approximation (see equation (3.3-3.4).

28

References

Related documents

Submitted to IEEE Transaction on Energy Conversion, March 2007 • Paper XI: Simulation of a Linear Generator for Wave Power Absorption - Part II: Verification [72] In this paper the

WEC Wave Energy Converter MPC Model Predictive Control PTO system Power Take-Off system LTI Linear Time Invariant...

In this project, I present a design of a scale model of a linear generator (LG) similar to a full size Wave Energy Converter (WEC) being developed at Uppsala University since 2002

In this thesis a Linear and a Nonlinear Model Predictive Controller have been developed with the goal to maximize the energy output of a Wave Energy Con- verter. As a

electromag- netic waves and wavelets for impulses (outbursts), can be used to model the above regional distribution of semantic content [5, 6], also demon- strated for the visible

If distant shadows are evaluated by integrating the light attenuation along cast rays, from each voxel to the light source, then a large number of sample points are needed. In order

Other renewables Climate change (Kg of CO 2 equivalent) (Pehnt, 2006).. Figure 4.3 shows a representation of different renewable energies Kg CO 2 per kWh including wave energy at

"Body is an Experiment of the Mind" is an intimation of corporeality; a thought that describes the self as a socially and environmentally vulnerable concept of body and