**1. Introduction**

The application of traditional energy resources (e.g., fossil fuels) causes serious environmental pollution. Renewable energy resources, such as wind, sunlight, ocean currents, waves, and tidal current, have been introduced as alternatives for achieving a low-carbon society. Ocean current energy is a potential power source that must be developed. Various forms of ocean energy are being investigated as potential sources of power generation [1–4]. For example, the Kuroshio Current—a strong current passing through the east of Taiwan—is expected to be an excellent energy resource. It has a mean velocity of 1.2–1.53 m/s near the surface and the potential electricity capacity of the Taiwan Current is approximately 4 GW [5]. This ocean energy source is stable and abundant, with the potential to be developed and utilized.

Tidal current—which can be extracted from the rise and fall of sea levels under the gravitational force exerted by the Moon and Sun as well as Earth's rotation—is one of the most valuable resources. Moreover, tidal current energy is more predictable than wind and wave energies [6]. Tidal current turbine (TCT) can be categorized into horizontal- and vertical-axis tidal turbines [7,8]. These inventions can harness the kinetic energy of tides and principally convert it to electricity. Horizontal-axis TCTs are the most common device, with their rotation axis parallel to the current stream direction [9,10]. By contrast, vertical-axis TCTs rotate about a vertical axis perpendicularly to the current stream [11]. Chen and Lam [6] reviewed the survivability of tidal power devices used to harness tidal power.

Zhou et al. [7] presented up-to-date information on large tidal turbine projects with a power exceeding 500 kW as well as their achievements and development histories. Most industrialized marine current turbine (MCT) devices are horizontal-axis turbines, with the rotation axis parallel to the current flow direction. The main disadvantages associated with vertical-axis turbines include their relatively low self-starting capability, high torque fluctuations, and generally lower e fficiency than horizontal-axis turbines. The power of industrialized MCT devices rigidly fixed at a seabed below 80 m depth, such as Atlantis AR1000 turbine, Voith Hydro turbine, and GE-Alstom tidal turbine, is over 1 MW at a current speed of approximately 2.4–4 m/s. These are called seabed-mounted turbines. Some devices, such as the Scotrenewables SR250 turbine and Ocean Renewable Power Company (ORPC) turbine, are flexibly moored at deep seabeds.

A flexible moored device is an important tool for deployment in deep water. The majority of moored tidal current turbine developers agree that by using a flexibly moored system, the device will be automatically self-aligned to the direction of current flow [7,12–15]. A traditional design uses the gravity foundations or piles in deep water, which are complex and expensive. The flexible mooring lines and anchors can be deployed in deep water, where the other designs may by impractical [16]. Thus, it is important to develop a mathematical model for a flexible mooring system for ocean current energy systems. Muliawan et al. [17] determined the extreme responses in the mooring lines of a two-body floating wave energy converter with four catenary cables. Angelelli et al. [18] investigated the behavior of a wave energy convertor mooring system with four spread cables by using Ansys AQWA software. Chen et al. [4] investigated the wave-induced motions of a floating wave energy converter (WEC) with mooring lines by using the smoothed particle hydrodynamics method. Davidson and Ringwood [19] reviewed the mathematical models for wave energy convertor mooring systems. The marine energy developer Minesto [20] developed a floating subsea kite with flexible mooring.

Cribbs [7] proposed a conceptual design for the flexible mooring of a current turbine fixed to a 300-m-deep seabed. The system included a mooring chain, a mooring line, a flounder plate, two lines for the turbine and the platform, a marine turbine, and a rotating turbine using blade-estimated airfoils. However, this system has not been practically applied thus far. Chen et al. [4] successfully moored a 50 kW ocean current turbine supplied by Wanchi Company to an 850-m-deep seabed at the o ffshore area of Pingtung County, Taiwan. At a current speed of 1 m/s, the output power of the system was 26 kW. IHI and NEDO [21] conducted a demonstration experiment of the ocean current turbine located o ff the coast of Kuchinoshima Island, Kagoshima Prefecture, and obtained data for commercialization; the demonstration experiment was conducted for seven days. The turbine comprised a combination of three cylindrical floats, called pods, having a total length of approximately 20 m, width of approximately 20 m, and turbine rotor diameter of approximately 11 m. The turbine system was moored from the anchor installed on a 100-m-deep seabed.

Besides the current rotating-type turbine, the Wanchi Company is also developing an ocean current turbine with a translational blade. To produce more power from the ocean, it is configured with several matrix-array turbines. The system comprises a turbine, a buoyance platform, a traction rope, and a mooring foundation, as shown in Figure 1. The ocean current turbine is tethered to an ~900-m-deep seabed and the ocean current flows perpendicular to the turbine plane. This system might be more unstable than the horizontal rotational turbine. The system is composed of a turbine, a buoyance platform, a traction rope, and a mooring foundation. The system stability is important for the practical operation. So far, only a few studies have investigated the stability of the ocean current turbine system. No literature is devoted to the mathematical model of the system about the coupled heaven, surge and pitch motions. No analytical solution of the system is also presented. Because the system must be su fficiently stable under the e ffect of wave, the mathematical model is developed and the analytical solutions are presented in this study. Moreover, the e ffects of several parameters on the system stability are investigated.

In Section 2, a two-dimensional model for the motions of the ocean current turbine system and floater is developed. The turbine system and floater are treated as rigid bodies and the cables are divided into two sections. Theoretical solutions of the motions are presented in Section 3. As detailed in Section 4, a series of simulations were conducted for evaluating the dynamic stability of the system under various wave conditions. Section 5 summarizes the present study and provides some concluding remarks.

**Figure 1.** Current turbine system composed of turbine (**A**) (a: Generator; b: Float; c: Structure; d: Translational blade; e: Balanced weight; f: Cable; the current energy transferred by the blade to generator through the high-pressure oil); buoyance platform (**B**); traction rope (UHMWPE) (**C**); and anchor (**D**).

## **2. Governing Equations**

According to Figure 1, the analysis of the motion and dynamic stability of the ocean current turbine system with a flexible mooring cable is complex. Assume that the structures of turbine and carrier subjected the force due to wave and current are kept during the motion. Because the displacement and pitching motion of the system are significantly concerned, these two components are considered rigid bodies. Moreover, the inertia effect of the elastic cable is neglected. Therefore, the system is simulated as a discrete one. The overall system comprises an anchored mooring, a carrier, and a vertical ocean turbine system driven by a translational blade with a gravity anchor. A mooring system that comprises a long main cable with a sub-cable connecting the carrier and ocean turbine is deployed. Due to the complexity of system stability, the flow field is considered to be in steady state and the ocean current velocity is assumed to be constant and uniform. The horizontal force applied to the turbine and its structure during electricity generation is also assumed to be constant. The coupled motion of the system includes the horizontal, vertical, and pitching oscillations. As shown in Figures 2 and 3, owing to the wave fluctuation, the buoyance forces applied on the floating platform and turbine can be expressed as follows:

$$F\_{B1} = (H\_1 - X\_1)A\_1 \rho g \tag{1}$$

$$F\_{B2} = (H\_2 - X\_2)A\_2 \rho g \tag{2}$$

where the subscripts "1" and "2" denote the carrier and turbine, respectively, and *A*1 and *A*2 denote the corresponding hydrodynamic areas. The wave heights at the carrier and turbine are *H*1 = *H*0 sin Ω*t* and *H*2 = *H*0 sin(Ω*t* + <sup>Φ</sup>), where *H*0 is the wave amplitude. (*Hi* − *Xi*, I = 1~2) indicates the vertical displacement, Ω is the wave frequency, and Φ is the phase lag due to the propagation delay from the carrier to the turbine. On the basis of linear wave theory, the phase is expressed as Φ= −*kL*2, where *k* indicates the wave number, *L*2 is the horizontal distance between the carrier and turbine, ρ(=1025 kg · m<sup>−</sup>3) is the density of water, and g (=9.81 m · s<sup>−</sup>2) is the gravitational acceleration.

**Figure 2.** Concentrated mass model of the ocean current turbine system subjected to wave currents.

**Figure 3.** Coordinate and force distribution of the system.

*J. Mar. Sci. Eng.* **2020**, *8*, 687

According to the static equilibrium, the force in the horizontal direction at the joint o is expressed as follows:

$$T\_{wire} \cos \phi = T\_{turbine} \cos \phi \tag{3}$$

According to the dynamic equilibrium, the equations of motion of the floating carrier and turbine platform in the vertical direction are expressed as

..

..

$$M\_1 X\_1 - F\_{B1} + T\_{tur} \sin \phi - T\_{uire} \sin \theta = 0\tag{4}$$

and

..

..

$$M\_2 X\_2 - F\_{B2} - T\_{trr} \sin \phi = 0\tag{5}$$

where *X*1 and *X*2 are the vertical accelerations of the floater and turbine, respectively; *Ttur* is the tension of the wire between the turbine and floater; and *Twire* is that between the cable and anchor. If the vertical displacement is substantially smaller than the connecting cable between the floating and turbine platforms, the angle φ is very small, as shown in Figure 2, and can be approximated as follows:

$$
\sin \phi \approx \frac{X\_1 - X\_2}{L\_2} \tag{6}
$$

The cable made by polyethylene dyneema connecting the floating platform and mooring foundation is considered. The material is the ultra-high molecular weight polyethylene (UHMWPE). Because the material properties are grea<sup>t</sup> strength, light weight and flexible, the mooring cable is likely straight during the turbine subjected to the ocean current force. Moreover, because the vertical displacement is significantly smaller than the water depth, the following approximate relation can be obtained:

$$
\sin \theta \approx \frac{X\_0}{L\_1} \tag{7}
$$

Through substituting Equations (3), (6) and (7) into Equations (4) and (5), we obtain the coupled equations of motion in terms of vertical displacements {*X*1, *X*2} for the carrier and turbine:

$$\begin{aligned} &=T\_{\text{tur}}\frac{M\_1\ddot{X}\_1 + \left[A\_1\rho g + T\_{\text{tur}}\frac{1}{L\_2}\right]X\_1 - T\_{\text{tur}}\frac{1}{L\_2}X\_2}{\sqrt{L\_1^2 - X\_0^2}}\\ &=T\_{\text{tur}}\frac{X\_0}{\sqrt{L\_1^2 - X\_0^2}} + H\_1A\_1\rho g = T\_{\text{tur}}\frac{X\_0}{\sqrt{L\_1^2 - X\_0^2}} + A\_1\rho g H\_0 \sin\Omega t \end{aligned} \tag{8}$$

and

$$M\_2\ddot{X}\_2 - \frac{T\_{\text{tur}}}{L\_2}\mathbf{X}\_1 + \left(A\_2\rho g + \frac{T\_{\text{tur}}}{L\_2}\right)\mathbf{X}\_2 = A\_2\rho gH = A\_2\rho gH\_0\sin(\Omega t + \Phi) \tag{9}$$

Equations (8) and (9) can be expressed as the equation of vertical motion in the matrix form

$$
\begin{bmatrix} M\_1 & 0 \\ 0 & M\_2 \end{bmatrix} \begin{bmatrix} \ddot{X}\_1 \\ \ddot{X}\_2 \end{bmatrix} + \begin{bmatrix} K\_{11} & K\_{12} \\ K\_{21} & K\_{22} \end{bmatrix} \begin{bmatrix} X\_1 \\ X\_2 \end{bmatrix} = \begin{bmatrix} T\_{\text{Tur}} \frac{X\_0}{\sqrt{L\_1^2 - X\_0^2}} + A\_1 \rho g H\_0 \sin \Omega t \\\ A\_2 \rho g H\_0 \sin(\Omega t + \Phi) \end{bmatrix} \tag{10}
$$

where the first 2 × 2 matrix is the mass one, the second 2 × 2 matrix is the stiffness one, the last term is the forcing one due to the wave and the force of the turbine and

$$K\_{11} = A\_1 \rho g + T\_{\rm tur} \frac{1}{L\_2}, \; K\_{12} = K\_{21} = \frac{-T\_{\rm Tur}}{L\_2}, \; K\_{22} = A\_2 \rho g - \frac{T\_{\rm Tyr}}{L\_2},$$

where X0 is the depth of the seabed, {X1, X2} indicate the vertical displacements, Φ is the phase lag angle, and {M1, M2} are the masses.

Because the distance *eGB* between the centers of gravity and buoyance of the ocean turbine set is almost constant during the pitching motion, the righting moment is '*WeGB* sinψ' and the heeling angle curve can be expressed as '*eGB* sin ψ'. Based on the principle of dynamic equilibrium, the equation of horizontal translational motion for the turbine can be derived

..

$$M\_2 \varchi\_2 - F\_D + T\_{turbinu} \cos \phi = 0,\tag{11a}$$

where *FD* = *CD* 1 2ρ*<sup>A</sup> V* − . *Y*2 − . ψ*l* 2. Considering the velocity resulting from the oscillation . *Y*2 + . ψ*l* << *V*, the drag force becomes *FD* ≈ *CD* 1 2ρ*<sup>A</sup> V*<sup>2</sup> − 2*V* . *Y*2 + . ψ*l* .

Because the angle φ approaches zero and *CD* 1 2ρ*AV*<sup>2</sup> ≈ *Tturbine*, Equation (11a) becomes equivalent to the equation of horizontal motion:

$$M\_2\ddot{Y}\_2 + C\_D\rho A V \dot{Y}\_2 + C\_D\rho A V l \dot{\psi} = 0\tag{11b}$$

In the dynamic equilibrium equation, the pitching motion of the turbine is expressed as

$$I\ddot{\boldsymbol{\psi}} + \mathsf{W}\boldsymbol{e}\_{\mathsf{GB}}\sin\boldsymbol{\psi} + \mathsf{F}\mathsf{Det}\,\boldsymbol{\phi}\cos\boldsymbol{\psi} = (T\_{\text{Turbine}}\cos\phi)(\boldsymbol{\varepsilon}\mathbf{p}\cos\boldsymbol{\psi}) - (T\_{\text{Turbine}}\sin\phi)(\boldsymbol{\varepsilon}\mathbf{p}\sin\boldsymbol{\psi})\tag{12a}$$

Considering the pitching angle ψ to be small and based on Equation (6), the equation of the pitching motion becomes

$$
\dot{\Psi}\dot{\psi} - \mathbb{C}\_D \rho A V l e\_D \dot{\psi} + \left[ \mathcal{W} e\_{\rm GB} + T\_{\rm Turbine} e\_{\rm D} \left( \frac{X\_1 - X\_2}{L\_2} \right) \right] \dot{\psi} - \mathbb{C}\_D \rho A V e\_D \dot{\psi} = 0 \tag{12b}
$$

where *y* indicates the horizontal displacement, ψ indicates the pitch angle of the turbine, *A* is the area of drag, *CD* is the drag coe fficient, *eGB* is the distance between the centers of gravity and buoyance, *V* is the current velocity, and *l* is the radius of rotation.

#### **3. Theoretical Solutions**

#### *3.1. Solution of the Vertical Motion*

The solutions of vertical displacement of the turbine and carrier comprise static and dynamic components and can be expressed as

$$X\_1(t) = X\_{10} + X\_{11}(t), \\ X\_2(t) = X\_{20} + X\_{21}(t) \tag{13}$$

where {*X*10, *X*20} are the static displacements and{*X*11, *X*21} are the dynamic displacements. By substituting Equation (13) into Equation (10) and dividing it into the static and dynamic subsystems, we obtain the following:

Static subsystem:

$$
\begin{bmatrix} K\_{11} & K\_{12} \\ K\_{21} & K\_{22} \end{bmatrix} \begin{bmatrix} X\_{10} \\ X\_{20} \end{bmatrix} = \begin{bmatrix} T\_{Tur} \frac{X\_0}{\sqrt{L\_1^2 - X\_0^2}} \\ 0 \end{bmatrix} \tag{14}
$$

It should be noted that the Equation (14) demonstrates the relation between the vertical displacements of the turbine and carrier and the current velocity without the wave e ffect.

Dynamic subsystem:

$$
\begin{bmatrix} M\_1 & 0 \\ 0 & M\_2 \end{bmatrix} \begin{bmatrix} \bar{X}\_{11} \\ \bar{X}\_{21} \end{bmatrix} + \begin{bmatrix} K\_{11} & K\_{12} \\ K\_{21} & K\_{22} \end{bmatrix} \begin{bmatrix} X\_{11} \\ X\_{21} \end{bmatrix} = \begin{bmatrix} A\_1 \rho g H\_0 \sin \Omega t \\ A\_2 \rho g H\_0 \sin (\Omega t + \Phi) \end{bmatrix} \tag{15}
$$

*J. Mar. Sci. Eng.* **2020**, *8*, 687

It should be noted that the Equation (15) demonstrates the dynamic behavior subjected to the wave effect under the static equilibrium. From Equation (14), the theoretical solutions of the static equilibrium can be easily derived:

$$X\_{10} = \frac{T\_{Tur} \frac{X\_0}{\sqrt{L\_1^2 - X\_0^2}}}{\left(K\_{11} - K\_{12} \frac{K\_{21}}{K\_{22}}\right)}, \; X\_{20} = -\frac{K\_{21}}{K\_{22}} X\_{10} \tag{16}$$

Further, the harmonic response of the system under the wave effect will be obtained. The solutions of the dynamic state of Equation (15) are assumed as

$$
\begin{bmatrix} X\_{11} \\ X\_{21} \end{bmatrix} = \begin{bmatrix} X\_{11c} \\ X\_{21c} \end{bmatrix} \cos \Omega t + \begin{bmatrix} X\_{11s} \\ X\_{21s} \end{bmatrix} \sin \Omega t \tag{17}
$$

Substituting Equation (17) into Equation (15), and after some manipulations, the theoretical solutions of the dynamic state can be derived as follows:

$$
\begin{bmatrix} X\_{11c} \\ X\_{21c} \end{bmatrix} = \frac{\beta\_c}{|A|} \begin{bmatrix} -K\_{21} \\ K\_{11} - \Omega^2 M\_1 \end{bmatrix} \tag{18}
$$

and

$$
\begin{bmatrix} X\_{11s} \\ X\_{21s} \end{bmatrix} = \frac{1}{|A|} \begin{bmatrix} a\{K\_{22} - \Omega^2 M\_2\} - \beta\_5 K\_{21} \\ -aK\_{12} + \beta\_5 \{K\_{11} - \Omega^2 M\_1\} \end{bmatrix} \tag{19}
$$

where

$$|A| = \left| \begin{array}{c} K\_{11} - \Omega^2 M\_1 \\ K\_{21} \end{array} \begin{array}{c} K\_{12} \\ K\_{22} - \Omega^2 M\_2 \end{array} \right| , \alpha = A\_1 \rho g H\_0 \beta\_\delta = A\_2 \rho g H\_0 \sin \Phi\_\gamma \beta\_\delta = A\_2 \rho g H\_0 \cos \Phi. \tag{20}$$

As the determinant function of the frequency is equal to zero, |*A*| = 0, the resonance condition applies. The two resonant frequencies can be derived as

$$
\Omega\_{1,2}^2 = \frac{1}{2} \left[ \left( \frac{K\_{22}}{M\_2} + \frac{K\_{11}}{M\_1} \right) \pm \sqrt{\left( \frac{K\_{22}}{M\_2} + \frac{K\_{11}}{M\_1} \right)^2 - 4 \left( \frac{K\_{11}}{M\_1} \frac{K\_{22}}{M\_2} - \frac{K\_{12} K\_{21}}{M\_1 M\_2} \right)} \right] \tag{21}
$$

It is well known that if the wave frequency approaches to the natural frequency of the system, the resonance or instability will occur. However, based on Formula (21) one can tune the natural frequencies of system away the wave frequency to avoiding the instability of vertical vibration.

#### *3.2. Solutions of the Horizontal and Pitching Motions*

Consider the horizontal and pitching motions of the static and dynamic systems. The solutions of the horizontal displacement and pitching angle are expressed as

$$\begin{aligned} y &= y\_0 + y\_c \cos \Omega t + y\_s \sin \Omega t \\ \psi &= \psi\_0 + \psi\_c \cos \Omega t + \psi\_s \sin \Omega t \end{aligned} \tag{22}$$

where *y*0,ψ0! are the static horizontal translational and pitching displacements without the wave effect *yc*, *ys*,ψ*<sup>c</sup>*,ψ*s*! are the harmonic solutions. By substituting Equation (22) into Equation (11), we can obtain the following two equations:

$$
\rho\_1 y\_c + \rho\_2 y\_s + \rho\_3 \psi\_s = 0 \tag{23}
$$

$$p\_1 y\_5 + p\_2 y\_6 + p\_3 \psi\_6 = 0 \tag{24}$$

where

$$\begin{aligned} \rho\_1 &= -M\_2 \Omega^2, \rho\_2 = \Omega \mathbb{C}\_D \rho A V, \rho\_3 = \Omega \mathbb{C}\_D \rho A V l. \\\ p\_1 &= M\_2 \Omega^2, p\_2 = \Omega \mathbb{C}\_D \rho A V, p\_3 = \Omega \mathbb{C}\_D \rho A V l. \end{aligned} \tag{25}$$

By Substituting Equations (22) and (13) into Equation (12), after some manipulation, we obtain

$$\begin{aligned} & \mathfrak{W}\_{\rm GB}\psi\_{0} \\ &+ \begin{pmatrix} -\Omega^{2}I\psi\_{c} - \Omega\Gamma\_{D}\rho A V l e\_{D} \psi\_{s} + T\_{Tur}e\_{D} \Big( \frac{X\_{11} - X\_{21c}}{L\_{2}} \Big) \psi\_{0} \\ + \left( \mathbb{W}e\_{GB} + T\_{Tur}e\_{D} \frac{X\_{01} - X\_{02}}{L\_{2}} \right) \psi\_{c} - \Omega\Gamma\_{D}\rho A V e\_{D} \psi\_{s} \\ + \left( -\Omega^{2}I\psi\_{s} + \Omega\Gamma\_{D}\rho A V l e\_{D} \psi\_{c} + T\_{Tur}e\_{D} \Big( \frac{X\_{11c} - X\_{21c}}{L\_{2}} \Big) \psi\_{0} \\ + \left( \mathbb{W}e\_{GB} + T\_{Tur}e\_{D} \frac{X\_{01} - X\_{21}}{L\_{2}} \right) \psi\_{s} + \Omega\Gamma\_{D}\rho A V e\_{D} \psi\_{c} \\ + T\_{Tur}e\_{D} \Big( \frac{X\_{11} - X\_{21c}}{L\_{2}} \Big) \psi\_{c} \frac{1 + \cos 2\Omega t}{2} \\ + T\_{Tur}e\_{D} \Big( \frac{X\_{11} - X\_{21c}}{L\_{2}} \Big) \psi\_{s} \frac{1 - \cos 2\Omega t}{2} \\ + \left[ T\_{Tur}e\_{D} \Big( \frac{X\_{11} - X\_{21c}}{L\_{2}} \Big) \psi\_{c} + T\_{Tur}e\_{D} \Big( \frac{X\_{11c} - X\_{21c}}{L\_{2}} \Big) \psi\_{s} \right] \sin 2\Omega t = 0 \end{aligned} (26)$$

Equation (26) indicates that the horizontal and pitching displacements depend on the resonant frequencies Ω and vertical displacements {*X*1, *X*2}, which are time-dependent functions. Therefore, the characteristic equations can be derived using the orthogonality relation of {sin *<sup>n</sup>*Ω*t*, cos *m*Ω*t*} as follows:

First, integrating Equation (26) from 0 to period T (=2π/Ω), Equation (26) can be expressed as

$$q\_0\psi\_0 + q\_1\psi\_c + q\_2\psi\_s = 0\tag{27a}$$

where

$$q\_0 = We\_{\rm GB}, \\ q\_1 = \frac{1}{2} T\_{\rm Turn} \varepsilon\_D \left( \frac{X\_{11c} - X\_{21c}}{L\_2} \right), \\ q\_2 = \frac{1}{2} T\_{\rm Turn} \varepsilon\_D \left( \frac{X\_{11s} - X\_{21s}}{L\_2} \right). \tag{27b}$$

Second, multiplying Equation (26) with cos Ω*t* and integrating it from 0 to T, Equation (26) can be expressed as

$$r\_0\psi\_0 + r\_1\psi\_c + r\_2\psi\_s + r\_3\psi\_s = 0\tag{28a}$$

where

$$\begin{split} r\_0 &= -T\_{\rm tur} \varepsilon\_D (\frac{X\_{1lc} - X\_{2lc}}{L\_2})\_\prime \, r\_1 = -\left(-\Omega^2 I + W \varepsilon\_{\rm GB} + T\_{\rm tur} \varepsilon\_D \frac{X\_{1l} - X\_{2l}}{L\_2}\right) \\ r\_2 &= \Omega C\_D \rho A V l \varepsilon\_{\rm D} \, \ r\_3 = \Omega C\_D \rho A V \varepsilon\_{\rm D} .\end{split} \tag{28b}$$

Finally, multiplying Equation (26) with sin Ω*t* and integrating it from 0 to T, Equation (26) can be expressed as

$$s\_0\psi\_0 + s\_1\psi\_c + s\_2\psi\_s + s\_3\psi\_c = 0\tag{29a}$$

where

$$\begin{array}{c} s\_0 = T\_{\text{Tur}} e\_D (\frac{X\_{11i} - X\_{21}}{L\_2}) & s\_1 = \Omega C\_D \rho A V l \varepsilon\_{D\prime} \\ s\_2 = \left(-\Omega^2 I + W e\_{\text{GB}} + T\_{\text{Tur}} e\_D \frac{X\_{10} - X\_{20}}{L\_2}\right) & s\_3 = \Omega C\_D \rho A V \varepsilon\_D \end{array} \tag{29b}$$

Equations (23), (24) and (27)–(29) can be further expressed as

$$
\begin{bmatrix}
0 & 0 & o\_3 & o\_4 & o\_5 \\
0 & 0 & p\_3 & p\_4 & p\_5 \\
q\_1 & q\_2 & q\_3 & 0 & 0 \\
r\_1 & r\_2 & r\_3 & 0 & r\_5 \\
s\_1 & s\_2 & s\_3 & s\_4 & 0 \\
\end{bmatrix}
\begin{bmatrix}
\psi\_0 \\
\psi\_c \\
\psi\_s \\
y\_c \\
y\_s \\
\end{bmatrix} = 0 \tag{30a}
$$

The eigenvalues of Equation (30a) provide a measure of system stability. These values are characterized using a matrix, which is the Jacobian of the state of the system. The eigenvalues of the matrix are the characteristic roots of the state equation and can be determined using the roots of characteristic equations of the system.

$$
\begin{vmatrix}
0 & 0 & q\_3 & q\_4 & q\_5 \\
0 & 0 & p\_3 & p\_4 & p\_5 \\
q\_1 & q\_2 & q\_3 & 0 & 0 \\
r\_1 & r\_2 & r\_3 & 0 & r\_5 \\
s\_1 & s\_2 & s\_3 & s\_4 & 0
\end{vmatrix} = 0 \tag{30b}
$$

or

$$
\rho\_3(q\_1 r\_2 s\_4 p\_5) - \rho\_4(q\_1 r\_2 s\_3 p\_5 - s\_1 r\_2 q\_3 p\_5) + \rho\_5(q\_1 r\_2 s\_3 p\_4 - s\_1 r\_2 q\_3 p\_4 - r\_1 q\_2 p\_3 s\_4) = 0 \tag{30c}
$$

From the characteristic Equation (30c), we can determine the system's resonance under the effects of various environmental parameters.

#### **4. Numerical Results**

It is important to design the current turbine system to be stable in order to obtain good performance and control requirement. The theoretical solutions, including static and dynamic stability systems, are introduced in Section 3. Here, the behaviors of the turbines and floaters under the effects of wave and ocean currents are numerically analyzed. Figures 4–7 show the effects of wave frequency f, drag force *TTur*, and turbine area *A*2 on the amplitudes {*X*11, *X*21} of vertical vibration.

These amplitudes under the condition of "*TTur* = 40 tons, *A*1 = 29.9 m2, *A*2 = 80 m<sup>2</sup>#at various wave frequencies are shown in Figure 4, where the dynamic motions of the turbine and floater change with the wave frequency. A lower wave frequency causes a lager vertical displacement of the turbine and floater, and the excitation of the floater is larger than that of the turbine. The vertical vibrations have two peaks at the natural frequencies {0.751 Hz, 1.22 Hz}, which are almost the same for the turbine and floater. This is attributed to the resonance vertical motion. Figure 5 shows the effect of wave frequency on {*X*11, *X*21} under the condition of "*TTur* = 100 tons, *A*1 = 29.9 m2, *A*2 = 80 <sup>m</sup><sup>2</sup>#. The vertical vibrations have two peaks at the natural frequencies {0.761 Hz, 1.21 Hz}. The results in Figures 4 and 5 show that the pretension *TTur* has only a slight influence on the resonance condition.

**Figure 4.** Effect of the wave frequency f on the amplitudes {*X*11, *X*21} of vertical vibration. (*A*1 = 29.9 m2, *A*2 = 80 m2, *M*1 = 13.26 tons, *M*2 = 100 tons, *L*1 = 2780 m, *L*2 = 50 m, *X*0 = 850 m, *H*0 = 10 m, *TTur* = 40 tons).

**Figure 5.** Effect of wave frequency f on {*X*11, *X*21}. (*TTur* = 100 tons; other parameters are the same as those in Figure 4).

The effects of the hydrodynamic area of the turbine and floater on vertical excitation are shown in Figures 6 and 7, respectively. Figure 6 demonstrates the effect of wave frequency on {*X*11, *X*21} at {*TTur* = 40 tons, *A*1 = 29.9 m2, *A*2 = 120 m2}. The two natural frequencies at the resonance are {0.751 Hz, 1.49 Hz}. Compared to Figure 4, this figure shows that the turbine area has a substantial effect on the second natural frequency. However, the effect on the first natural frequency is negligible. The effect of the wave frequency on {*X*11, *X*21} at {*TTur* = 40 tons, *A*1 = 40 m2, *A*2 = 80 m2} is shown in Figure 7. The two natural frequencies are {0.871 Hz, 1.49 Hz}, which differ from the results shown in Figure 4. The resonance effect occurs at higher wave frequencies, which shows that the floater hydrodynamic area has a more substantial influence on the two natural frequencies. We can thus conclude that the larger the floater and turbine areas are, the higher the natural frequencies are. Thus, the effect of the floater area is more significantly greater than that of the turbine area.

It is well known that when dynamic stability is being considered, a larger- *eGB* between the gravity and buoyance centers is preferred. The critical distance *eGB* about the instability of pitching motion is investigated here. Figure 8 shows the variation in the critical distance *eGB* with the wave frequencies under three drag forces *TTur* at the moment of inertial *I* = 8.33 × 10<sup>8</sup>*kg* − *<sup>m</sup>*<sup>2</sup>. For the wave frequency over 0.03 Hz, the larger the drag force *TTur*, the longer the critical distance *eGB*,*critical*. In other word, for the larger the drag force *TTur* the longer distance is required for the dynamic stability. Moreover, for a wave frequency over 0.03 Hz, the critical distance, *eGB*, is less than one, except at a frequency of 1.9 Hz. Thus, if *eGB* is larger than 1 m, the instability will not occur at a wave frequency above 0.03 Hz. However, for a wave frequency under 0.03 Hz, *eGB* is larger than 2 m, and thus the ocean current system will be unstable. According to Figure 9, at *I* = 1.67 × 10<sup>6</sup>*kg* − *<sup>m</sup>*<sup>2</sup>, for the wave frequency under 1 Hz, the critical distance *eGB*,*critical* is larger than 2 m and the entire system cannot maintain dynamic stability. It is concluded from Figures 8 and 9 that the larger the moment of inertial *I*, the better the stability of system.

**Figure 6.** Effect of wave frequency f on {*X*11, *X*21}. (*A*2 = 120 m2; other parameters are the same as those in Figure 4).

**Figure 7.** Effect of the wave frequency f on {*X*11, *X*21}. (*A*1 = 40 m2; other parameters are the same as those in Figure 4).

**Figure 8.** Effect of the wave frequency f and drag force *TTur* on the critical distance between the centers of gravity and buoyance for resonance. (*CD* = 1.0, *V* = 1 m/s, *eD* = 1 m, = 1 m, *I* = 8.33 × 10<sup>8</sup> (kg − m<sup>2</sup>); other parameters are the same as those in Figure 4).

**Figure 9.** Effect of wave frequency f and drag force *TTur* on the critical distance between the centers of gravity and buoyance for resonance. (*I* = 1.67 × 10<sup>6</sup>(kg − m<sup>2</sup>); other parameters are the same as those in Figure 8).

## **5. Conclusions**

In this study, a mathematical model for a system comprising a turbine, buoyance platform, traction rope, and mooring foundation was developed. A theoretical solution of the dynamic stability analyses of the system is proposed. Based on the presented Formula (21) one can easily tune the natural frequencies of the system away from the wave frequency to avoid the instability of vertical

motion. In addition, the critical distance *eGB*,*critical* of pitching stability is investigated. Several trends about the stability due to wave excitation are obtained as follows:


In subsequent work the influences of nonlinear wave and current interaction will be investigated.

**Author Contributions:** Conceptualization, S.-M.L., Y.-Y.C., H.-C.H. and M.-S.L.; methodology, H.-C.H., S.-M.L. and Y.-Y.C.; software, H.-C.H. and M.-S.L.; validation, S.-M.L.; formal analysis, H.-C.H., S.-M.L.; investigation, H.-C.H.; resources, H.-C.H.; data curation, S.-M.L. and M.-S.L.; writing—original draft preparation, H.-C.H., S.-M.L. and M.-S.L.; writing—review and editing, H.-C.H.; visualization, S.-M.L.; supervision, Y.-Y.C.; project administration, H.-C.H.; funding acquisition, H.-C.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded bythe Ministry of Science and Technology of Taiwan, R. O. C., gran<sup>t</sup> number MOST 107-2221-E-110-077-MY3.

**Acknowledgments:** The authors would like to thank the referees for their helpful comments and suggestions. The support of the Ministry of Science and Technology of Taiwan, R. O. C., is gratefully acknowledged (MOST106- 3113-E-110-001-CC2, MOST107- 2221-E-110-077-MY3).

**Conflicts of Interest:** The authors declare no conflict of interest.
