Model and Simulation of a Floating Hybrid Wind and Current Turbines Integrated Generator System, Part I: Kinematics and Dynamics

Abstract

This initial publication is part of a series of publications that will appear soon, which pursue a final objective for the proposal of a fully integrated and controlled hybrid system composed of a floating wind turbine—type “OC3-Hywind”—and two marine current turbines with the aim of increasing the energy generated by the floating installation and, at the same time, use the set of turbines as actuators as part of an integral cooperative control system of the floating hybrid system to ensure the structural stability of the floating hybrid generator system (FHGS) in harsh weather conditions, which is a key issue in this type of floating systems. A specially designed tool to design, analyze, and control this type of FHGSs was developed using Matlab^®. In this tool, named Floating Hybrid Generator Systems Simulator (FHYGSYS), several tests were carried out on the structural stability of the system considering the interactive phase of the acting forces. Working in a programming environment like Matlab^® allows design freedom and the possibility of evaluating the system with different geometries, aerodynamic airfoils, and external meteorological conditions, and also including or eliminating certain elements, etc. This versatility will be helpful in future studies aimed at evaluating this system and maximizing the production of energy.

1. Introduction

Wind energy extraction techniques have evolved from the first wind turbines installed on land, allowing the increase in their size with the objective of generating more electricity. At the same time, the idea of using the potential energy from the sea to produce electricity is being consolidated in recent years. Additionally, different energy generation systems are currently under development to achieve the capacity of extracting energy from waves, tides, or marine currents [1,2].

Floating wind turbines (FWT) [3,4] are the alternative to offshore wind power generation for depths greater than 60 m. Operating at a distance from coasts involves using more important wind resources and reducing the visual and noise impact. Nevertheless, the floating nature of these devices involves stability requirements to account for the control system further than the power production.

The first prototype [5], Blue H, was installed in 2008 in a water depth of 113 m. Since then, researchers and firms have commissioned other prototypes.

Windfloat^® is a patented FWT developed by Principle Power Inc [6]. The Windplus consortium used this wind turbine (WT) to design the first floating wind farm in offshore Portugal [7]. The first turbine was installed in 2010, and, in July 2020, the floating wind farm was fully operational. WindFloat Atlantic is grid-connected to Portugal, generating up to 25 MW.

RWE Global is also working on FWT prototypes [8]. DemoSATH is a 2 MW turbine over a concrete structure [9]. The prototype has a single point of mooring that provides self-alignment to the current and wave direction. This project is installed on the north coast of Spain. TetraSpar is a 3.6 MW turbine over a tubular steel structure [10]. The prototype test site is located in Denmark, 10 km offshore with a depth of 200 m. Aqua Ventus is an 11 MW WT over a concrete semi-submersible structure [11]. This project is expected to be operative in 2024 in New England (USA).

A disrupting FWT system is being developed by X1wind [12]. This firm developed PivotBuoy^® [13], a system capable of self-orientating the floating turbine to maximize the generated power, thus reducing the weights, making FWT more competitive.

Real prototypes and scientific research [14,15] show that the control of FWTs faces challenges that cannot be overcome with conventional power control techniques. High structural loads or platform movement due to wave currents and tidal forces impose stability restrictions not considered in traditional WT control. This vital requirement has pointed to the need for advanced predictive control and condition monitoring techniques to preserve structural stability and integrity, as opposed to conventional power control techniques.

In a conventional wind turbine with structural support to the ground, the angular displacement of the tower is relatively small, even under challenging wind conditions. The axial force of the wind mainly causes bending moments in the tower. Under these conditions, the weight of the nacelle acts by compressing the tower, not bending it.

The dynamics of large horizontal axis wind turbines—with structural support to the ground—can be modeled using a five-degree-of-freedom model. The dominant modes [16] include:

• out-of-plane deflection of the blade flap rotor;
• in-plane deflection of the blade edge;
• fore-and-aft tower motions;
• powertrain roll and twist.

As suggested in Figure 1, the dynamics of the deformations associated with these degrees of freedom tend to be coupled.

For example, the tower fore–aft motion is strongly coupled to the blade flap motion, and the tower roll motion is strongly coupled to the blade edge and power train torsion.

In this context, large modern WTs allow the application of control techniques that make possible an independent adjustment of the pitch angles of each blade [17]. Individual pitch control extends the conventional objectives of pitch control to reducing fatigue loads, particularly by active damping tower oscillations.

In the case of complex floating systems with six degrees of freedom, the system behaves as a mass–spring–damper system affected by changing forces resulting from wind flows and hydrodynamic forces due to waves and marine currents. Under certain conditions of higher-than-normal wind speed, conventional pitch control techniques introduce negative damping in the movement of the floating tower. That causes an excitation of the natural frequency and may cause the floating structure to resonate by applying decreases in the wind opposition when varying the pitch angle of the blades to regulate the active power generated. This phenomenon was observed in tests carried out at the Ocean Basin Laboratory at Marintek in Trondheim [18].

In an FWT, the support platform moves freely, and the tower can experience angular displacements of several degrees. In this case, the weight of the nacelle is directly related to the bending of the tower. Of note in this phenomenon are the effects of its amplitude and frequency. In terms of frequency, due to constantly varying wind and wave loads, significant fatigue stresses occur in the structure.

These considerations lead to the subjection of special operating conditions worthy of a detailed analysis, which can only be performed with the help of simulation software tools. Specifically, the simulation tool FAST [19] is to be highlighted for this purpose.

Additionally, to perform such an analysis, it is necessary to have a coupled dynamic behavior model of the floating support base and the tower/nacelle. This model will help us to define which variables should be monitored in the context of a condition monitoring system of the FWT structural system.

Attention must be paid to the interactions between the mechanical effects due to inertia loads (rotor, nacelle, and tower) and the electrical effects (generator, control, and protection systems).

The controller’s goal should be to minimize turbine and platform motion while limiting mechanical wear on the generator and transmission. In this sense, some simulation tools, such as FAST [20], have system linearization tools that can be used to design controllers based on the linear-quadratic gain theory. Another option is the design of controllers based on the Lyapunov theory for the minimization of energy functions. The FAST simulation tool was specifically developed to carry out pre-study tests on the overloads that can occur, among others, on the blades and tower of the floating wind turbines. It allows the testing of different control strategies, such as Gain Scheduling PID, LQR with Collective Blade Pitch, and LQR with Individual Blade Pitch or H∞.

Studies were initially conducted with individual target controllers for rotor speed regulation using collective blade pitch [21] for different FWT systems [22]. In the last ten years, individual pitch control has been the main trend in this area [23]. In [24], an individual pitch control scheme was developed to deal with blade and pitch actuator faults in FWTs. It is shown that, with these faults, conventional pitch control techniques fail. Modern techniques [25] (such as sliding control) have also been applied to pitch control in an FWT, with promising results. The proposed controllers can accomplish better power regulation, reducing the platform pitch motion and the blade load.

From the point of view of control engineering, hybrid FWT and marine current turbines (MCT) foundations are promising generators that can achieve more excellent stability than FWT. These foundations are a hot topic in the research and study of control algorithms and marine power generation. A series of objectives is raised that aim to develop advanced control, monitoring, and diagnostic algorithms for the maximum use of marine generation systems to increase the performance and structural stability of the system and ensure its economic viability.

The traditional objective of control systems applied to hybrid generation systems has been to maximize the energy generated. In the case of generation systems resting on a floating offshore platform, other requirements, such as structural stability, may be equally or more important. Due to the difficulty of access, costly maintenance, and expensive commissioning involved in having generators offshore, ensuring the physical integrity of the generator is the priority.

The hybridization offers an opportunity to address the problem of controlling the structural stability of the floating system. Stability enhancement has been considered a major challenge since the early days of floating wind turbine design. With this objective, in this work, a specific solution is exposed, consisting of a floating hybrid system composed of a wind generation subsystem and a generation subsystem with two marine current turbines. This proposal allows the development of an integrated control system which simultaneously deals with the structural stability of the system and the optimization of the generation capacity.

A hybrid system capable of taking advantage of both wind energy and marine currents to produce electricity with the aim of maximizing the performance of a structure installed in the sea was designed in [26,27,28,29,30]. In [26], the first version of this type of hybrid system called HYWIKIM was presented. In [27], the first behavioral hypotheses of HYWIKIM were exposed. In [28], a real HYWIKIM prototype was presented, which was tested at the Real Club Náutico de Valencia. In [29], the results of the tests carried out with HYWIKIM were extended and presented. In [30], the first version of the mathematical model that allows the simulation of the behavior of this hybrid system concept was presented, using an “OC3-Hywind”-type floating wind turbine described in [31], to which two marine current turbines were coupled, as those described in [32].

The novelty of this proposal is that it introduces the concept of cooperative, integrated control of the two generation subsystems involved to counteract tendencies towards instability, which, if not avoided, reduce the useful life of the hybrid system. This issue is highlighted as especially important in [33]. This is an inherent problem with hybrid floating systems where, under certain circumstances, the applied pitch control can cause the floating system to resonate.

In order to effectively use tidal turbines together with wind turbines, the wind and current resources must coexist. In some cases, one—or both—of the natural resources can threaten the floating system’s stability. There are many places worldwide—generally in the straits—where wind resources [34] coexist with tidal currents capable of activating tidal turbines: wind speed above 10 m/s and tidal currents above 2 m/s. Examples of these kinds of locations (

Table 1. Current and wind speed in the straits.

Location	Maximum Current Speed	Depth *	Average Wind Speed
Banks Strait	2.3 m/s	26.8 m	>10 m/s
Strait of Gibraltar	2.5 m/s	40 m	15 m/s
Straits in Florida	1.5 m/s	20 m	6 m/s
Strait of Malacca	0.3 m/s	30 m	4 m/s
Dover Strait	2.5 m/s	18 m	9.5 m/s
Euripus Straits	2.8 m/s	5 m	7 m/s
Strait of Messina	3 m/s	20 m	7 m/s
Cook Strait	3.4 m/s	34 m	>10 m/s
Alas Strait	1.85 m/s	40 m	5 m/s
Bosporus Strait	1.9 m/s	31 m	6.5 m/s
Roosevelt Island	2.25 m/s	4.25 m	6.5 m/s
Ushant Island	3.5 m/s	20 m	9 m/s

* Depths correspond where data were collected or computed with the current profile.

) are the Banks Strait in Australia [35], the Strait of Gibraltar in Spain [36], the Straits in Florida in the United States [37], the Strait of Malacca in Malaysia [38], the Dover Strait in England [39], the Euripus Straits in Greece [40], the Strait of Messina in Italy [41], the Cook Strait in New Zealand [42], the Alas Strait in Indonesia [43], The Bosporus Strait in Turkey [44], the tidal strait near Roosevelt Island in United States [45], or the Ushant Island [46,47]. The current speed depth ranges from 25 to 40 m. Therefore, turbines must be placed at the correct depth in each case to take advantage of these currents. In the case of low-speed locations, a mechanical augmentation channel [27] can be used to achieve the proper speed.

Although the design of the floating hybrid system was made by placing the marine current turbines at a depth of 20 m, the simulator can be adapted to carry out simulations by placing the marine current turbines at different depths, depending on the characteristics of the simulation site.

Continuing with the work carried out in [26,27,28,29,30], the authors plan to publish a series of articles that broaden the research in this field. It is intended to expose, in an exhaustive way, the mathematical modeling of the floating hybrid system carried out, dividing the exhibition into two publications—this being the first of them. Using this model, the study of the behavior of the hybrid system will be addressed by applying different control strategies. The series of works will be closed with a final installment in which the hybrid system will be evaluated by changing the characteristics of the elements that compose it: length and number of blades of marine turbines, different aerodynamic profiles, different number of mooring lines, etc.

In this work, the first version of the FHYGSYS tool is introduced using the one-dimensional theory to compute the thrust of the turbines. This theory allows the obtainment of a good approximation to know which behavior the steady state response system will have. The operational capacity of the tool was validated by comparing the results with the certificated test of the OC3-Hywind calculated in FAST 8. This comparison is offered in Appendix E.

Why did we choose to develop our own tool like FHYGSYS? Mainly because, when the research line began, there was no simulation tool that would clearly allow the simulation of marine current turbines, offering the possibility of freely designing different control strategies. The main objective of the research is precisely that: to study the behavior of a hybrid floating system like the one described using different automatic control strategies.

On the other hand, this research has allowed the development of a deep understanding of the mechanics of this type of complex system. In addition, the developed modeling techniques can be applied not only to floating systems, but also to any mechanical system. This offers a relatively easy path for the macroscopic modeling of mechanical systems, for example, for modeling systems in automatic control engineering.

This paper is structured as follows: Section 2.1 and Section 2.2 introduce the initial assumptions for the development of the tool. Section 2.3 exposes the mathematical model, showing the methodology followed for both kinematic and dynamics calculations. Next, Section 2.4, describes the forces that were considered relevant for the development of the mathematical model, but leaving their exhaustive explanation to the second part of this publication. Section 3 explains the results obtained in the performed tests. Finally, Section 4 presents the discussion and future work for this line of research.

2. Materials and Methods

2.1. Floating System

There are different solutions to the problems involved in installing wind turbines in the sea—some of them are compared in [48]. The solution chosen in this work is the one used in phase IV in the “OC3” project [31,49]. There are multiple studies about this solution, and thus, it can be considered a standard that allows us to easily compare the results obtained. This configuration is called “OC3-Hywind” [31,49] and is composed of a 5MW wind turbine [50] mounted on a spar-buoy floating platform [31]. It has a mooring system consisting of three lines anchored to the seabed [31,51], which prevents the floating system from drifting.

Continuing with the work carried out in [26,27,28,29,30], a mathematical model of a hybrid system capable of taking advantage of both wind energy and marine currents to produce electricity is developed (Figure 2).

The hybrid system is equipped with two marine current turbines (Figure 3) based on the designs described in [32,52]. These turbines perform a double function:

• the generation of electricity taking advantage of the kinetic energy transmitted by marine currents: in this situation, the system operates in a cooperative mode as a wind turbine and a double marine generator;
• contribution to the stabilization of the floating structure in hard weather conditions: in this situation, the system works in an active stabilization mode using the set of existing turbines as a cooperative control for system stabilization.

2.2. Modeling Hypothesis

For the realization of the mathematical model of the system, several modeling hypotheses were considered.

The mathematical model does not consider the flexible behavior of the structural elements from the system. Therefore, the mechanical analysis is carried out considering a structure composed of a set of rigid solids that cannot deform. This permits the location of their position just by placing a point belonging to the system in the Euclidean space. In principle, this simplification does not represent a considerable error related to the movement of the floating platform [49] and implies that the internal behavior of the different bodies from the floating system is ignored.

The floating system is moored to the seabed through three flexible lines. Their characteristics are found in [31,51].

The geometric depiction of the model is conducted using orthogonal systems of Cartesian coordinates. Two coordinate systems were mainly considered: a fixed inertial one (X, Y, Z) and a mobile one (x, y, z) that move together with the rigid solids from the floating system. Unless otherwise indicated, when reference is made to specific coordinates in this work, they will always refer to the inertial coordinate system (X, Y, Z) [51]. In the initial position, the origin of the two-coordinate systems match, only separating when some external force acts on the floating structure.

To represent the position of the floating system, six degrees of freedom were used (Figure 4), with three translations—surge, sway, and heave—and three rotations on each of the axes of the coordinate system—roll, pitch, yaw [31,51].

The position of the inertial coordinate system (X, Y, Z) is as follows [51]: the X–Y plane is located at sea level, the Z axis is on the tower, and the floating platform axis of symmetry—when found in its initial position—in the opposite direction to the acceleration of gravity, as shown in Figure 4.

At each instant of time, the next position of the mobile coordinate system (x, y, z) is calculated as a function of the accelerations produced by the forces acting on the floating system. From each new position of the mobile coordinate system, the position of the rest of the elements is also obtained. In [51], the process of calculating the position of the mooring system and the rigid solids that constitute the floating structure is similar to the one described above.

It is assumed that the classical hydrodynamic problem can be linearly modeled [51]. This assumption allows us, in the first place, to simplify the mathematical model of the incident waves, which is not adequate when the sea is in extreme conditions [51]. Secondly, it splits the hydrodynamic issue into three simple problems [51]: radiation, diffraction, and hydrostatic problem. Finally, it allows the application of the superposition theorem.

Although it is considered that the system consists of a set of rigid solids, some of these can rotate with respect to the rest. This is the case of the wind turbine and the nacelle, which can rotate to be aligned with the direction of the main wind or that of the support and marine turbines, which can turn to be oriented with the incident current.

2.3. Mathematical Model

The hybrid system is mainly composed of four structural elements: a floating platform, a tower, a nacelle, and support for current turbines (Figure 2). The system is also equipped with one wind turbine consisting of three blades, two marine current turbines (with two blades each), and a mooring system.

Figure 4 shows the floating system in its initial position. The positions with respect to the inertial coordinate system of the points of the elements that are part of the system are mostly extracted from [31,32,50].

The dynamic features (mass, center of mass, moments, and products of inertia) of each element of the system were obtained using Solid Works^®, representing each of the elements from the geometry and the dynamic characteristics that were extracted from [31,32,50]. All these inertial data are found in Appendix A.

This is an important piece of information for the simulation of the behavior of any mechanical system; in this case, it will allow the modification of the mechanical characteristics of any of the bodies of the floating system, thus obtaining simulations in future works.

The mathematical model of the hybrid system is an adaptation of the approach used in robotics for the development of the mathematical model of an industrial robot arm [53]. Following this approach, the mathematical model may be divided in two fundamental parts:

• Kinematics allows us to place in the orthogonal space all the points of the solids from the floating system without considering the forces that act on it. From the initial position of the mobile coordinate system (x, y, z), the position of the other points of the floating system is calculated at each instance of time.
• Dynamics processes all the forces that act on the hybrid system and calculates them at the origin of the mobile coordinate system, which allows us to obtain the resultant forces and moments that act on this point [54,55,56]. From these data, the new linear and angular accelerations can be, respectively, calculated [57,58,59], which enables us to obtain the new position for the six degrees of freedom of the system.

The relationship between kinematics and dynamics is shown in Figure 5. The dynamics is calculated from the inertial data and the external forces acting on the floating system. The mathematical models of kinematics and dynamics are explained in this paper in Section 2.3.1 and Section 2.3.4 respectively. The inertial data of the floating system appears in Appendix A and its processing is explained in Section 2.3.2. A brief but complete exposition of the forces acting on the floating system is shown in Section 2.4, while an exhaustive explanation of these will be presented in the second part of this paper.

2.3.1. Kinematics

According to [53], there are different methods to represent positions and rotations, and the homogeneous transformation matrices stand out among them, which represent the position and the rotation together, entailing a reduced computational cost. The use of homogeneous transformation matrices implies the representation of positions in homogeneous coordinates. To represent the rotations, roll, pitch, and yaw angles were used; this set of rotations is more suitable when modeling ships or boats instead of using, for example, Euler angles [53,60,61].

To obtain the position of the floating system, the kinematic model obtains in each instant of time a homogeneous matrix: M_HT(t), which represents the necessary rotations and/or translations to place a point from its initial position to the one that occupies at each moment. This Matrix (1) is generally composed of three rotations and three translations, also represented by means of homogeneous transformation matrices. As the matrix product is not commutative, the order of operations must be considered, pre-multiplying the matrices that represent each action.

$$M_{HT}\left(t\right)=M_{HT}{}_{BODY}^{INERTIAL}=M_{T{q}_{123}}·M_T{}_{CoR}^{BODY}·M_{RZ{q}_6}·M_{RY{q}_5}·M_{RX{q}_4}·M_T{}_{BODY}^{CoR}$$

(1)

The position (in meters) and rotation (in radians) of the mobile coordinate system within the inertial system is represented by a vector q_i, where i represents each of the degrees of freedom of the floating system (1 = surge, 2 = sway, 3 = heave, 4 = roll, 5 = pitch, 6 = yaw). The mobile coordinate system is represented by the word body for brevity—because it is attached to the body of the floating system—and for consistency with the representation adopted in [61].

In Equation (1), the order of actions to perform is from right to left. Thus, if the center of rotation is the same for all three rotations, the translation matrices of Equation (1) will be written as in Equation (2), where CoR_k are the coordinates at the initial position of the corresponding center of rotation.

$$M_{HT}{}_{BODY}^{INERTIAL}=\left[\begin{array}{cccc}1& 0& 0& q_1\\ 0& 1& 0& q_2\\ 0& 0& 1& q_3\\ 0& 0& 0& 1\end{array}]·[\begin{array}{cccc}1& 0& 0& Co{R}_1\\ 0& 1& 0& {CoR}_2\\ 0& 0& 1& {CoR}_3\\ 0& 0& 0& 1\end{array}\left]·M_{RZ{q}_6}·M_{RY{q}_5}·M_{RX{q}_4}·\right[\begin{array}{cccc}1& 0& 0& -Co{R}_1\\ 0& 1& 0& -Co{R}_2\\ 0& 0& 1& -Co{R}_3\\ 0& 0& 0& 1\end{array}\right]$$

(2)

In the mathematical model, Equation (2) is used to rotate the points of the wind and marine current turbines, to rotate the blades by the precone angle, and to rotate the nacelle-hub-blade assembly of the wind turbine. It was also used in [30] to position the floating system, so the global center of mass was chosen as the center of rotation of the mobile system.

In this work, the origin of the mobile coordinate center (Figure 4) was chosen as the center of rotation, to compare the results more adequately in the model validation process with those obtained with FAST [51,62]. This choice allows us to simplify Equation (1) because, in its initial position, the origin of the mobile coordinate system coincides with the origin of the inertial coordinate system, which converts the translation matrices of CoR into identity matrices. So, when Simplifying Equation (1) to position the floating system at each instant of time, Equation (3) is used instead.

$$M_{HT}{}_{BODY}^{INERTIAL}=M_{T{q}_{123}}·M_{RZ{q}_6}·M_{RY{q}_5}·M_{RX{q}_4}$$

(3)

The translation matrix appeared in Equation (3) is the same as the one expressed in Equation (2). On the other hand, the rotation matrices that appear in Equation (1)–(3) are the classic rotation matrices around the X, Y, and Z axis, but expressed as homogeneous Matrix (4) [53,60,63].

$$M_{HT}{}_{BODY}^{INERTIAL}=M_{T{q}_{123}}·\left[\begin{array}{cccc}\cos q_6& -\sin q_6& 0& 0\\ \sin q_6& \cos q_6& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}]·[\begin{array}{cccc}\cos q_5& 0& \sin q_5& 0\\ 0& 1& 0& 0\\ -\sin q_5& 0& \cos q_5& 0\\ 0& 0& 0& 1\end{array}]·[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos q_4& -\sin q_4& 0\\ 0& \sin q_4& \cos q_4& 0\\ 0& 0& 0& 1\end{array}\right]$$

(4)

From Matrix (4), the position, in the inertial coordinate system, of the points of the floating system at each instance of time can be obtained Equation (5). This operation applies to all the points at all time steps, which implies a relatively high computational cost. To facilitate computation time, instead of multiplying each point separately in Equation (5), the largest number of points is put together in a 4 × n matrix—where n is the number of points—and then Equation (5) is applied.

$$\left[\begin{array}{c}{p}_{INERTIAL}{}_x \\ p_{INERTIAL}{}_y\\ p_{INERTIAL}{}_z\\ 1 \end{array}\right]=M_{HT}{}_{BODY}^{INERTIAL}·\left[\begin{array}{c}{p}_{BODY}{}_x\\ p_{BODY}{}_y\\ p_{BODY}{}_z\\ 1 \end{array}\right]$$

(5)

2.3.2. Inertial Data Processing

One of the reasons for designing this mathematical model is to be able to modify, add, or reduce any of the single bodies that make up the floating system. Each single body has inertial characteristics, such as mass, center of mass, and moments and products of inertia, all of which define its mechanical behavior [54,55,56,57,58,59]. This means that, to obtain the global inertial characteristics of the floating system, a method must be available to unify the individual inertial characteristics of each single body.

To apply the dynamics equations—this will be explained in the following sections—the inertial characteristics of the floating system must be expressed at the origin of the mobile coordinate system. Starting with the mass, as it is a scalar quantity, the total mass of the floating system (m_FS) is the sum of the masses (m_i) of each of the n single bodies Equation (6).

$$m_{FS}={\displaystyle \sum }_{i=1}^n{m}_i$$

(6)

From the masses (m_i) and centers of mass (CoM_i) of each of the n single bodies and the total mass of the floating system (m_FS), the center of mass (CoM) of the floating system is obtained Equation (7), [54,55,56]. This calculation is conducted by expressing the points in the mobile coordinate system, so the resulting center of mass is expressed in the same system. Each instance of time, the mathematical model obtains the new center of mass, which changes due to the rotation of the turbines of the floating system.

$$\left[\begin{array}{c}Co{M}_x\\ Co{M}_y\\ Co{M}_z\end{array}\right]=\left[\begin{array}{c}{x}_{CoM}\\ y_{CoM}\\ z_{CoM}\end{array}\right]=\frac{1}{m_{FS}}·{\displaystyle \sum }_{i=1}^n\left[\begin{array}{c}{m}_i·Co{M}_{xi}\\ m_i·Co{M}_{yi}\\ m_i·Co{M}_{zi}\end{array}\right]$$

(7)

The processing of the moments (Ixx, Iyy, Izz) and products (Ixy, Iyz, Izx) of inertia is conducted by expressing these as an inertia tensor matrix (M_IT): Equation (8) [57,58,59]. The inertia tensor data for each single body (I_{ij (SB)}) that make up the floating system appear in Appendix A.1, Appendix A.2, Appendix A.3, Appendix A.4, Appendix A.5, Appendix A.6, Appendix A.7 and Appendix A.8 of Appendix A.

$$M_{IT}=\left[\begin{array}{ccc}{I}_{xx \left(SB\right)}& -I_{xy \left(SB\right)}& -I_{zx \left(SB\right)}\\ -I_{xy \left(SB\right)}& I_{yy \left(SB\right)}& -I_{yz \left(SB\right)}\\ -I_{zx \left(SB\right)}& -I_{yz \left(SB\right)}& I_{zz \left(SB\right)}\end{array}\right]$$

(8)

To express the inertia tensor of a single body at the origin of the mobile coordinate system, two simple actions are applied: first, rotate the inertia tensor to align it with the mobile coordinate system, and second, move the inertia tensor from its initial position to the origin of the mobile coordinate system. Obviously, if the inertia tensor is initially aligned with the mobile coordinate system, it is only necessary to translate it.

• Rotation about the center of mass of the inertia tensor.

Considering the single body located in its final position with respect to the mobile coordinate system and having the inertia tensor obtained with respect to a specific point of the single body—usually the center of mass—and aligned with a coordinate system attached to the single body. To align the inertia tensor, all the rotations made by the body must be undone until reaching its definitive position.

These rotations can be undone one at a time by choosing the corresponding equation from Equations (9)–(11), depending on whether the rotation is about the X, Y, or Z axis. This operation yields three unit vectors (u_i, v_i, and w_i) rotated by an angle opposite to that between single body and mobile coordinate system.

$$\left[\begin{array}{ccc}{u}_1& v_1& w_1\\ u_2& v_2& w_2\\ u_3& v_3& w_3\end{array}\right]=M_{rot{X}_{\left(-\alpha \right)}}·\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(-\alpha \right)& -\sin \left(-\alpha \right)\\ 0& \sin \left(-\alpha \right)& \cos \left(-\alpha \right)\end{array}\right]·\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$$

(9)

$$\left[\begin{array}{ccc}{u}_1& v_1& w_1\\ u_2& v_2& w_2\\ u_3& v_3& w_3\end{array}\right]=M_{rot{Y}_{\left(-\beta \right)}}·\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]=\left[\begin{array}{ccc}\cos \left(-\beta \right)& 0& \sin \left(-\beta \right)\\ 0& 1& 0\\ -\sin \left(-\beta \right)& 0& \cos \left(-\beta \right)\end{array}\right]·\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$$

(10)

$$\left[\begin{array}{ccc}{u}_1& v_1& w_1\\ u_2& v_2& w_2\\ u_3& v_3& w_3\end{array}\right]=M_{rot{Z}_{\left(-\gamma \right)}}·\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]=\left[\begin{array}{ccc}\cos \left(-\gamma \right)& -\sin \left(-\gamma \right)& 0\\ \sin \left(-\gamma \right)& \cos \left(-\gamma \right)& 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$$

(11)

From the inertia tensor Matrix (8) and the three unit vectors obtained through Equations (9)–(11), applying Equation (12), the new inertia tensor (M_{IT (ROT)}) aligned with the mobile coordinate system is calculated: Equation (13). Equation (12) was deduced from the equations used to transform known moments and products of inertia about one coordinate system into others relative to a second coordinate system that has the same origin but is tilted relative to the first coordinate system [57,58,59]. Equation (12) is a compact matrix version of the equations used in [30] for the rotation of inertia tensors.

$$M_{IT \left(ROT\right)}=\left[\begin{array}{ccc} \left[\begin{array}{ccc}{u}_1& u_2& u_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]& -\left[\begin{array}{ccc}{u}_1& u_2& u_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]& -\left[\begin{array}{ccc}{w}_1& w_2& w_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]\\ -\left[\begin{array}{ccc}{u}_1& u_2& u_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]& \left[\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]& -\left[\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right]\\ -\left[\begin{array}{ccc}{w}_1& w_2& w_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]& -\left[\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right]& \left[\begin{array}{ccc}{w}_1& w_2& w_3\end{array}\right]·M_{IT}·\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right]\end{array}\right]$$

(12)

$$M_{IT \left(ROT\right)}=\left[\begin{array}{ccc}{I}_{xx \left(ROT\right)}& I_{xy \left(ROT\right)}& I_{zx \left(ROT\right)}\\ I_{xy \left(ROT\right)}& I_{yy \left(ROT\right)}& I_{yz \left(ROT\right)}\\ I_{zx \left(ROT\right)}& I_{yz \left(ROT\right)}& I_{zz \left(ROT\right)}\end{array}\right]$$

(13)

All of this process was validated by comparing the results with those obtained from the creation of different solids with Solid Works^®, rotated at different angles with respect to a given coordinate system.

Although the exposed methodology has been explained for the case of a single rotation, this method allows the composition of several rotations in the same calculation, simply by multiplying the matrices with the desired rotations successively. This increases the computation speed of the simulation. Appendix B shows an example in which successive rotations of the inertia tensor of a single body were made until it was aligned with the mobile coordinate system.

2.

Translation of the inertia tensor from one point to another.

The translation of the inertial tensor of a single body from its center of mass to the origin of the mobile coordinate system is a simpler process; this is achieved by applying the parallel axes theorem or Steiner’s theorem [54,55,56,57,58,59].

From the center of mass of the single body (initial point) and the origin of the mobile coordinate system (final point), the position vector of the inertia tensor is obtained: Equation (14).

$$\left[\begin{array}{c}itp{v}_x\\ itp{v}_y\\ itp{v}_z\end{array}\right]=\left[\begin{array}{c}f{p}_x\\ f{p}_y\\ f{p}_z\end{array}\right]-\left[\begin{array}{c}i{p}_x\\ i{p}_y\\ i{p}_z\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]-\left[\begin{array}{c}Co{M}_{ xi}\\ Co{M}_{ yi}\\ Co{M}_{ zi}\end{array}\right]$$

(14)

With the position vector calculated in Equation (14) and the inertia tensor known, applying the parallel axes theorem expressed in Equation (15), the new inertia tensor (M_{IT (TRANS)}) translated to the origin of the mobile coordinate system is obtained: Equation (16).

$$\begin{array}{c}{M}_{IT \left(TRANS\right)}\hfill \\ =\left[\begin{array}{ccc}{I}_{xx \left(ROT\right)}+m_i·\left({{itpv}_y}^2+{{itpv}_z}^2\right)& I_{xy \left(ROT\right)}+m_i·{itpv}_x·{itpv}_y& I_{zx \left(ROT\right)}+m_i·{itpv}_z·{itpv}_x\\ I_{xy \left(ROT\right)}+m_i·{itpv}_x·{itpv}_y& I_{yy \left(ROT\right)}+m_i·\left({{itpv}_x}^2+{{itpv}_z}^2\right)& I_{yz \left(ROT\right)}+m_i·{itpv}_y·{itpv}_z\\ I_{zx \left(ROT\right)}+m_i·{itpv}_z·{itpv}_x& I_{yz \left(ROT\right)}+m_i·{itpv}_y·{itpv}_z& I_{zz \left(ROT\right)}+m_i·\left({{itpv}_x}^2+{{itpv}_y}^2\right)\end{array}\right]\hfill \end{array}$$

(15)

$$M_{IT \left(TRANS\right)}=\left[\begin{array}{ccc}{I}_{xx \left(TRANS\right)}& I_{xy \left(TRANS\right)}& I_{zx \left(TRANS\right)}\\ I_{xy \left(TRANS\right)}& I_{yy \left(TRANS\right)}& I_{yz \left(TRANS\right)}\\ I_{zx \left(TRANS\right)}& I_{yz \left(TRANS\right)}& I_{zz \left(TRANS\right)}\end{array}\right]$$

(16)

After rotating and translating the inertia tensor of all the single bodies that make up the floating system, the moments and products of inertia evaluated at the origin of the mobile coordinate system and oriented with respect to it are obtained. The global inertia tensor of the floating system (M_{IT (FS)}) is conducted by adding each of these inertia tensors Equation (17). This result may be seen in

Table A1. Inertial properties of the floating system.

Element	Floating Hybrid System
Density (kg/m3)	-
Mass (kg)	8,138,259
Ixx (kg·m2)	67,434,701,761
Iyy (kg·m2)	67,398,679,463
Izz (kg·m2)	144,576,159
Ixy (kg·m2)	7.56248
Iyz (kg·m2)	−75.8696
Izx (kg·m2)	−12,699,829
Center of Mass
Surge (m)	−0.0172219
Sway (m)	0
Heave (m)	−76.6108

of Appendix A.

$$M_{IT \left(FS\right)}={\displaystyle \sum }_{i=1}^n{M}_{IT i \left(TRANS\right)}=\left[\begin{array}{ccc}{I}_{xx}& I_{xy}& I_{zx}\\ I_{xy}& I_{yy}& I_{yz}\\ I_{zx}& I_{yz}& I_{zz}\end{array}\right]$$

(17)

Finally, all the calculation carried out in the inertial data processing is collected into Matrix (18), a rigid body matrix (M_RB), thus facilitating the use of these data in the calculation of dynamics [61].

$$M_{RB}=\left[\begin{array}{cc}\begin{array}{c}{m}_{FS}\\ 0\\ 0\end{array}& \begin{array}{c}0\\ m_{FS}\\ 0\end{array}\\ \begin{array}{c}0\\ m_{FS}·z_{CoM}\\ -m_{FS}·y_{CoM}\end{array}& \begin{array}{c}-m_{FS}·z_{CoM}\\ 0\\ m_{FS}·x_{CoM}\end{array}\end{array} \begin{array}{cc}\begin{array}{c}0\\ 0\\ m_{FS}\end{array}& \begin{array}{c}0\\ -m_{FS}·z_{CoM}\\ m_{FS}·y_{CoM}\end{array}\\ \begin{array}{c}{m}_{FS}·y_{CoM}\\ -m_{FS}·x_{CoM}\\ 0\end{array}& \begin{array}{c}{I}_{xx}\\ -I_{xy}\\ -I_{zx}\end{array}\end{array} \begin{array}{cc}\begin{array}{c}{m}_{FS}·z_{CoM}\\ 0\\ -m_{FS}·x_{CoM}\end{array}& \begin{array}{c}-m_{FS}·y_{CoM}\\ m_{FS}·x_{CoM}\\ 0\end{array}\\ \begin{array}{c}-I_{xy}\\ I_{yy}\\ -I_{yz}\end{array}& \begin{array}{c}-I_{zx}\\ -I_{yz}\\ I_{zz}\end{array}\end{array}\right]$$

(18)

The moments and products of inertia included in Equation (8)–(18) depend on time, the notation (I_ii (t), I_ij (t)) has not been included in the expressions for clarity. The variation with time given in these values is very small and can even be neglected without appreciating significant changes in the behavior of the floating system.

2.3.3. Added Mass Processing

A body moving in a fluid behaves as if it has more mass than it really does [64]. The processing of the added mass allows this effect to be considered. According to [65], the added mass can be interpreted as a particular volume of fluid particles that are accelerated with the body.

The calculation of the added mass depends only on the body shape and the density of the fluid [64]; this implies that the added mass depends on the area of the floating system that is submerged at each instant of time. For this reason, the calculation process starts by deducting the submerged volume of the floating system Equation (19). The single body that will be partially submerged is the floating platform—the inertial characteristics of the body and the volume that is submerged when the floating platform is in the initial position appear in Appendix A.1. The support of the marine current turbines and their hubs and blades are totally submerged—in Appendix A.6, Appendix A.7 and Appendix A.8 the inertial characteristics of the bodies and their submerged volumes are shown.

$$V_{SUM}\left(t\right)={\displaystyle \sum }_{i=1}^n{V}_i\left(t\right)$$

(19)

$$m_{SUM}\left(t\right)=V_{SUM}\left(t\right)·{\rho }_{SEA WATER}$$

(20)

The calculation of Equation (19) is made by dividing the floating platform into three regions: a lower cylindrical zone with a diameter of 9.4 m and an upper one, also cylindrical with a diameter of 6.5 m, these zones are joined by a linearly tapered conical region—in [49], all these features are described at greater length. The calculation of the submerged volume of the floating platform is conducted in a simple way by obtaining the volume of the solids in these three regions. When the platform is tilted, the volume of the upper cylinder is obtained by considering the volume of a truncated cylinder.

The calculation of the volume Equation (19) is completed by adding the volume of the floating platform to the volume of the support for the marine current turbines, their hubs, and their blades.

Using Equation (20), the mass of the submerged volume is obtained by multiplying it by the density of seawater—1025 kg/m³.

To calculate the added mass, it is necessary to know the position of the center of buoyancy of the floating system Equation (21). This is the center of gravity [66] of the submerged volume explained in previous paragraphs found at each instant of time. The calculation is performed in a similar way to the calculation of the center of mass of the floating system Equation (7).

$$\left[\begin{array}{c}Co{B}_x\left(t\right)\\ Co{B}_y\left(t\right)\\ Co{B}_z\left(t\right)\end{array}\right]=\left[\begin{array}{c}{x}_{CoB}\\ y_{CoB}\\ z_{CoB}\end{array}\right]=\frac{1}{V_{SUM}\left(t\right)}·{\displaystyle \sum }_{i=1}^n\left[\begin{array}{c}{V}_i\left(t\right)·Co{B}_{xi}\\ V_i\left(t\right)·Co{B}_{yi}\\ V_i\left(t\right)·Co{B}_{zi}\end{array}\right]$$

(21)

The volumes (V_i (t)) that appear in Equation (21) correspond to the volume of the lower cylinder, the upper cylinder, and the conical region of the floating platform, as well as the volumes of the support of the marine current turbines, their hubs, and their blades. The centers of buoyancy (CoB_i) correspond to those of mentioned volumes—when the platform is tilted, the center of buoyancy of the upper cylinder of the floating platform is obtained by considering the center of buoyancy of a truncated cylinder.

In addition to the volume and center of buoyancy of the zones indicated in the previous paragraph, for the processing of the added mass, it is necessary to know the resulting inertia tensor of all these zones expressed at the origin of the mobile coordinate system.

The inertia tensor of the floating platform is obtained—in the same way that it was conducted for the calculation of the volume and the center of buoyancy—by dividing the floating platform into the same three regions: lower cylinder, upper cylinder, and conical region. The calculation of the moments and products of inertia of these zones are obtained in a simple way, since they are very common solids [54,55,56,57,58,59].

The inertia tensors of the submerged volumes of the support of the marine current turbines, their hubs and their blades are indicated in Appendix A, Appendix A.6, Appendix A.7 and Appendix A.8.

Once the inertia tensors of all the submerged volumes are known, the necessary rotation and translation operations are carried out—in a similar way to that explained in Section 2.3.2. Equations (8)–(17)—to finally obtain the inertia tensor expressed at the origin of the mobile coordinate system Equation (22).

$$M_{IT \left(SUM\right)\left(FS\right)}\left(t\right)=\left[\begin{array}{ccc}{I}_{xx \left(SUM\right)}\left(t\right)& I_{xy \left(SUM\right)}\left(t\right)& I_{zx \left(SUM\right)}\left(t\right)\\ I_{xy \left(SUM\right)}\left(t\right)& I_{yy \left(SUM\right)}\left(t\right)& I_{yz \left(SUM\right)}\left(t\right)\\ I_{zx \left(SUM\right)}\left(t\right)& I_{yz \left(SUM\right)}\left(t\right)& I_{zz \left(SUM\right)}\left(t\right)\end{array}\right]$$

(22)

Knowing that the added mass must be interpreted as a particular volume of fluid particles that are accelerated with the body [65], the movement of this volume of particles depends on the shape of the body [64]. For this reason, the values of the submerged mass Equation (20) and the inertia tensor Equation (22) cannot be used directly without considering the shape of the body.

This consideration is complicated to carry out for the real shape of the submerged volume of a body; in [61,64], the process for calculating the constants called added mass derivatives is exposed.

In this work, we chose to carry out a simplification that allows us to process the added mass with a negligible error and a low computational cost. This simplification was made by resembling the submerged volume of the floating system—to calculate the constants that describe the shape of the body—to that of an ellipsoid as shown in Figure 6.

The values of the semi-axes of the ellipsoid chosen are the following: a and b were assigned the value of the radius of the lower cylinder of the floating platform (4.7 m) and c was assigned half the length of the floating platform that is submerged at each instant of time (60 m in its initial position).

In [64], the equations to obtain the added mass derivatives of an ellipsoid and a prolate spheroid—when a = b and c > a—are exposed. In these equations appears the mass and moments of inertia of the ellipsoid multiplied by a factor that considers the shape of the ellipsoid. The calculation of these factors is made from three parameters, ${\alpha }_0$, ${\beta }_0$ and ${\gamma }_0$ Equations (24) and (25), which depend on the values of the semi-axes of the ellipsoid. The parameters ${\alpha }_0$, ${\beta }_0$ and ${\gamma }_0$ are obtained from the eccentricity of the meridian elliptical section Equation (23) [64,67].

$$e=\sqrt{1-\frac{a^2}{c^2}}$$

(23)

$${\alpha }_0={\beta }_0=\frac{1}{e^2}-\frac{1-e^2}{2·e^3}·\log \frac{1+e}{1-e}$$

(24)

$${\gamma }_0=\frac{2·\left(1-e^2\right)}{e^3}·\left(\frac{1}{2}·\log \frac{1+e}{1-e}-e\right)$$

(25)

In [67], the mathematical development of the ${\alpha }_0$, ${\beta }_0$ and ${\gamma }_0$ factors are described in greater detail until reaching the approaches expressed in [64].

Using Equations (26)–(28), the mass of the floating system corrected for the shape of the submerged volume is obtained. They are the same as indicated in [64], but substituting the mass of the ellipsoid for the mass of the submerged volume of the floating system.

$$X_{\dot{u}}=-\frac{{\alpha }_0}{2-{\alpha }_0}·m_{SUM}\left(t\right)$$

(26)

$$Y_{\dot{v}}=-\frac{{\beta }_0}{2-{\beta }_0}·m_{SUM}\left(t\right)$$

(27)

$$Z_{\dot{w}}=-\frac{{\gamma }_0}{2-{\gamma }_0}·m_{SUM}\left(t\right)$$

(28)

In the same way, with Equations (29)–(31) the moments of inertia of the submerged volume of the floating system corrected according to the shape of this volume are obtained. Equations (29)–(31) are based on those indicated in [64] for an ellipsoid but substituting the moments of inertia of the ellipsoid for those of the submerged volume of the floating system. For the corrections expressed in Equations (26)–(31), the submerged volume is considered to be an ellipsoid.

$$K_{\dot{p}}=-\frac{{\left(b^2-c^2\right)}^2·\left({\gamma }_0-{\beta }_0\right)}{2·\left(b^4-c^4\right)+{\left(b^2+c^2\right)}^2·\left({\beta }_0-{\gamma }_0\right)}·I_{xx \left(SUM\right)}\left(t\right)$$

(29)

$$M_{\dot{q}}=-\frac{{\left(c^2-a^2\right)}^2·\left({\alpha }_0-{\gamma }_0\right)}{2·\left(c^4-a^4\right)+{\left(c^2+a^2\right)}^2·\left({\gamma }_0-{\alpha }_0\right)}·I_{yy \left(SUM\right)}\left(t\right)$$

(30)

$$N_{\dot{r}}=-\frac{{\left(a^2-b^2\right)}^2·\left({\beta }_0-{\alpha }_0\right)}{2·\left(a^4-b^4\right)+{\left(a^2+b^2\right)}^2·\left({\alpha }_0-{\beta }_0\right)}·I_{zz \left(SUM\right)}\left(t\right)=0$$

(31)

Finally, all the calculation carried out in the added mass processing is collected into a matrix Equation (32)—added mass matrix—(M_AM) thus facilitating the use of these data in the calculation of dynamics. The added mass matrix was created in a similar way to the rigid body Matrix (18).

$$M_{AM}=-\left[\begin{array}{cc}\begin{array}{c}{X}_{\stackrel{·}{u}}\\ 0\\ 0\end{array}& \begin{array}{c}0\\ Y_{\stackrel{·}{v}}\\ 0\end{array}\\ \begin{array}{c}0\\ m_{SUM}·z_{CoB}\\ -m_{SUM}·y_{CoB}\end{array}& \begin{array}{c}-m_{SUM}·z_{CoB}\\ 0\\ m_{SUM}·x_{CoB}\end{array}\end{array} \begin{array}{cc}\begin{array}{c}0\\ 0\\ Z_{\stackrel{·}{w}}\end{array}& \begin{array}{c}0\\ -m_{SUM}·z_{CoB}\\ m_{SUM}·y_{CoB}\end{array}\\ \begin{array}{c}{m}_{SUM}·y_{CoB}\\ -m_{SUM}·x_{CoB}\\ 0\end{array}& \begin{array}{c}{K}_{\stackrel{·}{p}}\\ -I_{xy\left(SUM\right)}\\ -I_{zx\left(SUM\right)}\end{array}\end{array} \begin{array}{cc}\begin{array}{c}{m}_{SUM}·z_{CoB}\\ 0\\ -m_{SUM}·x_{CoB}\end{array}& \begin{array}{c}-m_{SUM}·y_{CoB}\\ m_{SUM}·x_{CoB}\\ 0\end{array}\\ \begin{array}{c}-I_{xy\left(SUM\right)}\\ M_{\stackrel{·}{q}}\\ -I_{yz\left(SUM\right)}\end{array}& \begin{array}{c}-I_{zx\left(SUM\right)}\\ -I_{yz\left(SUM\right)}\\ N_{\stackrel{·}{r}}\end{array}\end{array}\right]$$

(32)

The mass and moments and products of inertia appearing in Equations (22) and (32) are time dependent, the notation in the expressions has not been included for clarity.

In relation to the simplifications made by resembling the shape of the submerged volume to that of an ellipsoid, the behavior of the floating system was compared with numerous simulations performed under the same conditions using the FAST software, obtaining a high coincidence. The result of this validation is offered in Appendix E.

2.3.4. Dynamics

The calculation of the dynamics corresponds to the one explained in [61], which is based on the application of Newton’s second law to a floating system. First, the acceleration vector Equation (33) is calculated from the rigid body matrix, the added mass matrix and the resulting vector of forces acting on the floating system expressed in the mobile coordinate system [61]. The resultant vector of forces is a vector of six components, three forces, and three moments taken about the X, Y, and Z axes of the mobile coordinate system. So, the acceleration vector is another vector also with six components, three linear and three angular accelerations expressed in the mobile coordinate system.

$${\overrightarrow{av}}_{i \left(BODY\right)}={\left(M_{RB}+M_{AM}\right)}^{-1}·{\overrightarrow{F_i}}_{\left(BODY\right)}^{ TOTAL}={\left(M_{RB}+M_{AM}\right)}^{-1}·{\left[F_x F_y F_z M_x M_y M_z\right]}^T$$

(33)

Equation (33) is not considered correct expressed in this way, since Newton’s second law is only valid when applied to inertial systems [68]. To validate Equation (33), the fact that the resultant of the forces has been expressed in a non-inertial system, such as the mobile coordinate system, must be compensated for.

The compensation is performed by calculating the value of the acceleration vector expressed in the inertial system from the acceleration expressed in the mobile system Equation (34) [61].

$${\overrightarrow{av}}_{i \left(INERTIAL\right)}=\left[\begin{array}{cc}{M}_{BODY}^{INERTIAL}& 0_{3x3}\\ 0_{3x3}& M_{ANG}{}_{BODY}^{INERTIAL}\end{array}\right]·{\overrightarrow{av}}_{i \left(BODY\right)}+\left[\begin{array}{cc}\frac{d}{dt}{M}_{BODY}^{INERTIAL}& 0_{3x3}\\ 0_{3x3}& \frac{d}{dt} M_{ANG}{}_{BODY}^{INERTIAL}\end{array}\right]·{\overrightarrow{vv}}_{i \left(BODY\right)}$$

(34)

Through Equations (35) and (36) the linear velocity transformation matrix and the angular velocity transformation matrix are obtained, respectively [61]. These matrices consider the successive rotations of the linear and angular components of the acceleration vector between the mobile and inertial coordinate systems. The derivatives of these matrices also appear in Equation (34), these are obtained—element by element—applying the numerical differentiation method of Richardson’s extrapolation [69]. The equations used to apply this method are shown in Appendix C.1.

$$M_{BODY}^{INERTIAL}=M_{RZ{q}_6}·M_{RY{q}_5}·M_{RX{q}_4}=\left[\begin{array}{ccc}\cos q_6& -\sin q_6& 0\\ \sin q_6& \cos q_6& 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{ccc}\cos q_5& 0& \sin q_5\\ 0& 1& 0\\ -\sin q_5& 0& \cos q_5\end{array}\right]·\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos q_4& -\sin q_4\\ 0& \sin q_4& \cos q_4\end{array}\right]$$

(35)

$$M_{ANG}{}_{BODY}^{INERTIAL}=\left[\begin{array}{ccc}1& \sin q_4·\tan q_5& \cos q_4·\tan q_5\\ 0& \cos q_4& -\sin q_4\\ 0& \frac{\sin q_4}{\cos q_5}& \frac{\cos q_4}{\cos q_5}\end{array}\right]$$

(36)

Integrating the acceleration vector expressed in the inertial coordinate system, the velocity vector expressed in the same system is obtained Equation (37). Integrating again, the vector with the new position of the floating system is obtained Equation (38). The integrals are performed by applying the numerical integration method of the Romberg algorithm [69]. The equations used to apply this method are explained in Appendix C.2.

$${\overrightarrow{vv}}_{i \left(INERTIAL\right)}={\displaystyle \int }{\overrightarrow{av}}_{i \left(INERTIAL\right)}·dt$$

(37)

$${\overrightarrow{q}}_{i \left(INERTIAL\right)}={\displaystyle \int }{\overrightarrow{vv}}_{i \left(INERTIAL\right)}·dt$$

(38)

The velocity vector expressed in the mobile coordinate system is calculated by applying Equation (39) [61]; these data are needed for the calculation of Equation (34) and also allows us to obtain the kinetic energy that reaches the floating system Equation (40) and the added mass Equation (41) [61].

$${\overrightarrow{vv}}_{i \left(BODY\right)}={\left[\begin{array}{cc}{M}_{BODY}^{INERTIAL}& 0_{3x3}\\ 0_{3x3}& M_{ANG}{}_{BODY}^{INERTIAL}\end{array}\right]}^{-1}·{\overrightarrow{vv}}_{i \left(INERTIAL\right)}$$

(39)

$$T_{RB}=\frac{1}{2}·{\left({\overrightarrow{vv}}_{i \left(BODY\right)}\right)}^T·M_{RB}·{\overrightarrow{vv}}_{i \left(BODY\right)}$$

(40)

$$T_{AM}=\frac{1}{2}·{\left({\overrightarrow{vv}}_{i \left(BODY\right)}\right)}^T·M_{AM}·{\overrightarrow{vv}}_{i \left(BODY\right)}$$

(41)

The processing of the kinematics, the inertial data, the added mass, the dynamics, and the external forces is completed for each instant of time, being related as shown in Figure 5. The total time to be simulated is divided into time increments, performing the processing described in each of these time instants. This value must be small enough so that errors made both in the numerical differentiation and in the numerical integration can be considered negligible [69]. The chosen value with which acceptable results and a low computational cost are obtained is 0.1 s.

2.4. Influential Forces on the Floating Hybrid System

The simulation of the behavior of the floating system focuses on solving the problem of aero-hydro dynamics coupling [70]. Obtaining the resulting forces acting on the floating system is based on the fact that those forces do not change too much at each instant of time. This assumption implies that the accelerations obtained from these forces must not vary significantly in the same period. If this assumption is given, it will be possible to apply the superposition theorem [51], calculating each of the force vectors separately, adding them vectorially to obtain the resulting force vector.

The calculation of the different forces acting on the floating system is generally easier if it is conducted by expressing the force vector in the mobile coordinate system. Obtaining the resultant of forces in this coordinate system and processing the dynamics as indicated in the previous section, adequate results are reached [61]. In this work, this procedure was adopted for simplicity and coherence with the methodology explained in [61].

$${\overrightarrow{F_i}}_{\left(BODY\right)}^{ TOTAL}={\overrightarrow{F_i}}_{\left(BODY\right)}^{ WIND TURBINE}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ GRAVITY}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ HYDRODYNAMICS}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ MOORING SYSTEM}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ CORIOLIS}$$

(42)

The calculation of the forces acting on the floating system can be summarized in two actions: first, obtain the vector of forces and the point of application of the force and second, calculate the force and the equivalent moment at the origin of the mobile coordinate system. Once all the forces are expressed at this point, the resultant vector Equation (42) is calculated, which is used in Equation (33) to calculate the acceleration vector and with it complete the processing of dynamics. The detailed explanation of the calculation of each of these forces will be presented in the second part of this paper.

The resultant vector of forces Equation (42) is obtained from the forces caused by the wind turbine, gravity, hydrodynamics, the mooring system, and Coriolis acceleration.

The vector of gravitational forces is calculated by applying Newton’s law to the floating system, using the total mass of it (m_FS) calculated with Equation (6) and the gravity acceleration (

Table 2. Simulation environment constants.

Property	Value	Symbol
Gravity Acceleration 1	9.80665 m/s2	g
Density of Air 1	1.225 kg/m3	ρAIR
Density of Seawater 1	1025 kg/m3	ρSEA WATER

¹ Values used by FAST8.

). Taking the center of mass of the floating system (CoM) obtained with Equation (7) as the point of application of the force, the moment vector of the force [54,55,56] expressed at the origin of the mobile coordinate system is found.

The vector of forces of the mooring system is obtained by calculating the force exerted by each mooring line separately and then adding them: Equation (43). The properties of the mooring system appear in [49].

$${\overrightarrow{F_i}}_{\left(BODY\right)}^{ MOORING SYSTEM}={\overrightarrow{F_i}}_{\left(BODY\right)}^{ LINE 1}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ LINE 2}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ LINE 3}$$

(43)

The vertical and horizontal force exerted by each mooring line is calculated by applying the methodology described in [51,71], from which the resultant force of each line is obtained. Taking the fairlead of each line (p_Fi) as the point of application of the force (

Table 3. Mooring system significant values.

Element	Symbol	Position (m)
Fairlead Mooring Line 1	pF1	(5.2, 0, −70)
Anchor Mooring Line 1	pA1	(853.87, 0, −320)
Yaw angle with respect to inertial X axis	φLine1	0 deg
Fairlead Mooring Line 2	pF2	(–2.6, 4.503, −70)
Anchor Mooring Line 2	pA2	(–426.9, 739.5, −320)
Yaw angle with respect to inertial X axis	φLine2	120 deg
Fairlead Mooring Line 3	pF3	(–2.6, −4.503, −70)
Anchor Mooring Line 2	pA3	(–426.9, −739.5, −320)
Yaw angle with respect to inertial X axis	φLine3	240 deg

shows these values for the initial position of the floating system), the moment vector of the force expressed at the origin of the mobile coordinate system is found.

The vector of Coriolis forces is calculated from the velocity vector obtained with Equation (39), applying the criteria explained in [61].

The force vector of wind turbine is obtained from the thrust generated by the turbine; this is calculated by applying the one-dimensional theory [72,73].

Taking the center of mass of the turbine (CoM_WTurbine(t)) as the point of application of the force (this is obtained by applying Equation (7) from the centers of mass of the hub (

Table A7. Inertial properties of the wind hub.

Element	Wind Hub (Values Expressed in Local Coordinate System)	Wind Hub
Density (kg/m3)	6677	6677
Mass (kg)	56,781.9	56,781.9
Ixx (kg·m2)	186,851	188,149
Iyy (kg·m2)	357,732	357,732
Izz (kg·m2)	357,732	356,434
Ixy (kg·m2)	0	0
Iyz (kg·m2)	0	0
Izx (kg·m2)	0	−14,836.6
Center of Mass
Surge (m)	−2.28302	−5.00102
Sway (m)	0	0
Heave (m)	0	90.0001

) and of the three blades of the wind turbine (

Table A8. Inertial properties of the wind blade.

Element	Wind Blade (Values Expressed in Local Coordinate System)	Blade 1	Blade 2	Blade 3
Density (kg/m3)	599.55	599.55	599.55	599.55
Mass (kg)	17,740.2	17,740.2	17,740.2	17,740.2
Ixx (kg·m2)	3,118,006	3,108,994	3,100,328	3,101,268
Iyy (kg·m2)	3,106,302	3,106,302	800,161	803,206
Izz (kg·m2)	30,740.4	39,751.7	2,354,560	2,350,574
Ixy (kg·m2)	3278.72	3346.11	199,425	−202,771
Iyz (kg·m2)	1616.54	1471.99	1,317,325	−1,318,797
Izx (kg·m2)	35,997.0	170,397	118,667	113,465
Center of Mass
Surge (m)	0.0881023	−4.14516	−7.01647	−7.01551
Sway (m)	−0.00636023	−0.00636023	−19.0142	19.0205
Heave (m)	20.4750	111.968	79.1493	79.1603

)—both tables can be found in Appendix A), the moment vector of the force expressed at the origin of the mobile coordinate system is found. To calculate the thrust of the turbines, the data from

Table 4. Wind turbine properties.

Property	Value	Symbol
Hub Height	90 m	zREF
Shaft Tilt	5 deg	shT
Precone	2.5 deg	preC
Gearbox Ratio	97:1	gearR
Electrical Generator Efficiency	0.944	genE
Blade Length	61.5 m	Lblade
Hub Radius	1.5 m	rhub

,

Table 5. Rotor speed and pitch angle as a function of wind speed.

Wind Speed (VREF)	Rotor Speed (Ω)	Pitch Angle (φPITCH)
3 m/s	6.97 rpm	0 deg
4 m/s	7.18 rpm	0 deg
5 m/s	7.51 rpm	0 deg
6 m/s	7.94 rpm	0 deg
7 m/s	8.47 rpm	0 deg
8 m/s	9.16 rpm	0 deg
9 m/s	10.3 rpm	0 deg
10 m/s	11.43 rpm	0 deg
11 m/s	11.89 rpm	0 deg
11.4 m/s	12.1 rpm	0 deg
12 m/s	12.1 rpm	3.82 deg
25 m/s	12.1 rpm	23.47 deg

and

Table 6. Generator power as a function of wind speed.

Wind Speed (VREF)	Generator Power (PGEN)
3 m/s	50 kW
5 m/s	400 kW
6.5 m/s	1000 kW
8.2 m/s	2000 kW
9.3 m/s	3000 kW
11.4 m/s	5000 kW
25 m/s	5000 kW

are used.

The data in Table 4 were extracted from [50]; those in Table 5 were extracted from [51], while those in Table 6 were deduced from the graphs of steady-state responses as a function of wind speed from [50].

The modeling of the wind was carried out in the uniform wind shear type [72], using the corresponding equations indicated in [74].

The hydrodynamic force vector is obtained by calculating the hydrostatic, additional damping, added mass, viscous drag, marine current turbines, and wave forces (44).

$${\overrightarrow{F_i}}_{\left(BODY\right)}^{ HYDRODYNAMICS}={\overrightarrow{F_i}}_{\left(BODY\right)}^{ HYDROSTATICS}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ ADDITIONAL DAMPING}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ ADDED MASS}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ VISCOUS DRAG}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ CURRENT TURBINE}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ WAVES}$$

(44)

The vector of hydrostatic forces is calculated by applying the Archimedes’ principle to the submerged volume of the floating system, using its total mass (m_SUM) calculated with Equation (20) and the gravity acceleration (Table 2). Taking the center of buoyancy of the floating system (CoB) obtained with Equation (21) as the point of application of the force, the moment vector of the force expressed at the origin of the mobile coordinate system is found.

The vector of forces of the additional damping is calculated from the velocity vector obtained with Equation (37), applying the additional damping matrix described in [31].

The vector of forces of the hydrodynamic added mass is included in Equation (33), through the added mass matrix (M_AM). If it is desired to obtain the value of the vector of forces and moments to analyze their behavior, they can be calculated as indicated in [61,64].

All the simulations represented in this work do not include the effect of the waves, so the vector of forces of the waves will be explained in the second part of this paper.

The viscous drag force vector is calculated along the entire floating system, that is, along the tower and the floating platform. To do this, the influence of the wind on the tower and the near and sub-surface current on the floating platform are calculated Equation (45).

$${\overrightarrow{F_i}}_{\left(BODY\right)}^{ VISCOUS DRAG}={\overrightarrow{F_i}}_{\left(BODY\right)}^{ WIND TOWER}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ NEAR-SURFACE CURRENT}+{\overrightarrow{F_i}}_{\left(BODY\right)}^{ SUB-SURFACE CURRENT}$$

(45)

To obtain the vector of viscous drag forces, the wind speed along the tower and the near and sub-surface current speeds along the floating platform are obtained by applying the equations described in [74].

Obtaining the speed of the floating system at the same points in which the speeds indicated in the previous paragraph are calculated, and applying the equations described in [75], the vector of viscous drag forces is obtained. Taking each of these points as the points of application of the force, the moment vector of the force expressed at the origin of the mobile coordinate system is found.

The force vector of marine current turbines is obtained, as in the wind turbine, from the thrust generated by the turbines. The thrust of the marine current turbines is calculated from the same reasoning used for the calculation of the thrust of the wind turbine.

Taking, in this case, the center of mass of the turbine (CoM_MCTurbine(t)) as the point of application of the force(this is obtained by applying Equation (7) from the centers of mass of the hub (

Table A11. Inertial properties of the marine current hub.

Element	Marine Current Hub (Values Expressed in Local Coordinate System)	Clockwise Hub	Counterclockwise Hub
Density (kg/m3)	8500	8500	8500
Mass (kg)	10,906.1	10,906.1	10,906.1
Ixx (kg·m2)	9399.47	9399.47	9399.47
Iyy (kg·m2)	10,295.1	10,295.1	10,295.1
Izz (kg·m2)	10,536.5	10,536.5	10,536.5
Ixy (kg·m2)	0	0	0
Iyz (kg·m2)	0	0	0
Izx (kg·m2)	0	0	0
Center of Mass
Surge (m)	−1.26718	−1.06718	−1.06718
Sway (m)	0	17.1	−17.1
Heave (m)	0	−20	−20

) and of the two blades of the marine current turbine (

Table A13. Inertial properties of the clockwise marine current blade.

Element	Clockwise Blade (Values Expressed in Local Coordinate System)	Blade 1	Blade 2
Density (kg/m3)	1578	1578	1578
Mass (kg)	3213.43	3213.43	3213.43
Ixx (kg·m2)	15,720.1	15,720.1	15,720.1
Iyy (kg·m2)	15,493.9	15,493.9	15,493.9
Izz (kg·m2)	298.420	298.420	298.420
Ixy (kg·m2)	44.7426	44.7426	−44.7426
Iyz (kg·m2)	77.0487	77.0487	77.0487
Izx (kg·m2)	124.715	124.715	−124.715
Center of Mass
Surge (m)	−0.0731756	−0.973176	−0.973176
Sway (m)	−0.105801	16.9942	17.2058
Heave (m)	3.46741	−15.5326	−24.4674

and

Table A15. Inertial properties of the counterclockwise marine current blade.

Element	Counterclockwise Blade (Values Expressed in Local Coordinate System)	Blade 1	Blade 2
Density (kg/m3)	1578	1578	1578
Mass (kg)	3213.43	3213.43	3213.43
Ixx (kg·m2)	15,720.1	15,720.1	15,720.1
Iyy (kg·m2)	15,493.9	15,493.9	15,493.9
Izz (kg·m2)	298.420	298.420	298.420
Ixy (kg·m2)	−44.7426	−44.7426	44.7426
Iyz (kg·m2)	−77.0487	−77.0487	−77.0487
Izx (kg·m2)	124.715	124.715	−124.715
Center of Mass
Surge (m)	−0.0731756	−0.973176	−0.973176
Sway (m)	0.105801	−16.9942	−17.2058
Heave (m)	3.46741	−15.5326	−24.4674

)—both tables can be found in Appendix A), the moment vector of the force expressed at the origin of the mobile coordinate system is found. To calculate the thrust of the turbines, the data from

Table 7. Marine current turbine properties.

Property	Value	Symbol
Water Depth	320 m	dTOTAL
Shaft Tilt	0 deg	shT
Precone	0 deg	preC
Gearbox Ratio	97:1	gearR
Electrical Generator Efficiency	0.944	genE
Blade Length	9 m	Lblade
Hub Radius	1 m	rhub

,

Table 8. Rotor speed and pitch angle as a function of marine current speed.

Marine Current Speed (VSWL)	Rotor Speed (Ω)	Pitch Angle (φPITCH)
0.5 m/s	3.37 rpm	0 deg
0.6 m/s	4.032 rpm	0 deg
0.7 m/s	4.694 rpm	0 deg
0.8 m/s	5.356 rpm	0 deg
0.9 m/s	6.018 rpm	0 deg
1.0 m/s	6.68 rpm	0 deg
1.1 m/s	7.35 rpm	0 deg
1.2 m/s	8.02 rpm	0 deg
1.3 m/s	8.69 rpm	0 deg
1.4 m/s	9.36 rpm	0 deg
1.5 m/s	10.03 rpm	0 deg
1.6 m/s	10.398 rpm	0 deg
1.7 m/s	10.765 rpm	0 deg
1.8 m/s	11.133 rpm	0 deg
1.9 m/s	11.5 rpm	0 deg
3 m/s	11.5 rpm	14.5 deg

and

Table 9. Generator power as a function of marine current speed.

Marine Current Speed (VSWL)	Generator Power (PGEN)
0.5 m/s	5.5 kW
1 m/s	199.964 kW
1.5 m/s	394.429 kW
1.9 m/s	550 kW
3 m/s	550 kW

are used.

In Table 7, the water depth is taken from [31], the shaft tilt, precone, blade length, and hub radius are derived from [32]. The same values are chosen for the gearbox ratio and the electrical generator efficiency as for the wind turbine. The data in Table 8 and Table 9 were derived from [32,52].

The modeling of the marine current, as in the case of the wind, was carried out in the uniform wind shear type [72], using the corresponding equations indicated in [74].

In this work, the implemented modeling of wind turbines and marine currents is a simple modeling that is based on one-dimensional theory. This theory allows us to obtain a good approximation in order to know which behavior the steady state response system will have.

If it is desired to obtain a mathematical model of the turbine that offers the evolution of the floating system during its transient state, the blade element momentum theory [72,73] should be used. In addition, by applying the blade element momentum theory, results are obtained from the characteristics of the aerodynamic airfoils and the dimensions and inclination in each section of the blades, the characteristics of the electric generator, and the pitch angle of each blade, allowing the implementation of individual pitch regulators, etc. Therefore, the output magnitudes obtained with turbines modeled using the blade element momentum theory offer much richer information. This modeling will be explained in the second part of this publication.

However, in the modeling of marine current turbines, more factors can be considered to obtain more complete results, for example, applying unsteady hydrodynamics [76] to the modeling of marine turbines or predicting the occurrence of cavitation [77,78].

3. Results

The tests carried out represent the natural behavior of the floating system, without an integrated control system working. The wind was modeled as a vector that affects the hub of the wind turbine and its entire area using a uniform wind shear model [72]. The marine current was modeled in the same way, affecting the entire area of the two marine turbines. In addition, the action of the marine current along the entire floating platform was modeled using a uniform current shear model.

All simulations were carried out with an 11.4 m/s wind speed and a 1.9 m/s marine current speed (maximum thrust conditions for turbines) and calm sea (without wave thrust). This is a theoretical situation, but it could occur in specific locations (Table 1). These values of both the wind and marine current were chosen to limit the floating hybrid system and thus are able to be used to study its behavior in unfavorable situations. These values have also been used in the validation of the model in Appendix E.

In Figures 7, 9, and 11, the thrust vectors and their incidence in the floating system in each test are indicated, and in Figures 8, 10, and 12, the respective results of the tests are shown. In the Supplementary Materials section, you can find a video of the behavior of the floating system of each test.

3.1. Combination of Thrusts with a Phase Angle of 0 Degrees (Test 1)

In this test, the wind and the marine current have the same direction (Figure 7). This mainly causes a high displacement in surge (Figure 8a) and a heeling angle in pitch (Figure 8b).

However, the values reached for the rest of degrees of freedom are not relevant except for a certain displacement in heave (Figure 8a).

As the wind and marine current have the same orientation, all the thrust of the three turbines is concentrated in a single direction; this causes the floating system to experience a large amount of kinetic energy (Figure 8f), which becomes a considerable surge displacement (Figure 8a). This situation means that mooring lines 2 and 3 (Figure 2) must withstand a force greater than 3,000 kN in each of their fairleads (Figure 8g). The force that the lines of the mooring system would have to support would be even higher if both the wind and the marine current were oriented 180 degrees, since, in this situation, it would only be the mooring line 1 that would have to support all the force.

Although the radius of the marine turbines is 10 m [32] and that of the wind turbine is 63 m [50], the thrust generated by each type of turbine is similar due to the differences in density between air and water (Figure 8c). However, regarding torque, great differences are appreciated since the torque depends to a great extent on the radius of the turbine (Figure 8d). This fact explains the differences in the electrical power generated (Figure 8e), since torque and power are directly related.

3.2. Combination of Thrusts with a Phase Angle of 180 Degrees (Test 2)

In this case, the wind and the marine current are in opposite directions (Figure 9). This causes the effect of the thrusts to be counteracted, softening the behavior of the floating system.

As the sum of the thrusts of the two marine turbines is greater than the thrust of the wind turbine, the floating system also experiences surge displacement (Figure 10a).

By counteracting the thrusts between the wind turbine and the marine turbines, the floating system behaves as if only a minor thrust is produced in the 0 degrees direction. Therefore, the kinetic energy that the floating system reaches is less than what would be expected if this thrust compensation did not occur (Figure 10f). Mooring lines 2 and 3 must withstand a greater force than mooring line 1 on their fairleads, but less than would be expected without this compensation (Figure 10g).

A striking result is that, being the sum of the thrust of the marine turbines greater than thrust of the wind turbine, the angle of heel is produced with a negative pitch (Figure 10b); this is because the center of rotation of the floating system is close to its center of mass (Table A1 of Appendix A) and the distance of the wind turbine from the center of rotation is greater than the distance of the marine turbines from it, which causes a much greater angular momentum produced by the wind turbine than that produced by the marine turbines.

3.3. Case 3. Combination of Thrusts with a Phase Angle of 90 Degrees (Test 3)

In the latter case, the wind and the marine current are in perpendicular directions (Figure 11). This generates important increases in all the degrees of freedom of the floating system.

The perpendicular orientation of the wind and the marine current produces a predominant surge displacement (Figure 12a), it also means that the thrusts between the marine and wind turbines are not canceled, which leads rotations in two degrees of freedom, roll and pitch (Figure 12b); this causes an angle of heel resulting from the combination of these. Relevant unwanted oscillations are also produced that can have a negative effect on the different subsystems of the floating system (Figure 12). The combination of wind and marine current vectors cause, in the mooring system, an unequal distribution of the forces in the fairlead of the lines, with line 3 bearing a greater load (Figure 12g).

3.4. Interpretation of Results and Experimental Conclusions

With these tests, the behavior of the floating system is exposed in three extreme combinations of wind vectors and marine currents (phase angles of 0, 90 and 180 degrees) with values of wind and marine current that generate the maximum thrust in each turbine (11.4 m/s wind speed and a 1.9 m/s marine current speed).

To analyze the results obtained in detail, it is interesting to compare the simulations together. In the first place, the angle of heel (Figure 13) is analyzed as it is one of the most important magnitudes to know, since factors as important as the generation of energy and the structural fatigue of the floating system depend on it. Appendix D shows the calculation method used for the direction and angle of heel.

Comparing the heeling angles of the three tests (Figure 13a), a greater heeling angle is observed when the thrust of the wind turbine and that of the marine turbines are added (Test 1), and a lower heeling angle when the thrusts of the turbines counteract each other (Test 2). This allows us to define, for future work, what can be considered a favorable and unfavorable combination of operation of the floating system, to correct or take advantage of the situation through an automatic control strategy.

When the wind and marine current are perpendicular (Test 3), the floating system reaches an intermediate angle of heel. Although this situation seems more favorable than that of Test 1, the existing oscillations in the response can negatively affect the integrity of the floating system. These oscillations also appear in the responses of the rest of the magnitudes of Test 3; minimizing them is another issue to consider in the design of the control strategy.

Another interesting comparison is knowing the distance that the floating system reaches from its initial position (Figure 14a). This was calculated by obtaining the module of the position vector whose caomponents are the displacements in surge, sway, and heave. These displacements are related to the forces that the mooring system must withstand. If there is more separation from the initial position, more elongation of the mooring lines occurs.

Figure 14b shows the responses of line 3 for the three tests; this line presents the greatest effort in each test. Comparing the responses in Figure 14a,b, it is verified that a greater distance from the initial position produces a greater stress in the mooring lines.

Analyzing the responses carefully, a high stress is observed in the mooring line of Test 3. The explanation is as follows: In Tests 1 and 2, the wind and the marine current influence in such a way that, in the floating system, there is always 2 mooring lines sharing the efforts—mooring lines 2 and 3 (Figure 2, Figure 8g and Figure 10g)—if the responses of Tests 1 and 2 in Figure 14 are observed, it is verified that the differences—when the responses reach the final value—between the distances to the initial position and the forces of the mooring system can be considered quasi-proportional. On the other hand, the combination of the wind and the marine current in Test 3 means that the efforts are not distributed between two mooring lines, but line 3 supports a greater force than the other two (Figure 12g). This situation must also be considered as an unfavorable situation, when all the effort to retain the floating system falls on one of the mooring lines. This must be considered when designing the strength of the mooring lines and can also be considered when designing the automatic control strategy.

The most important comparison is that of electrical power production (Figure 15a) since this is the reason for the existence of the floating hybrid system. In Test 2, it is when more electrical power is generated, and in Test 1, it is when less electrical power is generated, Test 3 being an intermediate case. Figure 15b allows us to compare the kinetic energy developed by the floating system in each case; these data are related to the movement of the floating system until it has reached a final position, for a determined orientation of the wind and marine current.

From the results, it can be deduced that the ideal operating situation arises when the wind and the marine current have opposite directions and the orientation between the two fluids causes the stresses to be distributed on two of the mooring lines. On the contrary, the most unfavorable situation is when the wind and the marine current have the same direction and their orientation means that one of the three mooring lines must retain the floating system.

Intermediate situations between these two extremes produce moderate situations, but it is interesting to analyze them from the point of view of frequency since there were oscillations in the similar response analyzed (Test 3) that could cause structural problems. In these situations, the most common, the design of an adequate control strategy contributes to minimizing damage to the floating system, maximizing energy production.

4. Discussion

The main objective of the development of the mathematical model is to analyze the stability and energy production of the floating system under different conditions to develop automatic control strategies that improve the behavior of the system, from the point of view of structural fatigue and the maximum energy generation.

In this paper, a large part of the developed mathematical model is explained, leaving for the second part of this paper the detailed exposition of the modeling of the influencing forces on the floating hybrid system. The modeling techniques explained in this work allow us to obtain, in a simple way, the modeling of any mechanical system. This is of great help in control engineering when reliable modeling is needed without going too deep into the mechanical analysis.

The results presented in Section 3 were obtained by applying the one-dimensional theory to model the thrust of wind and marine current turbines. In the second part of this paper, the modeling of turbines using One-dimensional theory and Blade element momentum theory will also be explained. The Blade element momentum theory will permit the evaluation of the floating hybrid system applying different geometries, aerodynamic airfoils, including or eliminating certain elements, etc., also thanks to the availability of the inertial characteristics of the single bodies that make up the floating system (Appendix A) and a simple procedure (explained in Section 2.3.2) to process these characteristics and obtain the global inertia tensor of the system.

Therefore, the tests carried out in Section 3 represent the natural behavior of the floating system when it is subjected to the action of the wind and the marine current, without an integrated control system working to minimize structural stress while maximizing energy production. In future work, these control strategies and the results obtained will be exposed.

The programming of the FHYGSYS was conducted in Matlab^® as it is an environment oriented to matrix calculation, although it could have been conducted with any other programming language that would allow a similar ease of matrix calculation. The programming was organized applying the clean code techniques exposed in [79]. Without the application of these programming techniques, the realization of FHYGSYS would have been much more complicated.

From the analysis of the tests carried out in Section 3, it follows that the smaller the angle of heel, the greater the energy production. As this is the main objective of the floating hybrid system, the control system will have to use the angle of heel as the main variable that it will have to minimize. The orientation of the wind and the marine current should also be considered to sacrifice energy production, if necessary, when the hybrid system approaches the worst-case situation described in the previous section. In this sense, the set of turbines can work in active stabilization mode to help, as a cooperative control, for the stabilization of the system. Regarding the possibilities of actuation, the torque of the turbines must also be considered—which can be controlled for low fluid speeds—and the pitch angle that can be used in any range of fluid speeds.

As a complement to the control system, to help increase the stability of the floating hybrid system, the possibility of increasing the mooring lines can be considered.

As future work in a forthcoming publication, cooperative control techniques that make compatible the achievement of maximum structural stability and optimal performance of the generation of the floating system will be addressed.

While the proposed model uses a horizontal WT, using a vertical WT is a possibility that could be explored in the future. In terms of power generation, horizontal WT works better [80], and it is more efficient because there are few wind angle changes. This is the most probable scenario in the locations described in the introduction with high wind speed and high marine currents. On the other hand, vertical WT has half the weight for the same power, which is an advantage for stability, and it works better with the wind changing directions. Future works will explore the use of vertical WTs in order to test the benefits of this kind of turbine.