Ship Towed by Kite: Investigation of the Dynamic Coupling

Abstract

This paper presents a series of dynamic simulations for a ship towed by kite. To ensure time efficient computations, seakeeping analysis with forward speed correction factors is carried out in the frequency domain and then transformed in the time domain by convolution. The seakeeping modeling is coupled with a zero-mass kite modeling assuming linear dependence of aerodynamic characteristics with respect to turning rate. Decoupled (segregated) and coupled (monolithic) approaches are assessed and compared in different environmental conditions. Results show that in regular beam waves, strong interactions between the kite and the ship motions are captured by the monolithic approach. Around the wave frequency, especially for the lower one tested (0.4 rad/s), a kite lock-in phenomenon is revealed. It is concluded that the mean kite towing force can be increased whereas the ship roll amplitude can even be decreased compared to a non-kite assisted ship propulsion configuration.

1. Introduction

This work is part of the beyond the sea^® research program launched by Yves Parlier and managed in partnership with French Ministry of Defence with the support of French Environment and Energy Management Agency (ADEME). The project attempts to develop a tethered kite system as an auxiliary device for the propulsion of merchant ships. The existing knowledge on ships towed by kite has demonstrated great prospects for this technology in terms of CO₂ emissions and fuel savings. However, studies on the limitations of this concept to guarantee the safety and the integrity of the ship are very limited.

Considering the mean kite towing force, Leloup et al. [1] and Naaijen et al. [2] solved the horizontal balance equations of a ship towed by kite to determine the fuel savings. Ran et al. [3] studied the contribution of a kite to the mean ship thrust, drift angle and rudder angle. All these previous studies neglected the interactions between the kite and the ship. The kite force is imposed as a predefined external force to the ship. The motions of such a system are highly dynamic since a kite experiences a periodic dynamic flight. In Bigi et al. [4], the influence of the kite attachment point on the deck is investigated on a fishing vessel equipped with a kite. This study is limited to horizontal ship motions. A maneuvering modeling is implemented with a monolithic coupling approach between the ship and the kite. The water is supposed to be calm and Bigi et al. [4] did not take into account the effect of the radiated waves on the ship motions. Thus, the influence of the kite excitation frequency on the added mass and damping of the ship is neglected. Since the hydrodynamic added mass and damping may depend strongly on the frequency of the motion [5,6,7,8], this assumption is questionable.

A ship sailing in waves is commonly studied through seakeeping codes based on the potential flow theory under the assumption of linear response of the ship to a given perturbation on a mean path. These studies are usually carried out in the frequency domain [9,10]. However, since the kite and the ship may be strongly coupled, their interactions cannot be directly computed through a spectral description of the kite excitation. Such computations are therefore limited to weak coupling. To perform a strong coupling between the kite and the ship, a time domain formulation is required. The aim of this paper is to assess the importance of taking into account the coupling between the kite and the ship motions by a time domain method.

As highlighted by Skejic [11], fast time-domain methods able to compute the 6 dof combining horizontal and vertical motions of a ship are based on the concept of linear convolution [12,13,14]. Such methods are first introduced by Cummins [15] in 1962. Böttcher [12] developed unified theories that couple the linear manoeuvring and seakeeping. In specific, he extended the method to simulate large amplitude motions taking into account nonlinear effects with Froude–Krilov effort and nonlinear roll damping for instance. To enable a fast computation of the linear convolution product, the use of state-space systems is introduced [13,16,17]. The identification of the state-space systems is detailed in [18] for the zero forward speed case. In this paper, this method is extended to the forward speed case.

Kite modeling is formulated under the zero mass assumption. Behrel et al. [19] presents experimental results comparing zero mass and point mass modeling. Differences are about few percent. Actually, both modelings give the same results when their coefficients identification is consistent. If taking into account the mass can be important for control issue, this is clearly not the case for performance assessment. Studies on kite performance assessment as ship propulsion device are based on zero-mass assumption (e.g., Wellicome and Wilkinson [20], Naaijen et al. [2], Dadd et al. [21] and Leloup et al. [1]). In addition, a linear evolution of the kite aerodynamics with the radius of trajectory curvature is proposed in this paper.

Section 2 presents the considered modeling of a ship towed by kite. Section 3 presents the results. First, the kite characteristics identification is detailed, then Section 3.2 defines the case study on the surface vessel combatant DTMB 5512. Section 3.3 and Section 3.4 investigate respectively through a calm water case and a regular beam wave case, the coupling between ship and kite. A general discussion on the methods and results is provided in Section 4 before to conclude.

2. Ship Towed by Kite Modeling

2.1. Reference Frame and Parametrization

The Earth fixed frame $\mathbf{n}$, the ship fixed frame $\mathbf{s}$ and the hydrodynamic frame $\mathbf{h}$ are sketched in Figure 1. The earth fixed frame is centered on $O_n$ with ${\underline{z}}_n$ pointing downward. The ship fixed frame is composed of ${\underline{x}}_s$ pointing forward in the ship symmetry plane, ${\underline{z}}_s$ pointing downward and ${\underline{y}}_s$ completing the direct orthogonal basis. ${\underline{z}}_s$ is normal to free surface when the ship is at the equilibrium. The origin of $\mathbf{s}$ denoted by $O_s$ is in the ship symmetry plane at a longitudinal distance $l_x$ from the aft perpendicular and at a vertical distance $l_z$ from the baseline. For a point P, the expression of its coordinates in a frame $\mathbf{f}$ is ${\underline{P}}^{\left(\mathbf{f}\right)}={\left[p_x^{\left(\mathbf{f}\right)},p_y^{\left(\mathbf{f}\right)},p_z^{\left(\mathbf{f}\right)}\right]}^T$ or ${\left[p_x,p_y,p_z\right]}_{\mathbf{f}}^T$. The vectors of a frame $\mathbf{f}$ are ${\underline{x}}_f$, ${\underline{y}}_f$ and ${\underline{z}}_f$ and the origin $\mathbf{f}$ is denoted $O_f$. The generalized position vector of the ship denoted by $\underline{S}={\left[s_x^{\left(\mathbf{n}\right)},s_y^{\left(\mathbf{n}\right)},s_z^{\left(\mathbf{n}\right)},{\varphi }_s^{\left(\mathbf{n}\right)},{\theta }_s^{\left(\mathbf{n}\right)},{\psi }_s^{\left(\mathbf{n}\right)}\right]}^T$ is the assembly of the position of $O_s$ and the ship Euler’s angles with respect to the $\mathbf{n}$ frame. The generalized ship velocity at $O_s$ expressed in $\mathbf{s}$ with respect to the $\mathbf{n}$ frame is denoted by ${\underline{V}}_s={\left[u_s,v_s,w_s,p_s,q_s,r_s\right]}_{\mathbf{s}}^T$, where the first three components are the linear velocities and the last three components are the turning rates. The transformation of a vector expressed in the $\mathbf{s}$ frame denoted by ${\underline{n}}^{\left(\mathbf{s}\right)}$ can be expressed in the $\mathbf{n}$ frame with ${\underline{n}}^{\left(\mathbf{n}\right)}={\underset{=}{T}}_{\mathbf{s}}^{\mathbf{n}}{\underline{n}}^{\left(\mathbf{s}\right)}$, where ${\underset{=}{T}}_{\mathbf{s}}^{\mathbf{n}}$ is the direct cosine matrix expressed as:

$${\underset{=}{T}}_{\mathbf{s}}^{\mathbf{n}}=\left[\begin{array}{ccc}\mathrm{c}{\psi }_s\mathrm{c}{\theta }_s& -\mathrm{s}{\psi }_s\mathrm{c}{\varphi }_s+\mathrm{c}{\psi }_s\mathrm{s}{\theta }_s\mathrm{s}{\varphi }_s& \mathrm{s}{\psi }_s\mathrm{s}{\varphi }_s+\mathrm{c}{\psi }_s\mathrm{c}{\varphi }_s\mathrm{s}{\theta }_s\\ \mathrm{s}{\psi }_s\mathrm{c}{\theta }_s& \mathrm{c}{\psi }_s\mathrm{c}{\varphi }_s+\mathrm{s}{\varphi }_s\mathrm{s}{\theta }_s\mathrm{s}{\psi }_s& -\mathrm{c}{\psi }_s\mathrm{s}{\varphi }_s+\mathrm{s}{\theta }_s\mathrm{s}{\psi }_s\mathrm{c}{\varphi }_s\\ -\mathrm{s}{\theta }_s& \mathrm{c}{\theta }_s\mathrm{s}{\varphi }_s& \mathrm{c}{\theta }_s\mathrm{c}{\varphi }_s\end{array}\right]$$

(1)

where, $\mathrm{c}$ and $\mathrm{s}$ denote respectively the cosine and the sinus functions. The turning rates and the time derivatives of the ship Euler’s angles satisfy the following relationship:

$$\left[\begin{array}{c}{p}_s\\ q_s\\ r_s\end{array}\right]={\underset{=}{R}}_{\mathbf{n}}^{\mathbf{s}}\left[\begin{array}{c}{\dot{\varphi }}_s\\ {\dot{\theta }}_s\\ {\dot{\psi }}_s\end{array}\right],$$

(2)

where,

$${\underset{=}{R}}_{\mathbf{n}}^{\mathbf{s}}=\left[\begin{array}{ccc}1& 0& -\mathrm{s}{\theta }_s\\ 0& \mathrm{c}{\varphi }_s& \mathrm{c}{\theta }_s\mathrm{s}{\varphi }_s\\ 0& -\mathrm{s}{\varphi }_s& \mathrm{c}{\varphi }_s\mathrm{c}{\theta }_s\end{array}\right].$$

(3)

The hydrodynamic frame is used by the seakeeping theory. It can be considered as Galilean since $\mathbf{h}$ moves at the constant mean ship forward speed $U_h$ on a straight path with respect to the earth fixed frame. The generalized position vector of the ship expressed in the $\mathbf{h}$ frame is denoted by $\underline{\xi }={\left[{\xi }_1,{\xi }_2,{\xi }_3,{\xi }_4,{\xi }_5,{\xi }_6\right]}_{\mathbf{h}}^T$, where the first three components are the distance between $O_h$ and a point H fixed to the ship. The three last components are the Euler’s angle of the ship with respect to the $\mathbf{h}$ frame. It is assumed that H is in the symmetry plane of the ship and $\underline{O_{s}H}={\left[d_x,0,d_z\right]}_{\mathbf{s}}^T$ with respect to the $\mathbf{s}$ frame.

Figure 2 illustrates the notations used for the kite modeling. The relative wind coordinates system moves at the tether attachment point velocity ${\underline{U}}_A$ with respect to the Earth fixed frame. The relative wind velocity defined in Equation (4) is the difference between the true wind velocity and the velocity of the tether attachment point.

$${\underline{U}}_{rw}={\underline{U}}_{tw}-{\underline{U}}_A$$

(4)

${\underline{x}}_{rw}={\underline{U}}_{rw}/∥{\underline{U}}_{rw}∥$ and ${\underline{y}}_{rw}$ is defined as: ${\underline{y}}_{rw}={\underline{z}}_n\times {\underline{x}}_{rw}/∥{\underline{z}}_n\times {\underline{x}}_{rw}∥$. In order to define a direct orthogonal coordinates system, ${\underline{z}}_{rw}$ is defined as follows: ${\underline{z}}_{rw}={\underline{x}}_{rw}\times {\underline{y}}_{rw}$. It should be noticed that the relative wind frame depends on the considered altitude due to the wind gradient. Consequently, the relative wind frame is defined with respect to a given point P and is thus denoted $\mathbf{rw}\left(P\right)$. For instance, the relative wind frame drawn in Figure 2 is $\mathbf{rw}\left(K\right)$. At the kite location, the relative wind frame is obtained from:

$${\underline{U}}_{rw\left(K\right)}={\underline{U}}_{tw\left(K\right)}-{\underline{U}}_A$$

(5)

The apparent kite wind velocity is denoted by ${\underline{U}}_{aw}$ and can be expressed as follows:

$${\underline{U}}_{aw}={\underline{U}}_{\mathbf{rw}\left(K\right)}-{\underline{U}}_k$$

(6)

where ${\underline{U}}_k$ denotes the kite velocity with respect to the $\mathbf{rw}\left(K\right)$ frame; the kite reference frame centered on K is defined as follows: ${\underline{z}}_k=\underline{AK}/∥\underline{AK}∥$, ${\underline{y}}_k={\underline{z}}_k\times {\underline{z}}_{rw}/∥{\underline{z}}_k\times {\underline{z}}_{rw}∥$ and ${\underline{x}}_k={\underline{y}}_k\times {\underline{z}}_k$; the kite velocity is directed by ${\underline{x}}_{vk}$. The kite elevation angle is given by:

$${\theta }_k=\frac{\pi }{2}-arccos\left({\underline{z}}_k·{\underline{z}}_{rw}\right)$$

(7)

The kite azimuth angle is defined as follows:

$${\varphi }_k=arccos\left({\underline{y}}_k·{\underline{y}}_{rw}\right)$$

(8)

2.2. Ship Modeling

As introduced in Section 1, a time domain formulation is a priori mandatory to represent the complete interaction between kite and ship motions. For the ship motion’s modeling part, we implemented the very classical and efficient method described by Fossen [22]. Appendix A gives details on the state space technique used. Appendix A.1 describes the equations of motion for a ship in a seaway. Appendix A.2 shows how hydrodynamics are computed using strip method theory and how is carried out the computation of linear convolution terms. Appendix A.3 provides details on the calculation of wave forces.

2.2.1. Ship Motions

The equations describing the motion of the system can be transformed into a system of first order differential equations as in Equation (9). This system is obtained with the transformation of the generalized ship velocity vector from the $\mathbf{s}$ frame to the $\mathbf{n}$ frame and Equations (A9) and (A12).

$$\left\{\begin{array}{cc}\dot{\underline{S}}\hfill & =\left[\begin{array}{cc}{\underset{=}{T}}_{\mathbf{s}}^{\mathbf{n}}& \underset{=}{0}\\ \underset{=}{0}& {\underset{=}{R}}_{\mathbf{s}}^{\mathbf{n}}\end{array}\right]{\underline{V}}_s\hfill \\ \dot{\underline{V}}\hfill & ={\left[{\underset{=}{M}}_S+\underset{=}{\tilde{A}}\right]}^{-1}\left[\underline{F}-\left[\underset{=}{\tilde{B}}+{\underset{=}{B}}_{\varphi }+\underset{=}{D}\right]{\underline{V}}_s-\underline{\mu }-\underset{=}{C}\underline{S}\right]\hfill \\ {\underline{\dot{y}}}_{ij}\hfill & ={\underset{=}{A}}_{ij}^{{}^{\prime }}{\underline{y}}_{ij}+{\underline{B}}_{ij}^{{}^{\prime }}\delta V_{s,j},\forall i,j\in ⟦1;6⟧\hfill \end{array}\right.$$

(9)

Equation (9) represents 12 scalar equations for the ship and 75 scalars equations for the state-space systems. Assuming, for instance, that the order of each state space system is 5 and taking into account the ship symmetry. Thus, with the present modeling, the ship is described by 87 scalar equations depending on the state space modeling orders. This system of differential equations is numerically integrated with the 4th order Runge–Kutta scheme with a fixed time step.

2.2.2. Ship Modeling Validation

The present ship modeling was compared to Experimental Fluid Dynamic (EFD) data and to the STF strip theory results on the David Taylor Model Basin (DTMB) 5512. The model tested is a 1:46.6 scale model. The hull form and its characteristics are respectively plotted in Figure 3 and summarized in

Table 1. DTMB 5512 hull and full scale characteristics.

Parameter	Units	5512	Full Scale
Scale ratio	-	46.6	1
Length, $L_{pp}$	m	3.048	142.04
Beam, B	m	0.405	18.87
Draft, T	m	0.132	6.15
Weight	Kg - t	86.6	8763.5
LCG	m	1.536	71.58
VCG	m	0.162	7.55
Pitch radius of gyration, $k_5$	m	0.764	35.6

. The experimental data are provided by the University of Iowa [23] and are presented in Irvine et al. [24]. The EFD data concern the heave and pitch motions in regular head waves, with and without forward speed.

The computed ship motions were performed at zero forward speed and at a Froude number of 0.28, which corresponds to $U=1.53$ m.s⁻¹ and with frequency head waves from 1 rad.s⁻¹ to 7.5 rad.s⁻¹. $O_s$ was defined at the ship center of gravity, hence $l_x=LCG$ and $l_z=VCG$. Figure 4 and Figure 5 plot the heave and pitch transfer function obtained with the experimental data, with the STF strip theory and with the present modeling. The experimental data are obtained for different wave steepness $s_w=\left\{0.025,0.05,0.075\right\}$. The amplitude of the transfer function for heave motion (a) is directly the ratio of the heave amplitude motion to the wave amplitude. The amplitude of the transfer function for the pitch motion (c) is given by the ratio of the motion amplitude (in radians) to the wave amplitude multiplied by the wave number k. The phase angle of the present modeling is obtained by cross correlation between the free surface elevation and the ship motion time series.

Concerning the amplitude, an overall good agreement is found with (Figure 5) and without (Figure 4) forward speed between the EFD data, the STF strip theory, the time domain and the frequency domain modeling. For the considered waves, the influence of the wave steepness on the EFD data is not significant. As shown in Figure 4 and Figure 5, the STF strip theory and the time domain approach match with a very good accuracy for the amplitude. Very small differences can be observed in terms of phase angle, but these differences are probably caused by the post-processing method (Figure 4 and Figure 5). The very small differences with the STF strip theory are due to the approximations performed with the identification method of the transfer functions $H_{ij}$. As a conclusion for the heave and pitch motions, the transformation of the equation of motion into the ship reference frame and the state-space modeling identification method are consistent and accurate enough.

The roll motion is much more difficult to predict with a linear modeling. Full-scale results are presented in Figure 6. Results by Ikeda et al. [25] are considered as a reference. In the present method, the extra roll damping was set to $1.53·{10}^8$ kg.m.s⁻¹ in order to obtain the same roll amplitude at the natural frequency. As shown in Figure 6, for frequencies greater than 1 rad.s⁻¹, the roll motion amplitude is well predicted despite a slight difference in phase. Nevertheless predictions appear to be less accurate for the lower frequencies, except for the natural roll frequency, ${\omega }_{roll}=0.56$ rad.s⁻¹, where, by identification, both modelings match.

2.3. Kite Modeling

2.3.1. Short Literature Survey

Kites are mostly studied as a power generation system since 1980 [26]. Comparing the weight of the kite on the one hand and the traction provided on the other, Loyd considers that it is relevant to neglect the mass of the kite. Thus, he developed the so-called zero-mass modeling. Taking the skysails system, the most advanced industrial example, gives us a good illustration of the relevance of this assumption. Paulig et al. [27] computes a pulling force to mass ratio of 50daN/kg in the case of a 320 m² kite with 360 m tether. Wellicome extended the zero-mass modeling defining a figure-8 trajectory for the dynamic mode [20,28]. Argatov et al. [29] proposed more recently a simpler trajectory definition. In addition, Argatov definition allows a much faster calculation and its formulation avoids discontinuities in kite acceleration. The size of the trajectory can easily be modified thanks only to azimuth and elevation amplitude parameters. In 2006, Naaijen et al. [2] computes the fuel savings on a full range of true wind angles for the British Bombardier, a 50000 DWT tanker. In 2011, Dadd et al. [21] improved the calculations with more realistic trajectories using the Wellicome formalism. In 2014, Leloup et al. [1,30] proposed a fully analytical solution for the zero-mass modeling that allows faster and more reliable computations avoiding optimization process failures. Leloup takes into account adequately the wind gradient effect and takes advantage of the Argatov trajectory definition for even faster computations.

2.3.2. Zero-Mass Kite Model

According to the zero-mass kite modeling, the masses of the tether and the kite are neglected and the tether is assumed to be straight and of constant length $L_t$. Consequently, the kite velocity with respect to the $\mathbf{rw}\left(K\right)$ frame is normal to ${\underline{z}}_{k0}$ and for any configuration the tether tension is opposed to the aerodynamic kite force. Assuming that the kite flies with a constant lift to drag ratio and that the apparent wind velocity is in its symmetry plane, Leloup et al. [1] expressed the kite velocity according to a given kite velocity direction ${\underline{x}}_{vk}$ in Equation (10):

$${\underline{U}}_k=U_{rw}\left[\left({\underline{x}}_{vk}·{\underline{x}}_{rw}\right)+\sqrt{{\left({\underline{x}}_{vk}·{\underline{x}}_{rw}\right)}^2+{\left(\frac{{\underline{z}}_{k0}·{\underline{x}}_{rw}}{sin{ϵ}_k}\right)}^2-1}\right]{\underline{x}}_{vk}$$

(10)

The kite velocity is a real number if the following condition is satisfied:

$${\left({\underline{x}}_{vk}·{\underline{x}}_{rw}\right)}^2+{\left(\frac{{\underline{z}}_{k0}·{\underline{x}}_{rw}}{sin{ϵ}_k}\right)}^2-1⩾0.$$

(11)

Then, the tether tension is given by the following formula:

$${\underline{T}}_k=-\frac{C_{lk}{\rho }_a{A}_k{U}_{aw}^2}{2cos{ϵ}_k}{\underline{z}}_k$$

(12)

The generalized tether force vector acting on the ship at $O_s$ with respect to $\mathbf{s}$ frame is expressed as follows:

$${\underline{F}}_k={\left[\begin{array}{cc}{\underline{T}}_k^{\left(\mathbf{s}\right)}& {\underline{O_{s}A}}^{\left(\mathbf{s}\right)}\times {\underline{T}}_k^{\left(\mathbf{s}\right)}\end{array}\right]}^T$$

(13)

In order to represent the wind friction with the sea, the true wind velocity ${\underline{U}}_{tw}$ is function of the altitude according to the wind gradient law recommended by the ITTC [31].

$${\underline{U}}_{tw}\left(k_z^{\left(\mathbf{n}\right)}\right)={\underline{U}}_{10}{\left(\frac{k_z^{\left(\mathbf{n}\right)}}{10}\right)}^{1/7}$$

(14)

2.3.3. Kite Control According to a Trajectory

The kite velocity direction ${\underline{x}}_{vk}$ is controlled in order to follow a trajectory denoted by $\mathcal{C}$ (Figure 7). The kite velocity direction is defined by the target point $\tilde{K}$ expressed as follows:

$$\tilde{K}=\mathcal{C}\left(\lambda +∥{\underline{U}}_k∥dt\right),$$

(15)

where $\lambda$ is the curvilinear abscissa of $C_{\lambda }$ the closest point of the trajectory from the current kite position K. Hence, the kite velocity direction ${\underline{x}}_{vk}$ is defined as follows:

$${\underline{x}}_{vk}=\frac{\left(\underline{K\tilde{K}}·\underline{x_k}\right){\underline{x}}_k+\left(\underline{K\tilde{K}}·\underline{y_k}\right){\underline{y}}_k}{∥\left(\underline{K\tilde{K}}·\underline{x_k}\right){\underline{x}}_k+\left(\underline{K\tilde{K}}·\underline{y_k}\right){\underline{y}}_k∥}$$

(16)

The point $C_{\lambda }$ is determined according to the following equation:

$$\underline{K{C}_{\lambda }}·{\underline{t}}_{\lambda }=0,$$

(17)

where ${\underline{t}}_{\lambda }$ is the tangent vector to the trajectory at curvilinear abscissa $\lambda$. Equation (17) is solved numerically with a Newton–Raphson algorithm.

2.3.4. Theoretical Lemniscate Trajectory

To perform a dynamic flight, the considered theoretical trajectory is the eight shaped Lissajous curve as used by Argatov et al. [29] and Leloup [32]. This type of trajectory has the advantage to avoid twists of the tether system. The Lissajous trajectories are defined in terms of elevation ${\theta }_C$ and azimuth ${\varphi }_C$:

$$\left\{\begin{array}{cc}{\theta }_C\hfill & =\Delta {\theta }_{8}sin\left(2\alpha \right)+{\theta }_8\hfill \\ {\varphi }_C\hfill & =\Delta {\varphi }_{8}sin\left(\alpha \right)+{\varphi }_8\hfill \end{array}\right.,$$

(18)

where, $\alpha \in \left[0;2\pi \right]$. ${\theta }_8$ and ${\varphi }_8$ are the elevation and the azimuth of the center of the trajectory denoted by $C_8$. The trajectory is defined with respect to $\mathbf{rw}\left(z_{ref}\right)$, the relative wind frame with respect to the wind measurement altitude. On a sphere of radius $L_t$, a point C of the trajectory is defined in terms of elevation ${\theta }_C$ and azimuth ${\varphi }_C$ with respect to $\mathbf{rw}\left(z_{ref}\right)$ frame as follows:

$$\underline{C}={\left[\begin{array}{c}{L}_{t}cos{\theta }_{C}cos{\varphi }_C\\ L_{t}cos{\theta }_{C}sin{\varphi }_C\\ -L_{t}sin{\theta }_C\end{array}\right]}_{\mathbf{rw}\left(z_{ref}\right)},$$

(19)

As shown in Figure 8, the lemniscate trajectories can be rotated by an angle ${\chi }_8$ around $\underline{C_{8}A}$. Any point of the trajectory must satisfy Equation (11) insuring the realness of the kite velocity. The benefit to defining the trajectory with respect to the relative wind frame is that the realness of the kite velocity along the trajectory does not depend on the ship motions. Nevertheless, the ship motions may modify significantly the shape of the trajectory with respect to the Earth fixed frame $\mathbf{n}$. Vertical motions are especially known to cause unwanted overloads and flight instabilities [27]. Consequently, two trajectory definitions are investigated further in the paper. The first definition takes into account all components of the tether attachment point velocity. The second definition takes only the horizontal components of the tether attachment point velocity with respect to the earth fixed frame $\mathbf{n}$ to compute the relative wind frame. This modified relative wind frame of trajectory definition is denoted by $\widetilde{\mathbf{rw}}\left(z_{ref}\right)$ and is defined from:

$${\underline{U}}_{rw}^{\left(\mathbf{n}\right)}={\underline{U}}_{tw}^{\left(\mathbf{n}\right)}\left(z_{ref}\right)-{\left[u_a^{\left(\mathbf{n}\right)},v_a^{\left(\mathbf{n}\right)},0\right]}^T$$

(20)

2.4. Two Coupling Methods

The tether tension induces motions to the ship. The expression of generalized kite force acting on the ship is expressed as follows:

$${\underline{F}}_k={\left[\begin{array}{cc}{\underline{T}}_k^{\left(\mathbf{s}\right)}& {\underline{O_{S}A}}^{\left(\mathbf{s}\right)}\times {\underline{T}}_k^{\left(\mathbf{s}\right)}\end{array}\right]}^T$$

(21)

In addition, according to the zero-mass kite modeling (cf. Section 2.3.2 and Equation (5)), the ship motions can modify the kite flight and the tether tension through the relative wind speed at the kite altitude with respect to the tether attachment point A:

$${\underline{U}}_{rw}\left(k_z^{\left(\mathbf{n}\right)}\right)={\underline{U}}_{tw}\left(k_z^{\left(\mathbf{n}\right)}\right)-{\underline{U}}_A$$

(22)

Two coupling approaches are investigated in this paper: a monolithic and a segregated approaches. The monolithic approach takes into account all the coupling terms between the kite and the ship modeling. Here, the considered segregated approach assumes a predefined kite force and then solves the ship equations of motion separately.

2.4.1. A Monolithic Approach

The monolithic equations of motion of a ship towed by kite are obtained with Equations (9) and (10) as follows:

$$\left\{\begin{array}{cc}\dot{\underline{S}}\hfill & =\left[\begin{array}{cc}{\underset{=}{T}}_{\mathbf{s}}^{\mathbf{n}}& \underset{=}{0}\\ \underset{=}{0}& {\underset{=}{R}}_{\mathbf{s}}^{\mathbf{n}}\end{array}\right]{\underline{V}}_s\hfill \\ {\dot{\underline{V}}}_s\hfill & ={\left[{\underset{=}{M}}_S+\underset{=}{\tilde{A}}\right]}^{-1}\left[\underline{F}-\left[\underset{=}{\tilde{B}}+{\underset{=}{B}}_{\varphi }+\underset{=}{D}\right]{\underline{V}}_s-\underline{\mu }-\underset{=}{C}\underline{S}\right]\hfill \\ {\underline{\dot{y}}}_{ij}\hfill & ={\underset{=}{A}}_{ij}^{{}^{\prime }}{\underline{y}}_{ij}+{\underline{B}}_{ij}^{{}^{\prime }}\delta V_{s,j},\forall i,j\in ⟦1;6⟧\hfill \\ {\dot{\underline{K}}}^{\left(\mathbf{n}\right)}\hfill & ={\underline{U}}_A^{\left(\mathbf{n}\right)}+U_{rw}\sqrt{{\left({\underline{x}}_{vk}·{\underline{x}}_{rw}\right)}^2+{\left(\frac{{\underline{z}}_{k0}·{\underline{x}}_{rw}}{sin{ϵ}_k}\right)}^2-1}{\underline{x}}_{vk}^{\left(\mathbf{n}\right)}\hfill \end{array}\right.$$

(23)

Equation (23) adds 3 scalar equations for the kite. With the monolithic approach, the fully coupled system between the ship motions and the kite motions is solved. This monolithic system of differential equations is numerically integrated with the 4th order Runge–Kutta scheme with fixed time step.

2.4.2. A Segregated Approach

By contrast to the monolithic approach, the segregated approach considers only the mean tether attachment point velocity on the ship. The ship motions are computed applying the time series of the kite towing force as an external forces. Thus, the ship equation of motion can be expressed as follows:

$$\left\{\begin{array}{cc}\dot{\underline{S}}\hfill & =\left[\begin{array}{cc}{\underset{=}{T}}_{\mathbf{s}}^{\mathbf{c}}& \underset{=}{0}\\ \underset{=}{0}& {\underset{=}{R}}_{\mathbf{s}}^{\mathbf{c}}\end{array}\right]{\underline{V}}_s\hfill \\ {\dot{\underline{V}}}_s\hfill & ={\left({\underset{=}{M}}_S+\underset{=}{\tilde{A}}\right)}^{-1}\left[{\underline{F}}_k\left(t\right)+{\underline{F}}^{{}^{\prime }}-\left(\underset{=}{\tilde{B}}+{\underset{=}{B}}_{\varphi }+\underset{=}{D}\right){\underline{V}}_s-\underline{\mu }-\underset{=}{C}\underline{S}\right]\hfill \\ {\underline{\dot{y}}}_{ij}\hfill & ={\underset{=}{A}}_{ij}^{{}^{\prime }}{\underline{y}}_{ij}+{\underline{B}}_{ij}^{{}^{\prime }}\delta V_{s,j},\forall i,j\in ⟦1;6⟧\hfill \end{array}\right.$$

(24)

where ${\underline{F}}^{{}^{\prime }}$ denotes the external forces such as rudder, propeller and windage forces excluding the kite forces applied as a time series ${\underline{F}}_k\left(t\right)$. This segregated approach could be very practical to study the motions of a ship towed by kite. Even if here, this approach is performed into the time domain, the segregated approach can be performed into the frequency domain applying the kite excitation spectrum directly in Equation (A1).

3. Results

3.1. Validation and Comments on the Kite Aerodynamic Characteristics

On the basis of the experimental data obtained by Behrel et al. [19,33], the zero-mass modeling is compared to the experimental data. The zero mass kite modeling is dependent of two parameters, the kite lift coefficient $C_{lk}$ and the glide angle ${ϵ}_k$. These coefficients must be adapted in order to fit the data. The onshore full scale trials [19,33] were performed with a classical kite Cabrinha^® Switchblade2016 of 5 m² designed for kite-surfing. The tether length was 80 m long. During the run, the kite performed eight shaped trajectory controlled by an autopilot based on the algorithm proposed in [34,35]. The experimental kite position is determined with a 3D load cell assuming that the tethers are straight. The evolution of the wind velocity with the altitude was identified thanks to a SOnic Detection And Ranging (SODAR). The experimental results presented here correspond to a phase averaging post-processing of a 5-minute kite flight run. For the presented case, the wind velocity was interpolated from the SODAR data with the following linear function:

$$U_{tw}=3.16-0.035k_z^{\left(\mathbf{n}\right)}$$

(25)

The wind gradient function in Equation (25) is used for the following presented simulation results.

The zero-mass modeling presented in [1] assumes a constant glide angle and a constant lift coefficient. The kite velocity depends only on its direction and position in the wind window and on the glide angle ${ϵ}_k$. Consequently, the present zero-mass modeling results are obtained with a glide angle enabling that the simulated kite trajectory has the same period than the experimental data, which is 5.86 s. The kite towing force is dependent of the kite lift coefficient. Hence, the kite lift coefficient is adjusted in order to obtain the same mean towing force. The kite trajectory used for the simulation is the same than the kite trajectory of the experimental measurements. The position of the kite is integrated with the 4th order Runge–Kutta scheme with a time step of 0.1 s.

Figure 9 shows the evolution of the magnitude of the kite velocity (a) and the evolution of the magnitude of the tether tension at A on the average trajectory (b). The experimental velocity is obtained by finite differentiation, which induces some noise in the signal sampled at 50 Hz. This noise in terms of speed and tension for the simulated results is still present since the simulation uses the average kite trajectory of the experimental data. However, this noise is barely visible since the sampling frequency of the simulation is only 10 Hz. The kite glide angle and the lift coefficient used for the zero-mass modeling are respectively 12.45° and 0.855. The main difference between the experimental data and the results concerns the tether tension and the kite velocity amplitude. The simulated amplitude is slightly underestimated for the velocity and is particularly underestimated for the tension. The tether tension given by Equation (12) is a quadratic function of the apparent wind velocity. Therefore, for a relative kite velocity error, the relative tension error is twice. In order to perform an eight trajectory, the back tethers are steered. The difference between the back tethers lengths leads to important modifications of the kite flying shape. Hence, it can be expected a modification of the kite lift coefficient and of the glide angle leading to the significant differences observed between the zero-mass modeling with constant aerodynamic characteristics and the measurements. Leloup [36] and Duport [37] demonstrated the variation of the aerodynamic characteristics with respect to the turning rate. Even without considering tip vortices, the yaw rate induces an evolution of the local inflow velocity along the kite span. According to Equation (10), the kite velocity is independent of its area. Consequently, for a given radius of curvature of the trajectory, the evolution of the local inflow along the span and the effect of the flying shape deformation are more consequent with larger kites. Here it is assumed that the evolution of the kite glide angle and lift coefficient are proportional to the ratio of a characteristic length of the kite with the radius of curvature of the trajectory. For a given aspect ratio, the characteristic length of the kite can be $\sqrt{A_k}$. Denoting the trajectory radius of curvature $R_{\mathcal{C}}$, the glide angle and the lift coefficient could written as follows:

$$\left\{\begin{array}{cc}{ϵ}_k\hfill & ={ϵ}_0+{\kappa }_{ϵ}^{{}^{\prime }}\frac{\sqrt{A_k}}{R_{\mathcal{C}}}\hfill \\ C_{lk}\hfill & =C_{l0}+{\kappa }_l^{{}^{\prime }}\frac{\sqrt{A_k}}{R_{\mathcal{C}}}\hfill \end{array}\right.$$

(26)

where, ${\kappa }_{ϵ}^{{}^{\prime }}$ and ${\kappa }_l^{{}^{\prime }}$ are two coefficients. Equation (26) can be rewritten in terms of heading rate according to the relationship between the radius of curvature of the trajectory and the kite velocity. Furthermore, according to Equation (10), the kite velocity is proportional to the relative wind speed, thus the correction proposed can be rewritten as functions of the time derivative of the heading and the relative wind speed as follows:

$$\left\{\begin{array}{cc}{ϵ}_k\hfill & ={ϵ}_0+\frac{\sqrt{A_k}}{U_{rw}}{\kappa }_{ϵ}\left|\dot{\gamma }\right|\hfill \\ C_{lk}\hfill & =C_{l0}+\frac{\sqrt{A_k}}{U_{rw}}{\kappa }_l\left|\dot{\gamma }\right|\hfill \end{array}\right.$$

(27)

The determination of the coefficient ${ϵ}_0$, ${\kappa }_{ϵ}$, $C_{l0}$ and ${\kappa }_l$ can be evaluated by comparison with the experimental data. First, ${ϵ}_0$ and ${\kappa }_{ϵ}$ are identified in order to obtain respectively the same maximum and minimum speed than the experiments. Then, $C_{l0}$ and ${\kappa }_l$ are identified in order to obtain respectively the same maximum and minimum tether tension than the experiments. According to this method, the following coefficient values are identified:

$$\left\{\begin{array}{cc}{ϵ}_0\hfill & =0.2013\mathrm{rad}\hfill \\ {\kappa }_{ϵ}\hfill & =0.0422\mathrm{s}\hfill \\ C_{l0}\hfill & =0.9856\hfill \\ {\kappa }_l\hfill & =-0.3718{\mathrm{s}.\mathrm{rad}}^{-1}\hfill \end{array}\right.$$

(28)

These values are consistent with the ones by Behrel et al. [19], who did the identification with a quite different and more straightforward method. The results in terms of kite velocity and tether tension are plotted in Figure 9. With the modification of the glide angle and the lift coefficient, the noise in the kite velocity and the tether tension time series is more significant. Indeed, the computation of the kite yaw rate requires two finite differentiations of the kite position. As expected with the identification method of ${ϵ}_0$, ${\kappa }_{ϵ}$, $C_{l0}$ and ${\kappa }_l$, the amplitude of the kite velocity and the amplitude of the tether tension are respected.

An overall good agreement is found in terms of velocity but a slight phase difference is observable between the simulation results and the experimental data. The linear modifications of the glide angle and of the lift coefficient with the kite yaw rate lead to a better agreement with the experimental results than the zero mass modeling. However, these modifications should be confirmed with more experimental cases. Since at the first order the ship motions are proportional to the amplitude of the excitation forces, these modifications enhancing the prediction of the kite excitation force amplitude are crucial.

Moreover, the minimum allowable radius of curvature of the trajectory can be estimated according to Equation (26). Indeed, during a kite flight $C_{lk}$ must be positive. This means that the trajectory radius of curvature must satisfy the following condition:

$$R_{\mathcal{C}}⩾-\frac{{\kappa }_l^{{}^{\prime }}}{C_{l0}}\sqrt{A_k}$$

(29)

As a numerical application, the minimum trajectory radius of curvature of the 5 m² Cabrinha^® Switchblade is $2.23\sqrt{A_k}\approx 5.0$ m, Which is a little less than twice the projected wingspan. Given that the lift coefficient must have a strictly positive value in order to fly, we are here outside the validity domain of the proposed modeling. However, the order of magnitude is confirmed by Fagiano who evaluated the minimum turning radius at 2.5 times the wingspan [38,39]. The coefficients ${ϵ}_0$, ${\kappa }_{ϵ}$, $C_{l0}$ and ${\kappa }_l$ are dimensionless quantities. Consequently, the presented modifications in Equations (27) and (28) are retained as formulated for the rest of the paper.

3.2. Study Case

In order to simplify the analysis, only vertical ship motions (heave, roll and pitch) of the DTMB 5512 at full scale were considered here. The analysis was focused on the roll motion. Thus, in the scope to observe significant roll motion, a true wind angle ${\beta }_{tw}=90°$ was chosen. A true wind speed of reference $U_{10}=$10 m.s⁻¹ corresponding to the high range of a fresh breeze from the Beaufort scale was considered. The ship speed was set to $U_h=$7.5 m.s⁻¹ since it corresponds to a common sailing speed condition of the world merchant ship fleet. A kite with an area of $A_k=$500 m² and with the aerodynamic characteristics determined in Equation (28) was used. The tether attachment point was in the center plane, 7.9 m above the water line and 25 m in front of the center of gravity.

The kite flight trajectory corresponded to a Lissajous trajectory as defined in Section 2.3.4. The amplitudes of the trajectory were arbitrarily set to $\Delta {\varphi }_8=20°$ and $\Delta {\theta }_8=8°$. According to the formulation of the lift coefficient in Equation (26), the tether length condition to obtain a positive lift coefficient was $L_t⩾360$ m. The kite excitation was studied by varying the tether length. Tether length between 360 m and 1000 m were investigated. The center of the trajectory $[{\varphi }_8,{\theta }_8]$ and the angle of the trajectory ${\chi }_8$ around the axis $\underline{C_{8}A}$ were determined optimizing of the longitudinal kite towing force with the same code used by Leloup et al. [1]. In this configuration, the kite provided about 16% of the propulsive power, provided the data from [40].

Figure 10a–c show respectively the evolution of the ${\varphi }_8$, ${\theta }_8$ and ${\chi }_8$ with the tether length. A calm water case and three regular beam wave cases are considered. A regular beam wave of 2.5 m high consistent with a fresh breeze. Three wave frequencies were investigated $\left\{0.4,0.56,0.8\right\}$ rad.s⁻¹. The 0.56 rad.s⁻¹ wave frequency corresponded to the natural roll ship frequency. For all following results, the simulation time was 1640 s with a time step of 0.3 s.

3.3. Calm Water Case

3.3.1. Kite Excitation Spectrum

Figure 11a,b show respectively the time series and the spectrum of the kite roll excitation moment obtained with the segregated approach for a tether length $L_t=500$ m. Only the varying part of the kite excitation is taken into account to compute the kite excitation spectrum. With the segregated approach, the kite flight is not modified by the ship motions. It can be noticed in Figure 11 that the roll excitation moment is mainly composed of several harmonics. For convenience, the harmonics are denoted ${\omega }_{ki}$ where i is a natural number. Only the first, the second and the fourth harmonics are significant. The whole spectrum is not represented but the other harmonics that occur at higher frequencies are not significant. The second and the fourth harmonics are the most powerful. The second harmonic should appear to be the most critical for the ship motions due to its proximity with the natural roll ship frequency.

3.3.2. Comparison of the Segregated Approach with the Monolithic Approach

Figure 12 shows the evolution of the roll amplitude, first kite excitation harmonic frequency and amplitude kite moment with respect to the tether length. Three methods are compared: the segregated approach in dashed-dotted line and the two monolithic approaches with the kite trajectories defined in $\widetilde{\mathbf{rw}}\left(A\right)$ and in $\mathbf{rw}\left(A\right)$ respectively in solid and dashed lines.

As expected, the first kite excitation harmonic ${\omega }_{k1}$ decreases with the tether length since the angular amplitude of the Lissajous trajectories is remained constant. No major difference can be noticed between the three approaches in terms of harmonics frequencies. The roll amplitude $\Delta {\varphi }_s$ and the amplitude of the kite roll moment $\Delta K_k$ predicted by the segregated are higher than those predicted by the two monolithic approaches. The two trajectory definitions in $\widetilde{\mathbf{rw}}$ and $\mathbf{rw}$ give almost the same results both in terms of roll amplitude and kite roll moment amplitude. For all approaches presented, the roll amplitude evolution is similar. The kite moment amplitude increase quasi linearly with the tether length. Four ruptures (mainly three) can be observed in the evolution of the kite moment amplitude. These ruptures corresponds to the evolution of the trajectory parameters with the tether length as shown in Figure 10. Due to the wind gradient, the longer the tether is, the larger the kite roll moment is. This increase in kite roll moment explains the continuous increase in roll amplitude. However, two important increases can be noticed between $L_t=360$ m and $L_t=500$ m and from $L_t=790$ m, which are not explained by the kite roll moment. In fact, these two increases are due to the proximity between the two most powerful kite roll excitation harmonics and the vessel natural roll frequency. Indeed, for $L_t=440$ m the second harmonic frequency is almost equal to the natural roll frequency of the ship. From $L_t=790$ m, the fourth harmonic frequency approaches the natural roll ship frequency. In case of a slower true wind speed, the increase of the kite roll moment with the tether length would be less significant. Thus, it could be observed a maximum of roll amplitude for tether length corresponding to a match between a kite harmonic frequency and the natural roll frequency of the ship. According to these results in calm water, the segregated approach is slightly conservative with respect to both monolithic approaches.

3.4. Regular Beam Wave Case

To investigate a case closer to a real ocean environment, the influence of regular beam waves is studied. The wave considered is a 2.5 m high at the frequencies ${\omega }_w$ of 0.4 rad.s⁻¹, 0.56 rad.s⁻¹ and 0.8 rad.s⁻¹. As for the calm water case, the frequency domain of kite excitation is scanned with different tether lengths from $L_t=360$ m to $L_t=990$ m with a tether length step of 10 m.

3.4.1. Ship Vertical Motion Influence on Kite Trajectory

This section focuses on the influence of ship vertical motion on kite flight trajectory. Paulig et al. [27] underlined the major disturbances caused by this motion in terms of flight control and overload. Although this is in line with the rare feedback from experiences at sea, this has never been simulated before, all published studies being carried out so far in cases of calm water. Figure 13 shows the trajectory defined neglecting the ship vertical motion in solid line and the trajectory defined in the actual relative wind frame in dashed line. The latter has considerably reduced turning radius. It is not in line with the relationship on the minimum kite flight radius highlighted by Fagiano nor with the experiments by Behrel presented in Section 3.1. Such a sharp trajectory could not be achieved in the real world. As a consequence, in the rest of the paper, all the trajectories will be defined in $\widetilde{\mathbf{rw}}\left(A\right)$ and will be called.

3.4.2. Interactions with Regular Beam Waves

Figure 14, Figure 15 and Figure 16 show respectively for three wave frequencies, the roll amplitude of the ship, the first kite harmonic frequency and roll moment. For each wave cases, the monolithic approach is compared to the segregated approach. The dashed-dotted line corresponds to the roll amplitude of the ship due to the wave excitation without kite. For all wave cases, the roll amplitude and the amplitude of the kite roll moment increase globally with the tether length as for the calm water case. In contrary to the calm water case, the amplitude of the kite roll moment obtained with the monolithic approach is globally larger than the one obtained with the segregated approach. However, it should be noticed some drops in the evolutions of the roll amplitude and amplitude of the kite roll moment. For instance, in Figure 15a, for the tether length $L_t=440$ m, a significant drop of the roll amplitude can be noticed in the segregated case. This phenomenon of interaction between the kite and the ship can explained the fact that the monolithic approach predicts a smaller roll amplitude despite a larger kite roll moment amplitude compared to the segregated approach. Indeed, the secondary kite harmonic at the wave frequency has a relative phase with the ship motion. As an example, for the case plotted in Figure 17, the difference in phase angle between the ship motion and the kite excitation at $\omega =0.4$ rad.s⁻¹ is 60.7°. For the case with the frequency wave ${\omega }_w=0.8$ rad.s⁻¹, the interaction between the kite and the ship is less significant since the wave frequency is far from the most powerful kite harmonic frequencies. Moreover, in contrary to the calm water case, the segregated approach is not necessarily conservative with respect to the monolithic approach as shown by Figure 16a.

Concerning the evolution of the first harmonic frequency, differences can be noted between the two approaches. Indeed, for certain ranges of tether length, the harmonic frequency obtained with the monolithic approach remains constant. For instance, on Figure 15b, the first harmonic frequency remains constant to half the wave frequency for tether lengths within the range $\left[390;470\right]$ m. A sharp evolution of the roll amplitude concerns both approaches. In the monolithic case, this phenomenon does not occur only when a harmonic frequency of the kite corresponds to the wave frequency. With a wave frequency of 0.4 rad.s⁻¹ and a tether length $L_t=840$ m, a drop of the ship roll amplitude can be noticed. Figure 17 shows that no principal kite harmonic frequency corresponds to the wave frequency. However, due to coupling between kite and ship motions, secondary harmonics appear. The secondary harmonic frequencies are denoted ${\omega }_{ki}^{{}^{\prime }}$. This phenomenon occurs when the frequency gap between the closest principal harmonic frequency and the wave frequency is a submultiple of the first kite harmonic frequency. For the case presented in Figure 17, ${\omega }_{k1}=3{\omega }_{k1}^{{}^{\prime }}$. Furthermore, it should be noted that these drops in roll amplitude are due to the interactions between the kite and the ship. With the monolithic approach, these drops are independent from initial conditions. The most important drops in roll amplitude occur when ${\omega }_{k1}={\omega }_{k1}^{{}^{\prime }}$ or ${\omega }_{k1}=2{\omega }_{k1}^{{}^{\prime }}$. The importance of the phenomenon occurring at ${\omega }_{k1}=n{\omega }_{k1}^{{}^{\prime }}$ decreases with the increasing value of the integer n. Moreover, the importance of the phenomenon decreases when the interaction with the wave concerns high harmonic orders.

3.4.3. Kite Lock-in Phenomenon

As shown in the previous section, an important interaction phenomenon between the kite and the ship occured. We observed the appearance of plateaus which were all the more pronounced as they were located near the first harmonics of the moment of excitation produced by the kite. It can be observed in Figure 14a and Figure 15a respectively at $L_t=650$ m and $L_t=470$ m, that the drop in roll amplitude predicted by the monolithic approach is important enough to lead to smaller roll amplitude than the case without kite. A particular attention is paid to the case where the wave frequency ${\omega }_w=0.4$ rad.s⁻¹ and $L_t=470$ m, because the interaction between the kite and the ship is win-win. Indeed, for this case the mean kite towing force predicted by the monolithic approach is 8.0% more important than the mean kite towing force predicted with the segregated approach. In addition, the ship roll amplitude is 1.4% weaker than without kite.

Figure 18 shows the spectrum of the roll motion (a) and the spectrum of the kite excitation (b). It can be noticed that the second harmonic of the kite roll moment is attracted towards the wave frequency. Since the kite harmonic frequency is increased towards the wave frequency, the kite performed the whole trajectory faster. The time to perform the trajectory is decreased by 5.7%. This leads to a higher apparent wind speed and therefore to a higher average kite towing force. Moreover, the difference in phase angle between the roll motion and the kite roll moment is 18.2°, leading to a reduction of the roll amplitude. This attraction towards the wave frequency can be noticed also on the secondary harmonic corresponding to the case ${\omega }_{k1}=2{\omega }_{k1}^{{}^{\prime }}$. However, the effects are less significant. This phenomenon is similar to the lock-in phenomenon happening with the vortex-induced vibrations [41,42]. Henceforth, this phenomenon is called the kite lock-in.

4. Discussion

A method coupling a time domain seakeeping and a zero-mass kite modeling was developed to investigate the dynamic motion of ship towed by kite. The zero-mass kite modeling neglects the inertial forces, its deformation and the dynamic behavior of the tether. The ship modeling was based on the linear seakeeping STF strip theory [9]. A monolithic approach coupling the kite and the ship modeling was compared to a segregated approach, which neglected the coupling terms. Both methods allowed fast computation. Implemented in Python, the monolithic approach computed faster than the real-time. The computation of the considered cases were five time faster than the real time on a processor Intel(R) Xeon(R) CPU E3-1220 v3 with a CPU frequency of 3.10 GHz.

The ship modeling provides satisfactory results in terms of heave and pitch motions. Since the roll motion is highly nonlinear, the linear prediction of the roll motion is less accurate. However, a similar evolution of the predicted roll motion was noticed between the proposed roll motion modeling and the modeling by Ikeda et al. [25] used as reference here. Consequently, a good confidence in the kite influence on the roll motion could be attributed in terms of evolution. Nevertheless, the value of the results in terms of roll motion should be considered with caution. Moreover, the ship motions were restricted to heave, pitch and roll motions and existing couplings, such as roll and sway coupling, were neglected.

A linear dependency of the kite aerodynamic characteristics was proposed to correct the kite velocity amplitude and the tether tension amplitude in order to get a better agreement between the experiments and the proposed kite modeling. The identification of the kite modeling was carried out with experimental data on a single run of 5 minutes. A more extensive validation should investigate the influence of the true wind speed and the influence of the trajectory position and orientation. Moreover, the determination of the kite aerodynamic coefficients was entirely based on experimental data. It could be interesting to be able to estimate by simulation the kite aerodynamic characteristics for a deeper understanding.

In case of a dynamic flight, the high tensions developed lead to small tether sags, justifying the assumption that the lines are straight and inelastic. However, the tether acts as an interface between the kite and the ship. The dynamic tether effect may have a significant influence on the interaction between the kite and the ship. In particular, additional theoretical research should be carried out to study the effect of the tether on the kite lock-in phenomenon.

Moreover, in order to validate the modeling approach proposed in this paper, several strategies could be implemented. Full scale experiments on a 6-meter long dedicated experimental boat towed by kite were carried out by Behrel et al. [33,43]. However, measurement of the environment is challenging at full scale. With towing tank tests, similarity issues are also difficult to solve, as Martin et al. [44] emphasize this in the case of offshore floating wind turbines. In order to get around the problem of similitude and to take benefit from the towing tank test conditions, the kite could be modeled with a hardware in the loop method such as proposed by Giberti and Ferrari [45] for classical sailing yacht.

5. Conclusions and Future Work

A dynamic modeling of a ship towed by a kite was developed by combining two modelings. One for the ship and the other for the kite. Each modeling was validated independently. Two coupling approaches between the kite and the ship were compared: a segregated approach where the kite force was applied as a predefined time series and a monolithic approach where all the coupling terms were taken into account.

The computation of the linear convolution term for the ship motion modeling was performed with a state-space approach in order to provide time efficient computations. The identification of the state-space modeling was revisited in order to fit the forward speed case requirements. Kite modeling was modified using linear correction terms for aerodynamic characteristics. These modifications improved significantly the predicted amplitude of kite velocity and tether tension.

In the calm water case, the coupling between the ship and the kite decreases the kite forces and consequently, the segregated approach was conservative compared to the monolithic one. Different tether lengths were tested with a constant angular trajectory size leading to an evolution of the kite excitation frequencies. For the considered case, the effect of the wind gradient led to a continuous increase in roll amplitude according to the tether length and excitation frequency. The vicinity between the kite harmonics and the ship natural roll frequency led to further local significant increases.

In case of regular beam waves, the interaction between the kite and the ship was more significant. The definition of the kite trajectory into the relative wind frame taking into account the vertical ship velocity led to unrealistic trajectories with very short radius of curvature. This problem was avoided by neglecting the vertical component of the vessel motion when defining the kite’s trajectory. Moreover, a kite lock-in phenomenon to the frequency of the waves was highlighted, which can be win-win. Indeed, for some configurations, the mean kite towing force was increased while the ship roll amplitude remained smaller than the case without kite.