Taming Dragons Improving power performance of tethered turbine tidal farms Master’s Thesis within the Sustainable Energy Systems programme ERIC SINGHROY Department of Energy and Environment Division of Power Electrical Power Engineering CHALMERS UNIVERSITY OF TECHNOLOGY Göteborg, Sweden 2017 MASTER’S THESIS Taming Dragons Improving power performance of tethered tidal turbine farms Master’s Thesis within the Sustainable Energy Systems programme ERIC SINGHROY SUPERVISOR: Prof. Yujing Liu EXAMINER Prof. Yujing Liu Department of Energy and Environment Division of Heat and Power Technology CHALMERS UNIVERSITY OF TECHNOLOGY Göteborg, Sweden 2017 Master’s Thesis within the Sustainable Energy Systems programme ERIC SINGHROY © ERIC SINGHROY, 2017 Department of Energy and Environment Division of Heat and Power Technology Chalmers University of Technology SE-412 96 Göteborg Sweden Telephone: + 46 (0)31-772 1000 Chalmers Reproservice Göteborg, Sweden 2017 I Master’s Thesis in the Sustainable Energy Systems programme ERIC SINGHROY Department of Energy and Environment Division of Heat and Power Technology Chalmers University of Technology ABSTRACT Wright a clear and consist abstract. 1 Contents ABSTRACT I CONTENTS 1 LIST OF FIGURES 3 OBJECTIVE 6 1 INTRODUCTION 7 2 TIDAL CURRENTS BACKGROUND 9 Tidal origins 9 Tidal energy source 10 Daily tidal fluctuation 11 Monthly tidal fluctuations 12 Other sources of tidal fluctuations 13 Relationship between tide height and tide current 13 3 TUSK PLANTS BACKGROUND 15 Plant farm array 15 Synchronous start up 15 Current research on TUSK tidal farms 16 4 MODELLING THE PLANT TRAJECTORY 18 Other trajectories 19 Bernoulli’s lemniscate described in three dimensions 19 Plant Trajectory 22 5 POWER CURVE MODELLING 23 Optimal power derivation 23 Power calculations 26 Definition of world frame of reference, 26 Definition of plant body frame of reference, 27 Rate of change of yaw, roll and pitch 28 Torque modelling 29 Modelling the tidal stream 29 Angle of attack 30 Lift and drag coefficient 30 2 5.9.1 Lift coefficient 31 5.9.2 Drag coefficient: 31 Affective velocity 32 Angle of inclination 32 Angle of attack and angle of inclination 33 Fixed angle of attack and variable angle of inclination 34 Relative velocity to the stream current 35 Variable angle of attack and fixed angle of inclination 36 Turbine velocity comparison 37 Power loss from pitch with angle of inclination torque 38 Final power curve 38 6 MITIGATING POWER FLUCTUATIONS 40 Fourier analysis 40 Fluctuation index 41 Signal shifting 41 Fourier analysis method 42 Iterative method 42 Harmonics analysis of basic curve 42 Two plants optimization 43 Three plants optimization 44 Two pairs 45 Benefits of reducing fluctuations for tidal energy 47 Future work 48 7 CONCLUSION 49 8 REFERENCES 50 3 List of Figures Figure 1 Tethered undersea kite system and components (not to scale) 7 Figure 2 Rendering of the Minesto Deep Green TUSK system [6] 7 Figure 3 Forces of gravity exerted by the moon onto the large bodies of waters 10 Figure 4: Accumulation of traction forces to create tidal bulges [9] 11 Figure 5 Mean Daily Tidal Range [10] 11 Figure 6 Diurnal Tides [11] 12 Figure 7 Semidiurnal Tides [11] 12 Figure 8 Mixed Semidiurnal [11] 12 Figure 9 Astral Alignment for the Spring and Neap Tides 13 Figure 10 Relationship between tidal heights and tidal currents [15] 14 Figure 11 Seven plant TUSK array, top view 15 Figure 12 Synchronous start-up of a three plant tidal farm, side view 16 Figure 13 Bernoulli's Lemniscate (a=1) 19 Figure 14 Example of unit circular trajectory 19 Figure 15 Bernoulli's Lemniscate compared with Viviani's Curve 20 Figure 16 Depiction of altitude angle 20 Figure 17 Depiction the intersection of a unit sphere and cylinder (left) and then applying the rotation matrix (right) 21 Figure 18 Orthogonal Projections of unit Viviani's Curve 21 Figure 19 3D simulation of base trajectory 22 Figure 20 Orthogonal projections of Viviani curve 22 Figure 21 Methodology map from a given trajectory to power curve output 23 Figure 22 Forces and velocity triangle on a crosswind TUSK (similar to [20]) 24 Figure 23 Optimal drag ratio for maximum power output 25 Figure 24 World Frame of Reference 27 Figure 25 Plant F.o.R. with respect to World F.o.R. 27 Figure 26 Viviani's curve with interpolation points 27 Figure 27 Example of a local coordinate system in 2-D 28 Figure 28 Depiction of changes in Pitch, Roll and Yaw 28 Figure 29 Changes in pitch, yaw and roll angles for one cycle 29 Figure 30 Vector field representing the tidal current stream 30 Figure 31 Angle of Attack α of the plant 30 Figure 32 Angle of Attack for NACA 0015 [29] 31 Figure 33 Angle of attack curve fit using equations (30) and (31) [18] 31 4 Figure 34 Relation between lift and drag co-efficient for a wide range of angles of attack 32 Figure 35 Angle of attack for one cycle 33 Figure 36 Lift to Drag Ratio plotted with comparison to lift efficient plot and drag coefficient plot. 33 Figure 37 Modelon model of angle of attack and how it varies along the wing [24] 34 Figure 38 Variation of the TUSK angle of attack over the trajectory 34 Figure 39 Adjusted angle of attack with adapted angle of inclination 35 Figure 40 Affective Velocity for one cycle 35 Figure 41 Mapping of plant orientation onto the trajectory 36 Figure 42 Adjusted angle of attack with fixed angle of inclination 36 Figure 43 Lift to drag ratio with the range of operation 37 Figure 44 Affective velocity for one cycle with a fixed angle of inclination 37 Figure 45 Comparison of turbine velocity and means 37 Figure 46 Effective pitch angular velocity 38 Figure 47 Power Curves progression in complexity 38 Figure 48 Base power curve for one cycle 39 Figure 49 Fourier transform from time domain to frequency domain (Modified from [30]) 41 Figure 50 Example power curves and properties 41 Figure 51 Power curve of a complete cycle 42 Figure 52 Harmonic analysis of base power curve 42 Figure 53 Power Curves of a synchronized pair compared with the optimal pairing 43 Figure 54 Optimal pairing power curve 43 Figure 55 Paired plants power curve analysis for each iteration 44 Figure 56 Frequency energy of the optimized pairing power curve 44 Figure 57 Power curves of synchronized triplets and optimized triplets 44 Figure 58 Optimum triplets off set 45 Figure 59 Power curve of triplets at optimum off set 45 Figure 60 Fourier analysis of the triplets optimized power curve (unlabled and labled) 45 Figure 61 Power curves of two plant pairs 46 Figure 62 Fluctuation index for shifting two plant pairs 46 Figure 63 Power curve of the optimum shift between two plant pairs 46 Figure 64 Fluctuation Index of each plant combination 47 Figure 65 Tidal farm with multiple plant pairs and a plant triplet 47 5 Figure 66 Tidal farm with 4 plants pairing and a plant triplet 47 6 Objective The objective of this thesis is to investigate and mitigate power fluctuations of tethered undersea kite (TUSK) system by modelling a TUSK plant power output for a given trajectory and suggest a co-ordination strategy for several plants in a farm. This involves synchronising the operations of two or more plants by time shifting a second plant to compliment the first plant. The shifting of a plant relative to another momentarily disturb power production. Multiple plants can be coupled to optimize the shifting required while minimising power disruption. Therefore, the collection of plants, known as a tidal farm, will produce a more regular power output during production on par with traditional marine turbines. 7 1 Introduction Ocean energy is a virtually untapped source of renewable energy. This is expected to change with the concentrated effort in development underway in many countries around the world [1]. Canada, France, Sweden, the United Kingdom's and the United States are all developing marine technologies that are maturing to full scale testing phase [2] [3]. These technologies can contribute to meet environmental policy goals against climate change and limit the adverse effects of fossil fuels use while producing power. One of these recent technological developments is led by Minesto AB, a company based in Gothenburg, Sweden. Minesto has been researching and testing Tethered Under-Sea Kite (TUSK) systems [4]. The TUSK system produces electrical power from a hydrokinetic turbine submerged in the tidal current. The power plant consist of a turbine and generator fixed to a large rigid winged kite. The plant is anchored either to the sea floor or to a floating barrage by a flexible tether. The tether and winged kite allow the plant to move in a cross current figure 8 pattern. Figure 1 is a simplified diagram of the TUSK system anchored to the sea bed. The plant experiences velocities 5 to 10 times greater than the tidal current. This increases in velocity is a distinct cost effective advantage for TUSK systems compared to fixed marine turbine technologies [5]. Figure 1 Tethered undersea kite system and components (not to scale) Figure 2 Rendering of the Minesto Deep Green TUSK system [6] Multiple TUSK plants are arranged as an array on the sea floor and are electrically grouped to exported power to a single point of connection with the local distribution grid. These tidal arrays are expected to produce power in the Megawatts range [4]. Tidal currents have natural fluctuation patterns that are predictable. This predictability gives it a distinct advantage compared to weather based renewable energy sources [1]. Yet, TUSK systems have an addition source of power fluctuations. As it progresses through the figure 8 path, the plant must slow down to make the turns and accelerates 8 during the straighter portion of the loop. Therefore, the plant experiences changes in speed while completing a loop. Since power output is related to speed, the plant will produce fluctuating power as a consequence. This phenomenon is modelled and will be explained in greater detail in Section 4: Methodology. The power fluctuations are compounded when multiple turbines are synchronized. Two or more plants are synchronized when they are at the same position on their paths and will experience the same accelerations at the same time. Therefore, the synchronized plants will exaggerate power peaks and dips, i.e. the variance in power production. This greater swing in power furthers the difficulties of balancing demands. 9 2 Tidal Currents Background Power generation from TUSK plants is not dispatchable and will fluctuate with the tides. Power fluctuation is the increase and decrease of the electrical current feeding into the power grid. Fluctuations in power supply are not desired and difficult to handle on grids with limited inertia [7]. Ideally, the power supply must be balanced with power demand. Unfortunately, the oscillating patterns found in tidal kite power generation do not follow demand load variations. The following section describes the sources of fluctuations for tidal power generation and the varying time scale and the magnitude. This is useful to understand to determine an appropriate solution and the feasibility of a steady power supply. TUSK tidal power generation will oscillate on time scale of minutes to yearly scales. Other sources of renewable energy also suffer from this variable power production. Yet, the important difference is that with sufficient understanding of the various sources of fluctuations described in this section, the power supplied from tidal generation can be predicted reliably [8]. Currents from the thermo-halin belts, river flows, floods and meteorological conditions can cause temporary currents that interfere with tidal current patterns. These currents can be important locally, but will not be expanded on in this thesis. Tidal origins As the rain forests are the lungs of the earth, the ocean’s tides are the pulse of the earth. The tides are a continuous force seen all over the world. The tides are the periodic rise and fall of sea levels of large bodies of water from the changes in gravitational forces of the Moon and the Sun upon the spinning Earth [9]. All bodies of water experience tides. But only largest bodies have enough collective continuous surface area to have measurable tides. As the tides rise and fall, vast amounts of water must travel around the world to follow the Moon and the Sun. The tidal currents are simply the transport of this water. It is important to distinguish the vertical movement of the tides and the horizontal movement of the tidal current, since the relationship between them is complex. Tidal plants operators are concerned with the tidal current behaviour for power production. The power extracted from the tidal current will vary over time, magnitude and placement. The largest source of power fluctuation is the oscillations of the tides themselves. As the tides rise and falls, the tidal current will reach its maximum velocity, slow down to a halt then accelerate to the same maximum velocity in the opposite direction. These tidal currents are known as the ebb and flood and will be correlated to the peak height of the tides. The shape of the sea floor, referred to as the bathymetry, can affect the magnitudes of the tides. These currents are intensified when water is funnelled through basins, between straights, narrows, around peninsulas and islands. In addition, the tides amplitude will fluctuate over a lunar month influenced by the relative position of the moon, the sun and the earth known as the astral alignment. The highest and lowest “high tides” are called the spring and neap tides respectively. Therefore, the nature of the fluctuations the plants will determine the location and the power production of the tidal plants. These types of fluctuations are predictable with enough knowledge of the site conditions. 10 Tidal energy source The tides are mainly caused by the gravitation pull of the moon. The sun exerts a less important gravitation force pull on to the tides. Although the moon is smaller in mass, it exerts a stronger force onto the tide than the sun due to its proximity to the Earth. The gravitational force exerted on the largest bodies of water are describe as follows: (1) 𝐹𝑑𝑚 = 𝐺𝑀𝑚𝑅𝑒 𝑑𝑚 3 [9] 𝐹𝑑𝑚 = Tide generating force of the Moon [N] Mm = Mass of the Moon[kg] Re = Radius of the Earth G = Gravitational Constant d = Distance between the centre of mass of the two objects (2) 𝐹𝑑𝑠 = 𝐺𝑀𝑠𝑅𝑒 𝑑𝑠 3 [9] 𝐹𝑑𝑠 = Tide generating force of the Moon [N] Ms = Mass of the Sun [kg] The above equations (1) and (2) are derived from Newton’s Law of Gravity and Newton’s Second Law of motion. These simplified equations show the cubic inverse relationship between distance and the deferential force. The force is a function of the distance between the centres of masses of the two bodies. The tide-generating force of the Moon is 2.18 times the tide generating force of the Sun due to fact that the moon is so much closer than the sun. Therefore, the tides will rise and fall with relations to the distance of the arterial alignments. The accumulation of the tangential force exerted on the surfaces on the oceans that squeezes to form the bulge of the oceans. This force is known as the traction force. Therefore, the tides only occur in large body of water with sufficient continuous surface. Figure 3 Forces of gravity exerted by the moon onto the large bodies of waters If the Earth was covered by a single ocean, the effective force of gravity of the Moon is seen both perpendicular and in line with the Earth-Moon alignment, as shown Figure 3. It is important to understand that the forces aligned with the Earth-Moon line is not sufficient to cause the tidal bulge. It is the traction force form the entire 11 continuous ocean surface that squeezes the ocean surface to form the tidal bulge. The differential traction force is demonstrated in Figure 4. The effects of the spinning earth cause the bulge to appear on opposite sides of the planet. Figure 4: Accumulation of traction forces to create tidal bulges [9] The continents and the ocean bathymetric can obstruct or intensify the timing and amplitude the tides. The figure below show a map of the tidal heights around the world. Figure 5 Mean Daily Tidal Range [10] The regions of high tides are not uniformly distributed around the world. This is due to the influences of the continents which produce complex geographic variations with the ocean basins [11]. Daily tidal fluctuation The tides are constantly rising and falling and are characterized by the number of peaks and lows in a lunar day of 24 hours and 50 minutes [12] . Figure 6 shows the tides that peak once a lunar day, known as diurnal tides. Semi-diurnal tides peak twice a lunar day (Figure 7). There are also mixed tides that are a combination of diurnal and semi- diurnal tides (Figure 8). 12 The timing and magnitude of tides are locally effected by the geographical location, the seabed depth and shape of bays and estuaries [13]. Funnel-shaped coastlines will exaggerate tidal rise. The opposite is true for narrow inlets. Shallow sea floor dissipate tides due to bottom drag [13]. Narrow and long bays, such as fjords, can exhibit a special property where the tides can be exaggerated in part because of the resonance frequency of the bay. As the sea level rises in the open ocean from the tides, a “tidal wave” enters the bay and propagates until it is reflected back toward the open ocean. If the timing is right, under the right conditions, this reflecting wave can overlap and interfere with the next tidal wave coming from the ocean at the next tidal cycle, creating a standing wave. The resulting the tides can “slosh” around like a water in a moving bath tube. In order for this phenomena to occur the bay inlet (L) must be less than twice the Rossby radius (𝐿𝑅). The Rossby Radius is a function the sea floor depth (H) and the regional tidal cycle frequency (f) as defined in equation (3) [14]. (3) 𝐿𝑅 = √𝑔𝐻 𝑓 When the geographic conditions satisfy inequality 𝐿 < 2𝐿𝑅 , the area can experience very large tides and fluctuate such as in the Bay of Fundy on Canada’s East coast. Monthly tidal fluctuations The height of the high tides will also fluctuate over the course of a lunar month: 29 and half days. When both the moon and the sun are in line, the amplitudes of the tides are exaggerated and are referred to as the spring tides. Yet when the two celestial bodies are perpendicular, the tidal effects are negated, and so high tide does not reach as high. These are called the Neap tides. The orientation of the tidal bulges are show in Figure 9 , for the moon phases in a lunar month. Figure 6 Diurnal Tides [11] Figure 7 Semidiurnal Tides [11] Figure 8 Mixed Semidiurnal [11] 13 Figure 9 Astral Alignment for the Spring and Neap Tides In addition to Spring and Neap tides, there are also Perigean and Apogean tides that are due to the elliptical orbit of the moon that effect semi-diurnal tides. Tropical and Equatorial tides are the result of the change in declination of the Sun and Moon and effect diurnal tides. These tides introduce further oscillation to consider. Other sources of tidal fluctuations Meteorological tides are the result of variation of local sea level due persistent strong winds and barometric pressure changes. These “tides” can cause sea level to rise on the same scale of lunar tides (several meters) in combination with storm surges. High pressure weather systems can lower these tides. Similarly, low pressure weather system can raise high tides [9]. Meteorological tides are neither cyclical nor predictable and are not considered in this study for energy extraction. Coastal and island communities can take advantage of the high-energy areas to power their communities. Giving the right geographical conditions and bathymetry, these currents can be funnelled up to 5.04 m/s (highest in North America [15]). This funnel effect can be found between straights and around peninsulas. These special locations with high energy content are more favourable locations for tidal energy production. Therefore, tidal power is variable but highly predictable; 35 days is the recommended time for sufficient data set to predict tides [16]. Relationship between tide height and tide current The ebb and flood currents accompany the rise and fall of the tide. Tidal currents are described by the predominant tidal cycle that cause the tide to rise, i.e. a spring current brings in the spring tides, etc. The following graph demonstrate different correlation 14 patterns between the height (in black) and the tidal current strength and direction (in red). Slack water are periods of stand still as the current transitions from ebb to flood and vice versa. Figure 10 Relationship between tidal heights and tidal currents [15] These correlation between the tidal peaks and tidal current will vary on a spectrum form standing wave relationship to progressive wave relationship depending on the coastal topography and location. This emphasizes the fact that there exists a relationship between tidal height and the intensity and direction of the tides but the timing of the peaks may be offset by up to quarter a cycle [11]. 15 3 TUSK Plants Background A single Minesto Deep Green TUSK plant is rated for 500 kW. The first commercial scale plant is expected to be installed in 2017. As of February 2017, Minesto will be the first to install a multi- TUSK plant array. Up to 80 MW of tidal power generation capacity will be connected to the Welsh electricity grid [17]. This means the array will be up to 160 plants. Plant farm array A large array of TUSK plants can be connected together in smaller collections of plants. The plants must be at least two tether lengths apart from the nearest plant. Figure 11 is a top view of a possible array layout of seven TUSK plants arranged in a flower pattern. The grouping can be an odd or even numbers of plants. This will affect the coordination strategy for operating the plant relative to each other. Figure 11 Seven plant TUSK array, top view Synchronous start up During a lunar cycle, the tidal current velocity will reduce to a stop in order to switch direction from ebb to flood. During this pause, the neutrally buoyant plant will float in slack. Then, as the tidal current ramps up, the plant will natural extend downstream to its farther point. Once the current is sufficiently strong, the plant will be being producing power and start on its trajectory. Therefore, it is assumed that all the plants in a tidal farm will naturally be synchronized since they all experience this start-up process. 16 Figure 12 Synchronous start-up of a three plant tidal farm, side view This means that without external control, the power output of a TUSK tidal farm will introduce large fluctuations into the local power grid or require further power quality and storage investment. For normal operations, the shifting will occur during or shortly after start-up. This reduces the impact of shifting onto the farm’s power output as plants may need to be slowed down or take a longer route to be appropriately off set relative to other plants in the array. Current research on TUSK tidal farms Tethered Under-Sea Kit turbines are a new technology. Only recently has there been major focus on these types of systems for off shore power [2]. Major breakthroughs in related fields and growing popular interest and political incentives have contributed in bringing this technology to light. Published works on TUSK systems are limited [18]. The idea for submersible plant using tendered kits for harvesting of tidal energy was introduce by Minesto’s patent in 2012 [19]. Recent research has been conducted on Airborne Wind Energy (AWE) begging by a paper from Miles Lyod on “Crosswind Kite Power” in the Journal of Energy back in 1980 [20]. The methodology for modelling AWE can be transferred to TUSK system. In 2015 and 2016, David Olinger and Yang Wang from Worcester Polytechnic Institute published the first technical papers on modelling and simulation of tethered undersea kites from their experience on airborne wind energy [18]. This research focus on the plant dynamics to simulate the maximum theoretical power and issues of cavitation on turbine blades. The models developed predict that the plant power curves are not constant nor follow a simple sinusoidal shape, i.e. are not smooth repetitive oscillations. Marcus Jansson’s master thesis on the hydrodynamic simulations of the TUSK was published by Lund University in 2013. The report included a detailed assessment of the force balance and modelling software. The results were limited to a circular trajectory with a fixed rudder angle [21]. Moodley et al. conducted an economic analysis on the Minesto Deep Green TUSK technology compared to other marine technologies to develop the Agulas Current off the South African Coast in 2012. This study suggested that the most cost effective solution in areas with low velocity currents is the TUSK system [22]. The European PowerKit Project is a European Union’s Horizon 2020 funded project to gather open-sea experience to ensure high survivability, reliability and performance for TUSK systems. Starting in early 2016, this project has 9 industry and academic partners 17 and publish an environmental report on the TUSK system. Two of their main deliverables include a power take-off system and a complete array and grid connection for the TUSK plants [23]. Modelon Inc. is a modelling specialist and has also worked on modelling the TUSK system closely with Minesto using Dymola software. Modelon has published a blog post of their work with simulating TUSK systems in late 2016 with an example of the different angles of attack along the kite wing going throw a complete loop [24]. 18 4 Modelling the Plant Trajectory A mathematical description of the trajectory the plant movements in space is needed to being obtain a power performance curve. The formal definition of a function describes the relation between an input and exactly one output. To trace a path in kinematics, an independent input known as a “parameter” is the input and three outputs are needed to define a position in Cartesian coordinates: x, y and z coordinates. A solution is to use parametric equations. Parametric equations are a set of equations that together can output a position for an independent input. These equations are a convenient way to represent curves and surfaces [25]. Once such parametric curve is the Lemniscate of Bernoulli. The Lemniscate of Bernoulli is a useful description of a sideways figure 8 shape that the trajectory of the plant resembles. This pattern is preferred as the crisscrossing figure 8 shape allows for a tether to complete a loop without twisting. Bernoulli’s Lemniscate is describes as the following equations: In 2D Cartesian coordinates: (4) (𝑥2 + 𝑦2)2 = 2𝑎2(𝑥2 − 𝑦2) a = Amplitude In Polar coordinates: (5) 𝑟2 = 2𝑎2 Cos(2𝜃) r = Radius θ = Degrees Using parametric equations in 2D Cartesian co-ordinates: The parameter t will complete a full cycle at every 4π interval. (6) 𝑥(𝑡) = 𝑎 𝐶𝑜𝑠(𝑡) 1 + 𝑆𝑖𝑛(𝑡)2 t = parameter (7) 𝑦(𝑡) = 𝑎 𝐶𝑜𝑠(𝑡)𝑆𝑖𝑛𝑡(𝑡) 1 + 𝑆𝑖𝑛(𝑡)2 (8) 𝜅(𝑡) = 3√3 − 𝐶𝑜𝑠(2𝑡)√1 + 𝐶𝑜𝑠(2𝑡) 𝑎(−3 + 𝐶𝑜𝑠(2𝑡)) κ(t) = curvature 19 Figure 13 Bernoulli's Lemniscate (a=1) Other trajectories The path can also resemble a circle or a race track. The path configuration depends on many factors such as the material limits of the tether and the zone where the plant should operate with-in the water column. A circular path will have a constant curvature as seen in Figure 14 but will still be susceptible to current changes with the water column depth. Figure 14 Example of unit circular trajectory Both circular and the race track path requires the tether to have a ball joint to avoid twisting. These joints may complicate the electrical connections and reduce the tension tolerance. The optimal trajectory balances the minimal curvature and maximum straighten flight path. Turning is a form of acceleration which required an applied torque derived from lift power. Therefore, the sharper the path curves, the more power is extracted to turn and less is available to generate forward speed. A trajectory optimised with respect to tidal current distribution in the water column, tether tension or twist tolerance on power fluctuations is not included in this work. Bernoulli’s lemniscate described in three dimensions Geometric distortion is unavoidable when a two-dimensional figure is projected onto a spherical surface according to Carl Friedrich Gauss's Theorema Egregium [26]. Therefore, the characteristics of Bernoulli’s lemniscate cannot be exactly duplicated in three dimensions. A Viviani’s curve is used as an approximation in three dimensions. Viviani’s curve is defined as the intersection between a cylinder and a sphere; 20 specifically, when a cylinder shares a tangent line with the sphere. The Viviani’s curve is described by the following parametric equations as a function of the sphere radius (R) and the radius of the cylinder (r). (9) 𝑥(𝑡) = 𝑎 + 𝑟 ∗ cos(𝑡) r = Cylinder radius [m] R = Sphere radius [m] a = R-r = Distance between the cylinder axis and the sphere origin t = Parameter (10) 𝑦(𝑡) = 2 ∗ √𝑎 ∗ 𝑟 ∗ 𝑠𝑖𝑛 ( 𝑡 2 ) (11) 𝑧(𝑡) = 𝑟 ∗ sin (𝑡) (12) 𝜅(𝑡) = √13 + 3 ∗ 𝐶𝑜𝑠(𝑡) 𝑎(3 + 𝐶𝑜𝑠(2𝑡)) 3 2 κ(t) = Curvature Figure 15 Bernoulli's Lemniscate compared with Viviani's Curve For every parameter step “t”, there exit a position vector describe as �̅�(𝑡) = [ 𝑥(𝑡) 𝑦(𝑡) 𝑧(𝑡) ] in the Cartesian coordinate system. The point of origin is the centre of the sphere and is the anchor of tether. Meaning point [ 0 0 0 ] is the center of the sphere bounded by the tether length “R”. The path of the plant occurs at an elevation from the sea floor. Figure 16 Depiction of altitude angle 21 The altitude angle “D” is found between the sea floor and chord line from the origin to the midpoint (the criss-cross) of the lemniscate to define the rotation matrix: (13) 𝑅𝑜𝑡𝑎𝑡𝑖𝑜𝑛 𝑀𝑎𝑡𝑖𝑥 = [ cos(𝐷) 0 − sin(𝐷) 0 1 0 sin(𝐷) 0 cos(𝐷) ] The following diagram is the result of applying the rotation matrix to the parametric function to describe the Viviani’s curve trajectory. Figure 17 Depiction the intersection of a unit sphere and cylinder (left) and then applying the rotation matrix (right) Therefore, the trajectory of the plant can be described by using the length of the tether, the radius of the cylinder, and the degree of elevation into the parametric equation of the Viviani curve with a rotation matrix. The Figure 18 shows the rotated Viviani curve from different orthogonal perspectives. This simplification of the plant trajectory is enough to give insight into the potential power output. Figure 18 Orthogonal Projections of unit Viviani's Curve The Viviani’s curve consist of two mirrored pedals. The circumference of a pedal is given by the following formula: This is similar to the circumference of circle where instead of π, the circumference and radius ratio is 𝛤 = 7.64005 and use the radius of the cylinder 𝑟𝑐. Therefore, the full linear (14) 𝐶𝑝𝑒𝑑𝑎𝑙 = 𝛤 ∗ 𝑟𝑐 𝐶𝑝𝑒𝑑𝑎𝑙 = Circumference of one pedal [m] 𝑟𝑐 = Radius of the cylinder [m] 𝛤 = Viviani’s pedal constant 22 distance travelled by the plant to complete one cycle is twice the circumference of a pedal. Given a cylinder radius of 25 meters, the total Viviani’s curve perimeter is equal to: (15) 𝐶𝑉𝑖𝑣𝑖𝑎𝑛𝑖 = 2 ∗ 𝛤 ∗ 𝑟𝑐 = 2 ∗ 7.64005 ∗ 25.0 = 382.0 meters 𝐶𝑉𝑖𝑣𝑖𝑎𝑛𝑖 = Viviani’s curve circumference Plant Trajectory The simulated plant followed a Viviani’s curve with a tether length of 100 meters and a rotation of 45 degrees. The trajectory has a cylinder radius of a quarter of the tether length of 25 m. These values are approximations of real world TUSK plants for modelling. The Viviani curve is projected onto a sphere and cylinder in Figure 19 and orthogonal projections of the interpolation points are shown in Figure 20. Figure 19 3D simulation of base trajectory Figure 20 Orthogonal projections of Viviani curve The relative position of each interpolation point of this trajectory is the foundations for further simulation of plant velocity and power output. 23 5 Power Curve Modelling The objective is to model the power curve of a plant given only the trajectory and the tidal current intensity. The trajectory is described in three dimensions using the appropriate parametric equations. These equations are derived further to give direction and curvature at any point within the trajectory. A local coordinate systems is overlaid onto the trajectory. This enables to record how a curve segment manifests itself in changing the pitch, roll, and yaw of the plant. Knowing the rotation about these axis gives the required torque need to navigate the plant to stay on path. Power output is reduced to compensate for the torque applied. Figure 21 Methodology map from a given trajectory to power curve output Figure 21 describes the approach to find an approximate power curve of a TUSK plant from simply knowing the path geometry and the tidal current speed. The power output of a hydrokinetic turbines is governed by basic power equation (16). (16) 𝑃 = 1 2 𝜌𝐴𝐶𝑝𝑉𝑐 3 𝑃 = Power output [Watts] 𝜌 = Salt Water density [1025kg/m3] 𝐴 = Turbine inlet area [1.77 m2] 𝐶𝑝 = Coefficient of performance 𝑉𝑐 = Turbine linear velocity [m/s] Since the generator and turbine are placed on the kite, the plant operates as ‘drag power’. This means that as the kite gains cross current speed, the turbine can extract more power from the current, yet induces greater drag. The following derivation will demonstrate that there exist an optimal speed to extract power before incurring a detrimental amount of drag. Optimal power derivation The power produced is modelled in equation (17) and is the product of the available current power density, in brackets, the kite wing area, the lift coefficient of the kite and the relative drag power coefficient. The relative drag power coefficient represent the portion of the available lift power that produces turbine power. 24 (17) 𝑃𝑜 = ( 1 2 𝜌𝑉𝑤 3)𝐴𝑤𝐶𝑙𝑃𝑑 𝑃𝑜 = Power produced from drag model [W] Pd = Cross current kite relative drag power Aw = Wing reference area of kite [m2] Cl = Lift Coefficient of the kite 𝑉𝑤= Tidal current velocity [m/s] (18) 𝐿 = 1 2 𝜌𝐶𝑙𝐴𝑤𝑉𝑎 2 𝐿 = Lift force of the kite [N] The apparent linear or forward velocity can be related to the lift and drag of the kite from similar triangle ratios. Figure 22 shows the force balance of a weightless TUSK and the velocity triangle. The TUSK is directly downstream of the tidal current, so the tether tension 𝑇 is co-linear with the current velocity 𝑉𝑤 . The crosscurrent velocity 𝑉𝑐 is perpendicular to the current velocity 𝑉𝑤 . The kite drag 𝐷𝑡 and kite lift are also perpendicular. Since drag 𝐷𝑡 and the apparent velocity 𝑉𝑎 are collinear, then the velocities and the forces create similar triangles. Figure 22 display the drag vector 𝐷𝑡 parallel with the apparent velocity vector 𝑉𝑎 . Figure 22 Forces and velocity triangle on a crosswind TUSK (similar to [20]) From similar triangles ratios yields equation (19) which can be rearranged for the crosscurrent velocity in equation (20). (19) 𝑉𝐶 𝑉𝑤 = 𝐿 𝐷𝑡 𝑉𝑐 = Kite velocity cross current [m/s] 𝑉𝑤 = Tidal Current velocity [m/s] (20) 𝑉𝑐 = 𝑉𝑤 ∗ 𝐿 𝐷𝑡 𝐿 = Lift of kite 𝐷𝑡 = Total Drag of kite The cross current kite velocity is approximately equal to the apparent velocity for large lift-to-drag ratio. (21) 𝑉𝑎 = 𝑉𝑤 ∗ 𝐿 𝐷𝑡 𝑉𝑎 = Apparent Kite velocity [m/s] 25 Having a turbine on board loads the kite with additional drag 𝐷𝑝. The power extracted results in this additional drag 𝐷𝑝 travelling at velocity 𝑉𝑎 assuming no turbine losses. (22) 𝐷𝑡 = 𝐷𝑝 + 𝐷𝑘 𝐷𝑡 = Total Drag 𝐷𝑝 = Turbin Drag 𝐷𝑘 = Drag of kite (23) 𝑃𝑜 = 𝐷𝑝 ∗ 𝑉𝑎 𝑃𝑜 = Drag power Therefore, substituting equation (22) into equation (21) results in (24). (24) 𝑉𝑎 = 𝑉𝑤 ∗ 𝐿 𝐷𝑝 + 𝐷𝑘 Combining equations (17), (18), (23) and (24) to give the following ‘relative drag power output’ Pd as a function of drag ratio 𝐷𝑝 𝐷𝑘 . (25) 𝑃𝑑 = ( 𝐿 𝐷𝐾 ) 2 ( 𝐷𝑝 𝐷𝑘 ) (1 + 𝐷𝑝 𝐷𝑘 ) 3 Normalized the relative drag power function and plot in Figure 23. (26) 𝑃𝑑 ( 𝐿 𝐷𝐾 ) 2 = ( 𝐷𝑝 𝐷𝑘 ) (1+ 𝐷𝑝 𝐷𝑘 ) 3 = 0.5 (1+0.5)3 = 0.5 3.375 = 0.1481 = 4 27 Figure 23 Optimal drag ratio for maximum power output The maximum drag power 𝑃𝑑 is at: 𝑃𝑑𝑚𝑎𝑥 = 4 27 ( 𝐿 𝐷𝑘 ) 2 and this occurs at a drag ratio of 𝐷𝑝 𝐷𝑘 = 0.5 as shown in Figure 23. So for optimal apparent velocity from equation (24): (27) 𝑉𝑎_𝑚𝑎𝑥 = 𝑉𝑤 𝐿 (2𝐷𝑘) + 𝐷𝑘 = 2 3 𝐿 𝐷𝑘 𝑉𝑤 26 Therefore, equation (27) describes the maximum forward velocity as a function of the lift and drag ratio and the current velocity when the kite is at optimal drag power extraction. Power calculations The power equation used for this model is as followed: (28) 𝑃 = 1 2 𝜌𝐴𝐶𝑝𝑉𝑎 3 − 𝐿𝑖𝛾�̇� P = Power output [Watts] ρ = Salt Water density [1025kg/m3] Cp = Coefficient of performance Va = Apparent turbine linear velocity [m/s] Li = Torque coefficient [Nm/degree] 𝛾�̇� = Angular velocities in Pitch, Yaw and Roll directions [degree/s] The Coefficient of performance for the turbine is set to the optimal Betz limits for un- ducted turbine. This formula reveals that the power output will fluctuate with respect to 2 terms: the velocity term Va which is a function of the angle of attack and the torque term 𝐿𝑖𝛾�̇� which reduces the power output when making sharper turns. Therefore, both terms in this power equation demonstrate the influence of the trajectory’s curvature and plant placement and orientation. Definition of world frame of reference, The World frame of reference (World F.o.R) is described by the simple Cartesian coordinate systems:[ 𝑥 𝑦 𝑧 ]. The centre point is defined by the plant anchor on the sea floor. The x-axis is in line with the tidal current flow and the z-axis indicate the elevation from the sea floor. Finally, the y-axis tracks the sided to side motion relative to the anchor. The parametric equations describing the plants trajectory is modelled from this frame of reference. Therefore, with a tether length of R, the initial point for an elevated path is:[ 𝑅 ∗ cos (𝐷) 0 𝑅 ∗ sin (𝐷) ]. 27 Figure 24 World Frame of Reference Definition of plant body frame of reference, The Plant frame of reference (Plant F.o.R.) is the local coordinate system is applied to each point interpolated along the path. Therefore the axis directions are changing relative to the World frame of reference, but seems fixed from the plant’s perspective. The origin of the body frame is the centre the plant. The three axis are defined by the plant position relative to the world axis’s origin. Figure 25 Plant F.o.R. with respect to World F.o.R. The three-dimensional orthogonal axes are e1, e2, and e3. Axis e3 is in line with the tether vector direction. The tether vector is found between the sphere origin and the plant position. Axis e3 extends beyond the plant pointing away (positive) from the world origin. Figure 26 Viviani's curve with interpolation points 28 The plant’s apparent velocity vector which points from the plant’s position towards the next interpolation point as shown on the Viviani Curve in Figure 26. The actual velocity vector is tangential to the curve and aligned with the chord connecting two adjacent points along the path. The e1 axis is defined as the velocity vector projected upon a plane that is tangential to the tether sphere. Thus, the e1 axis is always pointing forwards for the plant. Figure 27 is an example a local coordinate system on a two dimension for three interpolation points on a trajectory function F (t). The e2 axis is defined as the cross product between the e1 and e3 axis. This results in the e2 axis being co-linear with the plant wings. Figure 27 Example of a local coordinate system in 2-D Rate of change of yaw, roll and pitch The rate of change in Yaw, Roll and Pitch are related to the torque applied to turn the plant to re-orient it to travel toward the next interpolation point. This is the critical advantage of following a predetermined trajectory. The change in pitch angle is calculated as the velocity angle relative to the e1 axis on the e1 and e3 plane. Meaning a change of pitch of 0°, the plant is moving straight forward in the e1 direction. Since the plant is following the surface of a sphere and the e1 is co-linear with the tangent by definition, then the rate of change of the pitch angle should never be zero. The Yaw angle is measure from the velocity vector and the e1 axis on the e1 and e2 plane. The Roll angle is measured from the velocity vector and the e3 axis on the e2 and e3 plane. Figure 28 Depiction of changes in Pitch, Roll and Yaw 29 The change in pitch, yaw and roll, also known as the angular velocities, represent how the curves in the trajectory manifest itself from the plant’s perspective. Figure 29 presents the angles from the plant’s principle forward axis. Note, that the pitch angle is never zeros since the plant follows the curved surface of a sphere and will never goes in a straight line during operation. Figure 29 Changes in pitch, yaw and roll angles for one cycle Torque modelling The rate of change in pitch, yaw and roll reflects the how much the plant needs to turn to continue to follow the trajectory. Therefore, the sharper the turn, the larger the rate of change in angles and thus the greater the power need to apply the torques needed to make the turn. The relationship between the applied torque and the rate of change in angles is assumed to linear. Both computational fluid dynamic analysis and an intimate knowledge of the plant geometry would be needed to determine the non-linear torque functions. After consultation with Minesto, linear approximation of the torque behaviour is adequate. The torque coefficient Li is used to approximate the relationship. In addition, due to the hydrofoil shape, the torques needed to change the pitch, yaw and roll are not be equal and would differ on orders of magnitude. Therefore, the torque coefficient assigned to each angular velocity are the following: (29) 𝐿𝑖�̇� = [ 10𝑎 5𝑎 100𝑎 ] ∙ 2 ∗ [∆𝑃𝑖𝑡𝑐ℎ ∆𝑌𝑎𝑤 ∆𝑅𝑜𝑙𝑙] Let “a” be the scale appropriate for the plant size and operating environment. The following section is a series of power curves that progress in intricacy to result in a base simulation to describe a typical power curve for a TUSK plant. Each iteration adds another degree of complexion to the power curve until a satisfactory approximation is reached. Modelling the tidal stream The tidal steam is modelled as a uniform vector field. This means that at every point on the trajectory, the plant experiences the same stream magnitude and direction relative to the world frame of reference. 30 Figure 30 Vector field representing the tidal current stream In Figure 30, the tidal stream is represented as the vector field of equal magnitude and direction in blue. The trajectory is outline in black onto a hemisphere bounded by the tether length. This velocity profile does not consider “no-split condition” of the sea floor or surface wave interference. Angle of attack The angle of attack α is determined by the vector sum of the plant velocity 𝑉𝑐⃗⃗⃗ and the tidal stream direction 𝑉𝑤⃗⃗ ⃗⃗ ⃗. The angle of attack can also be defined as the angle between the where the wing is going and where it seems to be going. The resultant vector is the apparent velocity of the plant 𝑉𝑎⃗⃗ ⃗. The angle of attack is the angle between the apparent velocity and the cord length of the wing. The cord length can be appropriately pitched from the relative forward velocity to insure feasible angles of attack by adding an angle of inclination. In Figure 31, the chord line is co-linear with the actual velocity of the plant 𝑉𝑐⃗⃗⃗ . Figure 31 Angle of Attack α of the plant Lift and drag coefficient The plant is modelled as a single hydro-foil: NACA 0015. This hydro-foil was chosen because extensive studies has been done on this particular shape for a large range of angle of attack. According to Olinger and Wang [18], preliminary studies showed that the angle of attack for TUSK can exceed 20 degrees at low velocities and during sharp turns. In other TUSK modelling work, Marcus Jansson used NACA 012 [21]. Olinger also suggested that TUSK plant would more likely have a cambered foil, i.e. an 31 asymmetrical profile [27]. Later works of Olinger uses a NACA 0021 for numerical simulation of a TUSK plant in lift power production mode [28]. 5.9.1 Lift coefficient The 2-D lift co-efficient for NACA 0015 was determined by for a range of angle of attack from -90 degrees to 90 degrees [29]. Figure 32 Angle of Attack for NACA 0015 [29] Olinger and Wand used the following equations (30) and (31) to fit the curve presented in Figure 32. For an angle of attack (𝛼) between -20 degrees and 20 degrees [18]: (30) 𝐶𝐿 = −2.27 ∗ 10−4𝛼3 − 1.65 ∗ 10−19𝛼2 + 0.123𝛼 + 0.2 𝛼 = angle of attack [deg] For an angle of attack below -20 and beyond 20 degrees [18]: (31) 𝐶𝐿 = 5.16 ∗ 105 + 7.3 ∗ 10−24𝛼4 − 9.06 ∗ 10−6𝛼3 − 9.27 ∗ 10−20𝛼2 + 0.0405𝛼 + 0.2 Figure 33 Angle of attack curve fit using equations (30) and (31) [18] 5.9.2 Drag coefficient: The total wing drag co-efficient (𝐶𝐷𝑊) can be determined as a function of the lift co- efficient using the sum of parasitic drag and induced drag [18]. (32) 𝐶𝐷,0 = Parasitic drag coefficient = 0.05 [18] 32 𝐶𝐷𝑊 = 𝐶𝐷,0 + 𝐶𝑙𝑊 2 𝜋𝑒𝐴𝑅 𝐶𝑙𝑊 = Kite Wing total Lift Coefficient e = Oswald efficiency = 0.9 [18] AR = Aspect Ratio = 3 [18] The following Figure 34 shows the lift coefficient for a wide range of angle of attack from -90 degrees to 90 degrees. The drag coefficient is also displayed on the same axis as a function of the lift coefficient. It is to be noted that the drag coefficient is not valid after 40 degrees. Beyond this point, the drag coefficient is no longer a direct function of the lift coefficient and equation (32) is no longer a satisfactory approximation. Figure 34 Relation between lift and drag co-efficient for a wide range of angles of attack Affective velocity For this model, the plant is assumed to be able to operate at the optimal affective velocity for a given trajectory and current velocity. Per Loyd [20], the optimal cross current velocity for a tether system is as followed: (27) 𝑉𝑐 = 2 3 𝐶𝐿 𝐶𝐷 𝑉𝑊 𝑉𝑐 = Cross current Velocity [m/s] 𝑉𝑊 = Current Velocity [m/s] 𝐶𝐿 = Kite Lift Coefficient 𝐶𝐷 = Total Drag Coefficient See theory section for derivation of equation (27) from Loyd [20]. Therefore, the affective velocity is a function of the lift to drag ratio which is determined by the angle of attack. Angle of inclination Combing the lift (30) and drag (31) coefficient equations into the affective velocity equation (27) to give equation (33). The affective velocity is sensitive to the angle of attack. Equation (33) is for angles of attack α between 0 and 20 degrees: 33 (33) 𝑉𝑐 = 2 3 ∗ 𝑉𝑤 ∗ ( (−2.27 ∗ 10−4𝛼3 − 1.65 ∗ 10−19𝛼2 + 0.123𝛼 + 0.2) 𝐶𝐷,0 + (−2.27 ∗ 10−4𝛼3 − 1.65 ∗ 10−19𝛼2 + 0.123𝛼 + 0.2)2 𝜋𝑒𝐴𝑅 ) Therefore, an angle of inclination can be imposed to correct the angle of attack to the optimum value. This is done by pitching the wing relative to the plant. Pitching the wing also requires an applied torque. This is included in the torque calculations in section 5.17. Angle of attack and angle of inclination The angle of attack is calculated from the vector sum from the tidal current relative to the plant’s position and the affective velocity of the plant. The plant needs to operate at the optimal turbine drag and kite drag ratio and at the optimal lift to drag ratio to reach maximum affective velocity. Figure 35 Angle of attack for one cycle Without an angle of inclination, the angle of attack during operation varies between 8° and 14.5°. This range of angle of attack is found to be at a local minimum of the lift to drag ratio, highlighted in yellow in Figure 36. Yet, the local maximum for lift to drag ratio occurs at an angle of attack of 3.7°. Therefore, an angle of inclination is necessary to ensure that the angle of attack is reduced to optimal operations. Figure 36 Lift to Drag Ratio plotted with comparison to lift efficient plot and drag coefficient plot. The angle of attack is calculated for a single cross section of the TUSK wing. The angle of attack will vary along the wing during a cycle. Modelon also modelled the angle of 34 a TUSK plant shown in Figure 37. Modelon calculated the range of angle of attack at various points along a TUSK’s wing with a 12 meter wing span following a similar trajectory for 1.6 m/s tidal current. This figure was generated using Modelica and Dymola modelling software. This shows that calculated angle of attack in this study is approximately within the same magnitude and range of the angle of attack calculated by Modelon. Figure 37 Modelon model of angle of attack and how it varies along the wing [24] The angles of attack in Figure 37 are plotted for one side of a TUSK wing. Meaning that for half the cycle, the wing half is on the outside of the trajectory and the wing tip experiences higher velocities relative to the centre. The other half of the cycle, the wing half is inside the trajectory and experiences lower velocities relation to the centre. This is why the angle of attack is asymmetric in Figure 38 (left). Figure 38 Variation of the TUSK angle of attack over the trajectory In Figure 38 is the angle of attack of the TUSK plant modelled in this study. The angle of attack is assumed to be averaged over both sides of the TUSK wing. This is why the angle of attack is symmetric in Figure 38 (left). Also, the modelled trajectory has a wider breath than the Modelon trajectory. This shows that the angle of attack will vary for slightly different trajectories. Fixed angle of attack and variable angle of inclination To compensate for the variation in angle of attacked, the wing can be pitch relative to the plant by an angle of inclination presented in Figure 39. The new vector sum gives a steady angle of attack at the optimal angle of attack 3.7° for each position. Yet, the consequence is that the plant will still have to exert a torque to pitch the wing and will experience a decrease in velocity show in section 5.6. 35 Figure 39 Adjusted angle of attack with adapted angle of inclination With the wing appropriately pitch at each position, the affective velocity is presented in Figure 40. Since the optimum lift to drag ratio is set to 6.512, according to equation (34), the maximum cross current velocity is 4.3 times the current tidal velocity of 1.6 m/s by the following proof. (34) 𝑉𝑐 = 2 3 𝐿 𝐷𝑘 𝑉𝑤 = 2 3 (6.512)𝑉𝑤 = 4.341 𝑉𝑤 Substitute in the tidal current of 1.6m/s in equation (34): 4.341 (1.6 𝑚 𝑠 ) = 6.946 𝑚 𝑠 Figure 40 Affective Velocity for one cycle Therefore, the affective velocity revolves around the predicted maximum affective velocity. The fluctuations seen in Figure 40 is simply due to whether the plant is going against the current or with the tidal current during the cycle. This influence is investigating in greater detail the next section 5.14. Relative velocity to the stream current In the Figure 40, the affective velocity’s mean is below the predicted maximum affective velocity. This is due to the relationship between the plant’s orientation and the direction of the tidal current velocity. 36 Figure 41 Mapping of plant orientation onto the trajectory The due to the geometry of the trajectory, the plant does not spend the same amount to time going against (in yellow) and going with the flow (in green). The plant spends more time going against the tidal stream as seen in yellow in Figure 41. When the plant does along with the flow the drop in relative velocity is more substantial. Therefore, the oscillation in velocity will centre on a lower mean due to this more significant dips. Variable angle of attack and fixed angle of inclination Another strategy is to fix the angle of inclination in order to reduce the range of angle of attack to include the local maximum of lift to drag ratio. This strategic eliminates the need to constantly apply a torque to change the wing pitch, yet the affective velocity is not as consistent as the previous strategy. Figure 42 Adjusted angle of attack with fixed angle of inclination The angle of inclination was set to 5.0°. This result in an range of angel of attack between 3.7° and 9.53° degrees, highlighted in yellow in Figure 43. Since the relationship between the angle of attack and the power curve output is not linear, a clear optimal of a fixed angle of inclination is not evident and would require further study. 37 Figure 43 Lift to drag ratio with the range of operation The resultant affective velocity is displayed on the Figure 44 for one cycle. The affective velocity is then calculated to include the relationship with the tidal current direction. Figure 44 Affective velocity for one cycle with a fixed angle of inclination Turbine velocity comparison In Figure 45, both affective velocities relative to the tidal current are presented to compare the different curve characteristics. The plant velocity of the variable angle of inclination has a higher mean velocity and smaller amplitude range. Figure 45 Comparison of turbine velocity and means 38 Power loss from pitch with angle of inclination torque For the variable angle of inclination mode, the affective pitch is the sum of the angle of inclination and trajectory pitch. Figure 46 shows the combination of the two pitch requirements to result in the effective pitch. Figure 46 Effective pitch angular velocity From the power equation (28), the torque losses are a function of the angular velocity. Therefore, Figure 47 shows the power curve before considering torques and after to highlight the torques influence on the power curve. In addition, the power curve that includes the torque require to pitch the wings to the correct angle of inclination is also include for comparison. Figure 47 Power Curves progression in complexity The asterisks indicate the interpolation points to relate to the position on the trajectory for reference. As seen in Figure 47, the effects of torques reduces the overall power output and introduces irregularities. Even when considering the angle of inclination, the influence of the relative velocity on the angle of attack dominates the power curve fluctuations. Final power curve In conclusion, the final power curve is presented in Figure 48. The power curve is the result of the plant following the Viviani curve anchored by a 100 meter tether. It includes the influence of a 1.6 m/s tidal current, the cross-current speed, the plant’s 39 position and orientation and the torque dynamics. This power curve will be used as the base model for further analysis. A sensitivity analysis of the parameters onto the power curve’s shape is not included in this study. Figure 48 Base power curve for one cycle 40 6 Mitigating Power Fluctuations As discussed in the background, the tides will ebb and flow at least once a lunar day the time scale of hours at various magnitude around the world. A tidal plant power production will fluctuate according to tidal partners. This has shown, a TUSK tidal plant will suffer an additional layer of power fluctuation during operation. The total distance the plant travels for a cycle is 382.0 meters for a tether length of 100 meter. The apparent velocity of each interpolation point is used to approximate the average speed of each chord between the interpolation points. Therefore, by aggregating the quotients of the average speed of the chords and its respective length; the total time to complete a cycle was determined to be approximately 57 seconds. The time scale for the fluctuation is in seconds to minutes. Meaning within a minute, the power output will peak twice. This puts the fluctuations on a time scale significantly less than tidal flow (hours) yet an order of magnitude more than grid regulation issues (seconds). This time scale of fluctuations requires a solution other than common power quality control measures such as short term storage, ramping capacity or spinning reserves. Therefore, this work will investigate a control strategy to eliminate or reduced the impact of TUSK plant operations fluctuations. Fourier analysis The power curve produced by the model plant for one cycle can be considered as an oscillating signal. A Fourier analysis is a useful tool to study the properties of the power curve signal. A complex signal can be broken down into multiple simpler sinusoidal signals of different amplitude, frequencies and phases by using the following Fourier Transform: (35) 𝑓(𝜔) = ∫ 𝑓(𝑡)𝑒−2𝜋𝑖𝑡𝜔 ∞ −∞ 𝑑𝑥 The example power curve function 𝑓(𝑡) in the time domain is be transformed into the frequency domain to produce 𝑓(𝜔). This allows for the spectral energy to be identified in a signal. A Fourier analysis reveals the dominant simple sinusoid present in the power curve as shown in Figure 49. These signals are identified by their magnitude relative to the other signals. These constituent signals are known as the signal harmonics. These simple signals can be recombined using Fourier synthesis to reproduce the original complex signal. This method can be used to simulate an identical power curve of a second plant without having to re-simulate from first principals, see all of section 4 and section 5. 41 Figure 49 Fourier transform from time domain to frequency domain (Modified from [30]) Fluctuation index The quality of the power performance is the amount of fluctuations in the power curve. A fluctuation index is used to rate the power performance quality. Figure 50 Example power curves and properties Figure 50 is an example of two power curves for a cycle duration. The multiple turbine signal is the sum of two turbines and a mean twice that of a single turbine. This “Fluctuation index” is the measure of the magnitude of the spread relative to the mean power output in percent. (35) 𝐹𝑙𝑢𝑐𝑡𝑢𝑎𝑡𝑖𝑜𝑛 𝐼𝑛𝑑𝑒𝑥 % = 𝑆𝑡𝑎𝑛𝑑𝑎𝑟𝑑 𝐷𝑒𝑣𝑖𝑎𝑡𝑖𝑜𝑛 𝐴𝑟𝑖𝑡ℎ𝑚𝑒𝑡𝑖𝑐 𝑀𝑒𝑎𝑛 ∗ 100 This equation is also known as the “Relative Standard Deviation”. A lower fluctuation index, the less significant the fluctuation are and the better the power performance quality. Signal shifting Combining two identical power curves in phase reinforce one another producing a new signal with bigger peaks and valleys. Yet, if the signals are shifted out of phase so that the peaks correspond with the troughs of the counterpart, the signals will “cancel” each other out resulting in a new signal with reduced variations and a lower fluctuation index. 42 Fourier analysis method Fourier analysis of the complex signal reveals the phases and magnitude of the constituent signal. The larger the magnitude of the constituent signals, the more influence it has on the original signal’s shape. The phases of the dominant signals are recovered. So, a new signal can be reconstructed using this information and the phase shifted by half the phase of the dominant constituent signals to compliment the original signal. This ensures that the major peak of the new signal corresponds with the most major through of the original signal. This method may not eliminate all fluctuations in the power curve but can be used as a heuristic to avoid large computation of the iterative method for multiple plants. Iterative method The Iterative method conducts a Fourier analysis on the original signal and reconstruct a new identical signal. The new signal’s phase is then sifted slightly and the two signals are added together to create the “sum signal”. The sum signal’s fluctuation index is evaluated and the phase shift recorded. The 2nd signal is then shifted slightly again with a new sum signal’s fluctuation index. This iteration continues until the 2nd signal has been shifted by one half of its period. All the fluctuation indexes are compared. The phase shift of the sum signal with the lowest fluctuation index is the optimal phase shift. This method ensures that the optimum phase shift is chosen but requires more computation, especially for multiple plant combination. Harmonics analysis of basic curve Figure 52 is the result of the Fourier analysis of the base power curve plotting in Figure 51. This reveals the dominant sinusoid present in the power curve. The zeroth bin’s magnitude represent the mean of the power curve at 145.3 kW. It also reveals that 2nd harmonic and 4th harmonics dominate the signal, and yet there is still a noteworthy 1st and 3rd harmonic. The fluctuation index reveals the influence of the largest harmonics. The Fluctuation index for the base power is 39.45%. Note that in Figure 51 the power curve is still for one complete cycle, which is 2π in radians for spectral analysis. Figure 51 Power curve of a complete cycle Figure 52 Harmonic analysis of base power curve 43 Two plants optimization When two plants are paired, their power productions are off set relative to each other. One plant it kept in phase, referred to as the master. The second plant, referred to as the slave, is offset or phase shifted relative to the phase of the master plant. This allows for the combined output of the paired plants to cancel out some of the harmonics to reduce the relative fluctuations. The Figure 53 has the power curve of a synchronized pair, where the two plants are in phase. The fluctuation index of the synchronized pair is 39.45%. This fluctuation index is expected to be the same as a single plant since the index accounts for the spread of the variance relative to the power mean. Figure 53 Power Curves of a synchronized pair compared with the optimal pairing Using the iterative method, the optimal time shift is applied to the slave plant. The result is also shown as the Optimum output in Figure 53 to compare with the synchronized pair output. Figure 54 Optimal pairing power curve The optimal shifted pairs power output has a fluctuation index of 7.60%. This an 80% reduction in fluctuation. Figure 55 shows the fluctuation index of each iteration of the phase shifting. Note that the local minima is almost at ¼ cycle shift. 44 Figure 55 Paired plants power curve analysis for each iteration The benefit of pairing is revealed by conducting another Fourier analysis on the optimized power curve. As shown in Figure 56, only the fourth harmonic remains as the significant outlier. Figure 56 Frequency energy of the optimized pairing power curve Three plants optimization When combing three plants to smooth out the fluctuation, one plant is considered the master and is kept in line with its original phase and the other two plants are shifted with respect to the master plants, referred to as slaves. Figure 57 shows the power curve of the three plants in phase with each other and the optimized power curve, where the three plants are appropriately shifted. One master plant and two slave plants are referred to as triplets. The fluctuation index is reduced to 2.13% for the triplets. Figure 57 Power curves of synchronized triplets and optimized triplets 45 Figure 58 Optimum triplets off set Figure 58 shows the optimum off set of both slaves in the triplet. The optimisation was done using the iterative method. The resultant power curve is displayed in Figure 59. As the dominant harmonics are cancelled out, the high order harmonics are seen to influence the power curve’ shape. Yet, their overall impact is less noticeable relative to the power output’s magnitude. This phenomenon is evident after a spectral analysis on the optimized power curve in Figure 60. Figure 59 Power curve of triplets at optimum off set Figure 60 Fourier analysis of the triplets optimized power curve (unlabled and labled) In Figure 60, the spectral analysis relevels that the power output mean is approximately 434.8 kW, yet the second harmonic is a 40th fraction of the magnitude or 2.6%. Two pairs Finally, two pairs of plants are simulated to further reduce the fluctuations. Each pair has been optimally shifted and now the pair is shifted relative to another optimally shifted pair. The power curve of a pair is shown in Figure 54. The Figure 61 is the 46 power curve of the plant pairs in synchronous operation and optimally shifted using the iterative method. Figure 61 Power curves of two plant pairs The first optima shift is found at a 1/8th of the cycle. This confirms that the dominant harmonic for each pair was the fourth harmonics. This was shown in Figure 56 using spectral analysis. This two pair shifting reduces the fluctuation index to 0.81%. Figure 62 Fluctuation index for shifting two plant pairs Figure 63 is the power curve of the four plants combined in this fashion. Again, the high order harmonics are more influential on the shape of the power curve but have less impact on the magnitude of the fluctuation as seen in Figure 63. Figure 63 Power curve of the optimum shift between two plant pairs Therefore, the more plants are involved to coordinate shifting, the more harmonics can be cancelled and reduce the overall fluctuations. Simply combining the pants in synchronous does not reduce the relative fluctuation, since the fluctuations scale with respect to the number of combined plant power output. There is a diminishing return on 47 plant shifting as each pairing cancels out harmonics, the higher order harmonics have less impact on fluctuation magnitude. The trend of diminishing returns is displayed in Figure 64. Figure 64 Fluctuation Index of each plant combination Therefore, the operation strategy of shifting the plants can by applied for odd and even collections of plants. Plant pairs can be used to greatly minimize the power fluctuations of an array and plant triplets can be used to include any remaining plants. Below are two example strategy for coupling plants in a seven plant array to improve power output. Figure 65 Tidal farm with multiple plant pairs and a plant triplet Figure 66 Tidal farm with 4 plants pairing and a plant triplet The pairing strategy in Figure 66 has the advantage of targeting 3rd and 4th order harmonics directly if these are significant enough to merit the greater control intervention, otherwise the pairing strategy in Figure 65 require less control intervention and still reduce power fluctuations. Benefits of reducing fluctuations for tidal energy This work has shown that the power output from tidal energy has inherit fluctuations. This fluctuation may cause supply and demand imbalances know as contingencies. Large integration of tidal power in an electricity network will change the supply conditions and requires contingencies management. Considerable effort is spent by 0 5 10 15 20 25 30 35 40 45 0 1 2 3 4 5 FL U C TU A TI O N IN D EX PLANTS Fluctuation Index 48 power system operators to anticipate contingencies to maintain system frequency [31]. The grid can be constraints by limited ramp rates and limited energy supply of other power producers. The operation strategy suggested in this work can be a preventive alternative for operators. This operation strategy also allows for better forecasting of tidal energy which leads to better market reliability and thus a higher priority in an optimal dispatch order. Future work This study approximate the power curve for a given trajectory of TUSK plants using Matlab modelling. Simplifying assumptions used to approximate parameters and reduce computations such as neglecting tether dynamics and uniform tidal velocity profile. A sensitivity analysis on input parameters can refined the modelling of the power output. A stochastic study of the plant power output can be used to investigate the likely hood of plant synchronization of multiple plants. Further research and field tests or historical data are needed to develop a historical probability distribution of power output of TUSK tidal plants. From this stochastic analysis, an optimal number of plants can be coupled to reduce invasive controls operations. 49 7 Conclusion In conclusion, tidal energy production is a new renewable energy technology and has the potential to locally change power supply condition of electrical grids due to it intermittency. Tidal energy’s main advantage comes from its predictability compared to other renewable energy sources. This study highlighted the benefits of the predictability of tidal energy with sufficient understanding of the natural cycles and modelling. The power cure of the TUSK can be reasonably modelled given only the geometry, the tidal velocity and using simplified assumptions on the wing shape and tidal velocity profile. The power curve of a TUSK plant has been show to fluctuate on the time scale of approximately a minute during operation. The apparent velocity of the plant relative to the tidal current was shown to have the most influence on the power curve’s fluctuation for during optimal operations. Fourier analysis allowed to investigate the impact of various scales of fluctuations with in the power curves. Using phase shifting, the plants can be couples in multiples of pairs or in triplets. This operation strategy significantly reduces power fluctuation without incurring addition cost for contingency management. 50 8 References [1] International Energy Agency, “Renewables Subtopic Ocean: IEA,” International Energy Agency, January 2017. [Online]. Available: https://www.iea.org/topics/renewables/subtopics/ocean/. [2] International Renewable Energy Agency (IRENA), “Ocean Energy Technology readines, patents, deployment status and outlook,” August 2014. [Online]. Available: http://www.irena.org/DocumentDownloads/Publications/IRENA_Ocean_Energy _report_2014.pdf. [3] S. Gillman, “Tidal energy poised to turn commercial,” The EU Research & Innovation Magazine, 07 December 2016. [4] Minesto AB, “Technology Development: Minesto,” Minesto, January 2017. [Online]. Available: http://minesto.com/technology-development/. [5] M. N. S. C. S. P. C. R. Moodley, “A technical and economic analysis of energy extraction from the agulhas current on the east coast of South Africa,” IEEE Power and Energy Society General Meeting, pp. 1-8, 2012. [6] Tidal Energy Today, “TidalEnergyToday.com,” Minesto, 4 Jaunuary 2017. [Online]. Available: http://tidalenergytoday.com/2017/01/04/minesto-spreads-its- wings/. [Accessed 23 January 2017]. [7] S. TAKAO, S. TSUYOSHI, A. ATSUSHI and T. KASHIWAGI, “Reduction of Power Fluctuation by Distributed Generation in Micro Grid,” Electrical Engineering in Japan, vol. 163, no. 2, p. 14–20, 2008. [8] S. E. Anthony Lewis, “Ocean Energy,” Interngovernmental Pannel on Climat Change, 2010. [9] N. Bowditch, “Tides and Tidal Currents,” in The American practical navigator : an epitome of navigation, vol. 9, Paradise Cay Publications, 2002, pp. 129-150. [10] US Federal Energy Management Program, “Ocean Energy Technology: Overview,” p. 16, July 2009. [11] National Oceanic and Atmospheric Administration, “NOAA Ocean Service Education,” US govenment, 25 March 2008. [Online]. Available: http://oceanservice.noaa.gov/education/kits/tides/media/supp_tide07a.html. [Accessed 14 01 2017]. [12] O. Carlson, L. Hammar, Z. Norwood and E. Nyholm, Harnessing Energy Flows: Technologies for Renewable Power Production, B. Sandén, Ed., Gothenburg: Chalmers University of Technology, 2014, pp. 32-45. 51 [13] R. H. Charlier and C. W. Finkl, Ocean Energy, Tide and Tidal Power, Heidelberg: Springer, 2009. [14] S. H. K. K. LuAnne Thompson, “Tide Dynamics,” University of Washington, Seattle, 2005. [15] National Ocean Service, “FAQ - Tidal Current Predictions and Data,” NOAA, 15 October 2013. [Online]. Available: https://tidesandcurrents.noaa.gov/faq4.html. [Accessed October 2015]. [16] A. MacGillivray, “Ocean Energy: State of the Art,” SI OCEAN, Intelligent Energy Europe, 2013. [17] Minesto AB, “Press Releases: Minesto to expand the commercial roll-out of the Deep Green technology,” 9 February 2017. [Online]. Available: http://minesto.com/news-media/minesto-expand-commercial-roll-out-deep- green-technology. [18] Y. Wang and D. J. Olinger, “Modeling and Simulation of Tethered Undersea Kites,” in ASME 2016 Power and Energy Conference, Charlotte, 2016. [19] M. Landberg, “Submersible plant US 8246293 B2”. United States of America Patent PCT/EP2007/050924, 21 August 2012. [20] M. L. Loyd, “Crosswind Power,” Journal of Energy, vol. 4, no. 3, pp. 106-111, 1980. [21] M. Jansson, “Hydrodynamic Analysis and Simulation of a Tidal Energy Converter,” Lund University, Lund, 2013. [22] M. S. C. R. Moodley, “A Technical and Economic Analysis of Energy Extraction from the Agulhas Current on the East Coast of South Africa,” in IEEE Power and Energy Society General Meeting , 2012. [23] PowerKite Project, “PowerKite Project,” 2016. [Online]. Available: http://powerkite-project.eu/. [24] J. Malmqvist, “Modelon: Blog: Cross tidal kites,” Modelon, December 2016. [Online]. Available: http://www.modelon.com/blog/articles/cross-tidal-kites/. [25] C. Stover and E. W. Weisstein, “Parametric Equations: Wolfram Math World,” MathWorld--A Wolfram Web Resource. , 01 2017. [Online]. Available: http://mathworld.wolfram.com/ParametricEquations.html. [26] C. Rousseau, “Communication and Represntation,” in Mathematics of Planet Earth: Mathematicians Reflect on How to Discover Organize, and Protect Our Planet, SIAM, 2015, p. 109. 52 [27] D. Olinger, Interviewee, Insight on TUSK paper. [Interview]. 7 October 2016. [28] A. Ghasemi, D. J. Olinger and G. Tryggvason, “A Nonlinear Computational Model of Tethered Underwater Kites for Power Generation,” Journal of Fluids Engineering, vol. 138, no. 121401, p. 2, 2016. [29] R. E. Sheldahl and P. C. Klimas, “Aerodynamic Characteristics of Seven Symmetrical Airfoil Sections Through 180-Degree Angle of Attack for Use in Aerodynamic Analysis of Vertical Axis Wind Turbines,” Albuquerque, Sandia National Laboratories, 1981, p. 86. [30] A. D. Elster, “Fourier Transform,” Questions and Anwers in MRI , 2017. [Online]. Available: http://mriquestions.com/fourier-transform-ft.html. [31] M. R. H. Darryl R. Biggar, “Handling Contingencies: Efficient Dispatch in the Very Short Run,” in The Economics of Electricity Markets, Wiley, 2014, p. 229.