A Modified Free Wake Vortex Ring Method for Horizontal-Axis Wind Turbines

Abstract

A modified free-wake vortex ring model is proposed to compute the dynamics of a floating horizontal-axis wind turbine, which is divided into two parts. The near wake model uses a blade bound vortex model and trailed vortex model, which is developed based on vortex filament method with straight lifting lines assumption. By contrast, the far wake model is based on the vortex ring method. The proposed model is a good compromise between accuracy and computational cost, for example when compared with more complex vortex methods. The present model is used to assess the influence of floating platform motions on the performance of a horizontal-axis wind turbine rotor. The results are validated on the 5 MW NREL rotor and compared with other aerodynamic models for the same rotor subjected to different platform motions. The results show that the proposed method is reliable. In addition, the proposed method is less time consuming and has similar accuracy when comparing with more advanced vortex based methods.

1. Introduction

The unsteady aerodynamic loads experienced by floating offshore wind turbines (FOWTs) can be significantly different from those of bottom-fixed wind turbines [1]. Simulation tools that can accurately evaluate the relationship between the aerodynamic loads on the rotor and the platform motion are required especially for off-design conditions or novel concepts. Thus far, the available aero-hydro coupled analysis tools are almost exclusively based on the blade element momentum (BEM) theory, which describes the steady state behavior of a wind turbine but is unsuitable when the rotor interacts with its own wake, as is the case for FOWTs. An alternative is to use computational fluid dynamics (CFD) models that describe the flow physics around the rotor with higher fidelity. However, they are computationally expensive when dealing with the fully-coupled dynamics of FOWTs. Another approach, which is more accurate than momentum theory and less costly than CFD, is the so-called vortex methods. Vorticity-based methods have a number of different formulations, ranging from simple analytical models to more advanced numerical methods [2]. In addition, the computational cost of vortex models depends on the number of vortex elements and the interaction between them, ranging from relatively computationally less expensive models to more expensive ones [3,4]. Since the rotating wind turbines generate bounded circulations on the blades and release vorticity into the wake, the vortex based methods are naturally suitable for the simulation of wind turbine aerodynamics. The blade vorticity is assumed to be concentrated on lifting lines, with distributed vortex strengths representing the bound circulation, which convect to the wake as shed and trailed vorticity. The shed vorticity is the vorticity emitted in the wake due to time change of the bound circulation, while the trailed vorticity results from the spanwise change of bound circulation [2]. Within the scope of simplified vorticity models for rotors, the wake flow can be represented by vortex filaments such as the free wake vortex filament (FWVF) method of [5], vortex particles such as the nonlinear vortex lattice method (NVLM) method of [6], or a vortex ring model. As for the vortex ring methods, the wake is modeled with vortex rings with constant circulations determined by the aerodynamic loads on the rotor. Each vortex ring induces velocity on both the rotor and the other rings in the wake. The radius and spacing of the rings are time-varying, which is the so-called free wake vortex ring method. In this method, the analytical solution of the induced velocity of an axisymmetric vortex ring derived from Biot–Savart law [7,8] is used to calculate the induced velocities in the wake as well as on the rotor, which reduces computation cost while retaining the basic physical properties in the wake.

The simplest vortex ring methods available in the literature represent the rotor with an actuator disc, while the wake is modeled as vortex rings released at the tip of the rotor [9,10,11]. These models were used to calculate unsteady aerodynamic loads on both fixed and moving rotors. However, some aerodynamic features such as angle of attack, and lift and drag forces on the rotor could not be accurately reproduced due to the limitations associated with the actuator disc concept. A modeling method that represents the rotor blades with lifting lines, the near wake as a series of straight vortex lines and the far wake by vortex rings was initially suggested by Miller [12], relating the blade bound vortex strength with the vortex ring strength in the wake. In that case, the rotor loads were also dependent on the induced velocities from the vortex rings. Afjeh [13] and de Vaal [14,15] adapted and improved this method to the modeling of wind turbines, and good agreement was achieved when comparing the numerical results with measurement data.

The proposed free wake vortex ring method (FWVR) is based on the work of de Vaal [15] which is an aerodynamic model of a horizontal-axis wind turbine. It combines the Kutta–Joukowski lift theorem with the blade element theory, which can accurately and effectively predict the blade load using the time-evolution of the induced velocities. It splits the wake into near wake and far wake, leading to more realistic velocities on the blade and in the far-field. However, the application of de Vaal’s original method [15] can be challenging for some working states, for example the case of small or negative angles of attack at low wind speeds, or blade pitching motions at high wind speeds. In this paper, we present a few modifications of the original method to overcome these challenges. First, the propagation process of the vortex rings is modified and different nonlinear iteration methods are tested to solve the blade bound vortex strengths. Second, the near-wake trailed vortex model uses finite length of vortex segments instead of semi-infinite vortex lines, which is physically more realistic.

The proposed free wake vortex model was validated by analyzing the aerodynamic parameters on the NREL 5 MW wind turbine with both fixed-bottom and different types of platforms foundations, i.e., ITI barge, OC3-Hywind spar-buoy (OHS), and tension leg platform (TLP). The results were compared with those obtained from BEM theory, RANS method of [16] and two vortex based aerodynamic models—the FWVF method of [5] and the NVLM method of [6].

The paper is organized as follows. The theory of the modified free wake vortex ring method is described in Section 2. The numerical modeling method is described and validated in Section 3. The computational time evaluation is shown in Section 4. The wind turbine aerodynamic analysis is shown in Section 5. Finally, conclusions are drawn in Section 6.

2. Vortex Ring Theory

The proposed free wake vortex ring method consists of two parts: the near wake model and the far wake model. The induced velocities of the near wake model are calculated based on Biot–Savart law, as described in Section 2.1. To keep a low computational cost, the far wake model is represented by the axisymmetric vortex rings, as introduced in Section 2.2.

2.1. Velocity Induced by a Vortex Filament

The velocity ${\mathbf{V}}_P$ induced by a vortex filament with a constant strength $\Gamma$ in the field at a point P can be expressed by Biot–Savart law as [17]

$${\mathbf{V}}_P=-\frac{\Gamma }{4\pi }{\int }_{C(q)}\frac{{\mathbf{r}}_P-{\mathbf{r}}_q}{|{\mathbf{r}}_P-{\mathbf{r}}_q{|}^3}\times \frac{\partial {\mathbf{r}}_q}{\partial q}dq,$$

(1)

where $C(q)$ is the parametric curve which describes the path of the vortex filament, ${\mathbf{r}}_P$ is the position vector of the point P, ${\mathbf{r}}_q$ is the position vector of a point Q on the filament, $\mathbf{r}={\mathbf{r}}_P-{\mathbf{r}}_q$ is the vector pointing from point Q to point P, and $\frac{\partial {\mathbf{r}}_q}{\partial q}$ is the partial derivative of ${\mathbf{r}}_q$ with respect to the filament parameter. In Equation (1), the kernel of Biot–Savart operator is identified as

$$\mathbf{K}(\mathbf{r})=-\frac{1}{4\pi }\frac{\mathbf{r}}{{\left|\mathbf{r}\right|}^3},$$

(2)

which singular when the point P is on the vortex filament $C(q)$. To avoid the singularity, a smoothing method proposed by Hoydonck et al. [18] is used, which replaces $\mathbf{K}(\mathbf{r})$ with ${\mathbf{K}}_{\sigma }(\mathbf{r})$, defined as

$${\mathbf{K}}_{\sigma }(\mathbf{r})=-\frac{g(\sigma )}{{\left|\mathbf{r}\right|}^3}\mathbf{r},$$

(3)

where $\sigma$ is the non-dimensionalized length between the point P and a point on the filament, and $g(\sigma )$ is a three-dimensional velocity smoothing function [18]. Here, the Rosenhead–Moore velocity smoothing function is used, i.e.,

$$g(\sigma )=\frac{{\sigma }^3}{4\pi {({\sigma }^2+1)}^{\frac{3}{2}}}.$$

(4)

This leads to the following induction velocity,

$${\mathbf{V}}_P=\Gamma {\int }_{C(q)}{\mathbf{K}}_{\sigma }(\mathbf{r})\times \frac{\partial {\mathbf{r}}_q}{\partial q}dq.$$

(5)

Considering a straight line vortex filament, the parametric curve of the filament is given by $\mathbf{C}(q)={\mathbf{x}}_1+q({\mathbf{x}}_2-{\mathbf{x}}_1)$, with $0\le q\le 1$. Thus, the desingularized kernel ${\mathbf{K}}_{\sigma }(\mathbf{r})$ can be expressed as

$$\begin{array}{cc}\hfill {\mathbf{K}}_{\sigma }(\mathbf{r})& =-\frac{{\mathbf{x}}_P-\mathbf{C}(q)}{4\pi (|{\mathbf{x}}_P-\mathbf{C}(q){{|}^2+r_c^2)}^{\frac{3}{2}}}=-\frac{{\mathbf{x}}_P-(1-q){\mathbf{x}}_1+q{\mathbf{x}}_2}{4\pi (|{\mathbf{x}}_P-(1-q){\mathbf{x}}_1+q{\mathbf{x}}_2{{|}^2+r_c^2)}^{\frac{3}{2}}},\hfill \end{array}$$

(6)

where $r_c$ is the vortex core radius. Accordingly, the induced velocity ${\mathbf{V}}_P$ can be expressed as

$${\mathbf{V}}_P=-\Gamma {\int }_0^1\frac{{\mathbf{x}}_P-(1-q){\mathbf{x}}_1+q{\mathbf{x}}_2}{4\pi (|{\mathbf{x}}_P-(1-q){\mathbf{x}}_1+q{\mathbf{x}}_2{{|}^2+r_c^2)}^{\frac{3}{2}}}\times ({\mathbf{x}}_2-{\mathbf{x}}_1)dq.$$

(7)

When the point $x_2$ extends to infinity, the straight line vortex filament becomes a semi-infinite line. In this case, the upper limit of $q=\tilde{q}$ extends to infinity, and the induced velocity ${\mathbf{V}}_P$ becomes

$${\mathbf{V}}_P=-\Gamma \underset{\tilde{q}\to \infty }{lim}f(q),$$

(8)

$$f(q)={\int }_0^{\tilde{q}}\frac{{\mathbf{x}}_P-(1-q){\mathbf{x}}_1+q{\mathbf{x}}_2}{4\pi (|{\mathbf{x}}_P-(1-q){\mathbf{x}}_1+q{\mathbf{x}}_2{{|}^2+r_c^2)}^{\frac{3}{2}}}\times ({\mathbf{x}}_2-{\mathbf{x}}_1)dq.$$

(9)

2.2. Velocity Induced by an Axis-Symmetric Vortex Ring

Considering an ideal vortex ring with radius R, as shown in Figure 1, its parametric curve $C(\chi )$ with the angle $\chi$ as the parameterization parameter can be expressed as

$$C(\chi )=R(\cos \chi ,\sin \chi ,0),\chi \in [0,2\pi ].$$

(10)

In the local vortex ring coordinate frame $x^{{}^{\prime }}{y}^{{}^{\prime }}{z}^{{}^{\prime }}$ of an axis-symmetric vortex ring, the coordinates of P are $(x_P^{{}^{\prime }},0,z_P^{{}^{\prime }})$ and $C(\chi )$ is located in the $x^{{}^{\prime }}$$y^{{}^{\prime }}$-plane, which centered at the origin of the $z^{{}^{\prime }}$-axis. Thus, Taking the notation $x_P^{{}^{\prime }}=\eta R$ and $z_P^{{}^{\prime }}=\zeta R$, we have

$$\mathbf{r}(\chi )=R(\cos \chi ,\sin \chi ,0),$$

(11)

$${\mathbf{r}}_P=R(\eta ,0,\zeta ),\eta \in (0,+\infty ),\zeta \in (-\infty ,+\infty ).$$

(12)

The displacement vector $\mathbf{r}$ is thus given by

$$\mathbf{r}={\mathbf{r}}_P-\mathbf{r}(\chi )=R(\eta -\cos \chi ,-\sin \chi ,\zeta ).$$

(13)

The vortex ring induced velocity at P is obtained by substituting the above equations into Equation (5), which yields:

$$\mathbf{V}=-\frac{\Gamma }{4\pi R}{\int }_0^{2\pi }\frac{(\zeta \cos \chi ,\zeta \sin \chi ,1-\eta \cos \chi )}{{(1+{\eta }^2+{\zeta }^2+{\sigma }_c^2-2\eta \cos \chi )}^{\frac{3}{2}}}d\chi ,$$

(14)

where ${\sigma }_c=r_c/R$ is the non-dimensional vortex core radius; a constant value of 0.0116 is taken where the self induced velocity agree with ring-center velocity [15]. An analytical result of this integral has been derived as [7,8]

$$V_x=\frac{\Gamma }{2\pi R{C}_0}\left[-K(m)+\frac{1+{\eta }^2+{\zeta }^2}{C_1^2}E(m)\right]\frac{\zeta }{\eta },$$

(15)

$$V_z=\frac{\Gamma }{2\pi R{C}_0}\left[K(m)+\frac{1-{\eta }^2-{\zeta }^2}{C_1^2}E(m)\right],$$

(16)

with

$$C_0^2=1+2\eta +{\eta }^2+{\zeta }^2+{\sigma }_c^2,$$

(17)

$$C_1^2=1-2\eta +{\eta }^2+{\zeta }^2+{\sigma }_c^2,$$

(18)

and the functions $K(m)$ and $E(m)$ being the first and second types of complete elliptic integrals, respectively, with m defined as:

$$m=\frac{4\eta }{C_0^2}.$$

(19)

3. Numerical Model

As explained in Section 2, the present free wake vortex ring model consists of two parts: the near wake model and the far wake model. The near wake model includes the blade bound vortex model and the trailed vortex model, which are both represented by vortex line segments. The velocities induced by these vortex lines are calculated based on Biot–Savart law, as described in Section 2.1. By contrast, the velocity induced by the far wake model is given by Equations (15) and (16). Section 3.1 gives more insights into the discretization of blades using vortex elements, while Section 3.2 discusses how the far wake model is implemented and how vortex rings propagate in the wake.

3.1. Near Wake Models

3.1.1. Blade Bound Vortex Model

Figure 2 illustrates the blade bound vortex model for a rotor of $N_b=3$ blades with straight lifting lines assumption. Each blade is discretized with a series of vortex line segments. The radial endpoints on the vortex segments are denoted by $r_j^t$, while the control points are denoted by $r_i^n$, where $r_i^n=(r_j^t+r_{j+1}^t)/2$ ($i=1,\dots ,N$ and $j=1,\dots ,N+1$), N being the number of vortex segments. A local blade coordinate system $x^b{y}^b{z}^b$ is defined as right-handed and centered at the rotor root. The $x^b$-axis is along the pitch axis and pointing towards the tip of the blade; the $y^b$-axis points to the trailing edge of blade and is parallel to the chord line at the zero-twist blade section; and the $z^b$-axis is orthogonal to both $x^b$ and $y^b$. In the local blade coordinate system of a blade b, the rotor lies in the $x^b{y}^b$-plane, as shown in Figure 2. The coordinates of a control point $n_i$ and an endpoint $t_j$ are then expressed as

$${\mathbf{n}}_i=\left\{\begin{array}{c}0\\ r_i^n\\ 0\end{array}\right\},{\mathbf{t}}_j=\left[\begin{array}{ccc}\cos \Delta {\theta }_b& -\sin \Delta {\theta }_b& 0\\ \sin \Delta {\theta }_b& \cos \Delta {\theta }_b& 0\\ 0& 0& 1\end{array}\right]\left\{\begin{array}{c}0\\ r_j^t\\ 0\end{array}\right\},$$

(20)

where $\Delta {\theta }_b=(bi-b)2\pi /N_b$ is the angle between blade $bi$ and blade b. Assuming ${\mathbf{x}}_P={\mathbf{n}}_i$, ${\mathbf{x}}_1={\mathbf{t}}_1$ and ${\mathbf{x}}_2={\mathbf{t}}_2$, and substituting Equation (20) into Equation (7), it can be found that the velocity induced by the segment $x_1{x}_2$ at the control point ${\mathbf{n}}_i$ equals zero in both the $x_b$-direction and the $y_b$-direction, and in the $z_b$-direction is equal to

$$V_{ij}^{zb}={\Gamma }_j^b\frac{r_i^n\sin \Delta {\theta }_b\left[r_2(r_j^t-r_i^n\cos \Delta {\theta }_b)-r_1(r_{j+1}^t-r_i^n\cos \Delta {\theta }_b)\right]}{4\pi r_1{r}_2\left[r_c^2+{r_i^n}^2-{(r_i^n\cos \Delta {\theta }_b)}^2\right]}=A_{ij}^b{\Gamma }_j^b,$$

(21)

where $A_{ij}^b$ is the influence coefficient which only depends on the geometry and discretization of the rotor, and $r_1$ and $r_2$ are the lengths of the vectors ${\mathbf{x}}_1-{\mathbf{x}}_P$ and ${\mathbf{x}}_2-{\mathbf{x}}_P$, respectively, i.e.,

$$r_1=\sqrt{r_c^2+{r_i^n}^2+{r_j^t}^2-2r_i^n\cos \Delta {\theta }_b},$$

(22)

$$r_2=\sqrt{r_c^2+{r_i^n}^2+{r_{j+1}^t}^2-2r_i^n\cos \Delta {\theta }_b}.$$

(23)

The total induced velocity at a control point $n_i$ is the sum of the influences of all the bound vortex segments on all $N_b$ blades, which can be expressed as

$$V_i^{zb}=\sum _{b=1}^{N_b}\sum _{j=1}^N{V}_{ij}^{zb}=\sum _{b=1}^{N_b}\sum _{j=1}^N{A}_{ij}^b{\Gamma }_j^b.$$

(24)

Considering that the rotor is axis-symmetric and that all the blades’ geometry and discretization are identical, the total induced velocity at all N control points on a blade b can be simplified as

$$\left\{V^{zb}\right\}=\sum _{b=1}^{N_b}{\mathbf{A}}^b\left\{{\Gamma }^b\right\},$$

(25)

where ${\mathbf{A}}^b$ is the matrix of influence coefficients of the blade bound vortex model.

3.1.2. Trailed Vortex Model

In this section, a finite length vortex line model is introduced as opposed to the semi-infinite vortex line model proposed by de Vaal [14], as shown in Figure 3a. A local blade coordinate system $x^b{y}^b{z}^b$ is used, where the rotor lies in the $x^b{y}^b$-plane and the blade b lies along the $x^b$-axis. In the present model, the trailed vortex line begins from an endpoint $r_j^t$ on a blade and extends to a finite length $l_j^t$ as defined in Equation (26) in the direction normal to the blade in the $x^b{y}^b$-plane.

$$l_j^t=r_j^t{\theta }_t,$$

(26)

which equals to the arc length of $r_j^t$ times the angle ${\theta }_t$ swept, where ${\theta }_t$ is an input parameter.

The induced velocity at the control point $n_i$ on blade b can be determined using Equation (7). The coordinates of ${\mathbf{x}}_P$ and ${\mathbf{x}}_1$ are both given by Equation (20), while the coordinates of ${\mathbf{x}}_2$ can be written as

$${\mathbf{x}}_2=\left[\begin{array}{ccc}\cos \Delta {\theta }_b& -\sin \Delta {\theta }_b& 0\\ \sin \Delta {\theta }_b& \cos \Delta {\theta }_b& 0\\ 0& 0& 1\end{array}\right]\left\{\begin{array}{c}{r}_j^t{\theta }_1\\ r_j^t\\ 0\end{array}\right\}.$$

(27)

Substituting Equations (20) and (27) into Equation (7), evaluating the integral and simplifying, the velocity induced by $x_1{x}_2$ at ${\mathbf{n}}_i$ is derived as

$$V_{ij}^{zt}=\Delta {\Gamma }_j^t\frac{(r_j^t-r_i^n\cos \Delta {\theta }_b)\left[r_2^{{}^{\prime }}{r}_i^n\sin \Delta {\theta }_b-r_1(r_j^t-r_i^n\sin \Delta {\theta }_b)\right]}{4\pi r_1{r}_2^{{}^{\prime }}\left[r_c^2+{r_j^t}^2-2r_i^n{r}_j^t\cos \Delta {\theta }_b-{(r_i^n\cos \Delta {\theta }_b)}^2\right]}=A_{ij}^t\Delta {\Gamma }_j^t,$$

(28)

where $r_1$ and $r_2^{{}^{\prime }}$ are the lengths of the vectors ${\mathbf{x}}_1-{\mathbf{x}}_P$ and ${\mathbf{x}}_2-{\mathbf{x}}_P$, respectively, where $r_1$ is given in Equation (22) and $r_2^{{}^{\prime }}$ can be defined as

$$r_2^{{}^{\prime }}=\sqrt{r_c^2+{r_i^n}^2+{r_j^t}^2(1+{\theta }_1^2)-2r_i^n{r}_j^t(\cos \Delta {\theta }_b+{\theta }_1\sin \Delta {\theta }_b)}.$$

(29)

The total trailed model induced velocity at a control point $n_i$ is the sum of the induce of all the $N+1$ trailing vortex segments on all $N_b$ blades, which can be expressed as

$$V_i^{zt}=\sum _{b=1}^{N_b}\sum _{j=1}^{N+1}{V}_{ij}^{zt}=\sum _{b=1}^{N_b}\sum _{j=1}^{N+1}{A}_{ij}^t\Delta {\Gamma }_j^t,$$

(30)

which can be rewritten as

$$\left\{V^{zt}\right\}=\sum _{b=1}^{N_b}{\mathbf{A}}^t\left\{\Delta {\Gamma }^t\right\},$$

(31)

where ${\mathbf{A}}^t$ is the matrix of influence coefficients of the trailed vortex model.

3.1.3. Total Induced Velocity of the Near Wake

The total near wake induced velocity $V^{nw}$ is the sum up of the components of the blade bound vortex model in Equation (25) and the trailed vortex model in Equation (31). With non-zero values only in the $z_b$-direction, $V^{nw}$ is given in Equation (32), which only depends on the rotor geometry.

$$V^{nw}=\sum _{b=1}^{N_b}\left({\mathbf{A}}^b\left\{{\Gamma }^b\right\}+{\mathbf{A}}^t\left\{\Delta {\Gamma }^t\right\}\right).$$

(32)

3.2. Far Wake Model

The far wake model is formulated based on the assumption that the spiral vortices in the wake can be represented by a series of vortex rings, as illustrated in Figure 4. When all the blades together sweep a full revolution, a pair of new vortex rings $R^1$ is shed into the wake, the outer ring being shed from the outboard of the rotor and the inner ring being shed from the inboard of the rotor disc. The computation of the far wake induced velocity ${\mathbf{V}}_P^{fw}$ at a point P with $N_R$ pairs of vortex rings is given in Equation (33), with the theory as described in Section 2.2.

$${\mathbf{V}}_P^{fw}=\sum _{k=1}^{N_R}{\mathbf{V}}_{in,k}+\sum _{k=1}^{N_R}{\mathbf{V}}_{out,k},$$

(33)

where ${\mathbf{V}}_{in,k}$ is the induced by the kth inner ring and ${\mathbf{V}}_{out,k}$ is induced by the outer ring.

The Propagation of Vortex Ring

The far wake model is time-dependent, where two time scales are considered: (i) the time step for the vortex shedding is $\Delta T=T_P/N_b$, where $T_P$ is the period of the rotor rotation; and (ii) the time step used to compute the vortex rings propagation is defined as $\Delta t=\Delta T/n_{cst}$, where $n_{cst}$ is a constant integer. The calculation of the propagation is a cyclic process which involves three steps per $\Delta t$ at a time t:

• Identify the control points on the vortex rings and calculate the velocities (induced velocity and free stream velocity) on all the control points.
• Calculate the position of the control points both on the rotor and in the wake. The control points in the wake is given by Euler’s equation assuming an incompressible and inviscid fluid, as

  $$\mathbf{S}(t+\Delta t)=\mathbf{S}(t)+\mathbf{V}\Delta t,$$

  (34)

  where V is the speed of a control point in the global coordinate system.
• Update the position of the vortex rings based on the control points determined in step 2.

The geometry parameters that describe a vortex ring k are collected into the variable ${\mathbf{S}}_k$, as

$${\mathbf{S}}_k=\left\{\begin{array}{c}{x}_{O,k}\\ y_{O,k}\\ z_{O,k}\\ R_k\\ {\chi }_k\\ {\beta }_k\\ {\gamma }_k\end{array}\right\}k=1,\dots ,N_R,$$

(35)

where ${\mathbf{O}}_k$ is the vortex ring center, $R_k$ is the ring radius, and ${\gamma }_k$, ${\chi }_k$ and ${\beta }_k$ are, respectively, roll, pitch and yaw angles, all of which being functions of time. Here, we consider a local inertial coordinate system $\widehat{x}\widehat{y}\widehat{z}$ by translating the origin of the global coordinate coincide with that of a vortex ring k. At the upright position, when the ring lies in the $\widehat{x}\widehat{y}$-plane, ${\chi }_k$ and ${\beta }_k$ are both zero, and when it rotates anti-clockwise around the corresponding axis according to the right-hand law, the value of the angle is defined as positive. Due to the rotational invariance property of the ring, the roll angle ${\gamma }_k$ around the $\widehat{z}$-axis has no influence on the position of the ring itself but represents the azimuthal angle between the ring and the control point. In the following, the method to determine these parameters is introduced.

Let us define a local vortex ring coordinate system $x^{{}^{\prime }}{y}^{{}^{\prime }}{z}^{{}^{\prime }}$ by giving the local inertial coordinate system $\widehat{x}\widehat{y}\widehat{z}$ a pitch angle of ${\chi }_k$ and then a yaw angle of ${\beta }_k$, so that the vortex ring lies in the $x^{{}^{\prime }}{y}^{{}^{\prime }}$-plane. A vortex ring in the wake is represented by a series of control points distributed on the ring with the azimuthal angle $\Delta \theta$. With $N_c$ control points on a vortex ring, $\Delta \theta =2\pi /N_c$. Thus, the coordinate ${\mathbf{X}}_{i,k}^{{}^{\prime }}=(x_{i,k}^{{}^{\prime }},y_{i,k}^{{}^{\prime }},z_{i,k}^{{}^{\prime }})$ of a control point i on a ring k in the local vortex ring coordinate system is given as

$$\left\{\begin{array}{cc}& x_{i,k}^{{}^{\prime }}=R_k\sin [\Delta \theta (i-1)]\\ & y_{i,k}^{{}^{\prime }}=R_k\cos [\Delta \theta (i-1)]\\ & z_{i,k}^{{}^{\prime }}=0\end{array}\right.i=1,\dots ,N_c,k=1,\dots ,N_R.$$

(36)

which need to be expressed in the global coordinate frame denoted by ${\mathbf{X}}_{i,k}(t-\Delta t)$. At a time step t, the position of the control point denoted by ${\mathbf{X}}_{i,k}(t)$ is determined by Equation (34) as:

$${\mathbf{X}}_{i,k}(t)={\mathbf{X}}_{i,k}(t-\Delta t)+{\mathbf{V}}_{i,k}\Delta t.$$

(37)

Since the control points move independently, the new positions of the control points may no longer form an exact circle. This discrepancy is however very small if the time step $\Delta t$ is small enough. Nevertheless, a new vortex ring with the geometry parameters as shown in Equation (35) needs to be determined by the updated position of the azimuthal control points. Firstly, the vortex ring origin coordinate ${\mathbf{O}}_k$ can be expressed by taking the average of the coordinates of all the control points as

$${\mathbf{O}}_k=\frac{1}{N_c}\sum _{i=1}^{N_c}{\mathbf{X}}_{i,k}=\frac{1}{N_c}\sum _{i=1}^{N_c}(x_{i,k},y_{i,k},z_{i,k}).$$

(38)

Next, the vortex ring radius $R_k$ are determined by taking the average length of the segments ${\mathbf{s}}_{i,O}$ from each control point i to the center of the ring ${\mathbf{O}}_k$ as

$$R_k=\frac{1}{N_c}\sum _{i=1}^{N_c}\sqrt{{(x_{i,k}-x_{O,k})}^2+{(y_{i,k}-y_{O,k})}^2+{(z_{i,k}-z_{O,k})}^2}.$$

(39)

3.3. Characteristics of the First Rings Shed in the Wake

After a certain distance behind the rotor, the trailed vortex segments of the near wake roll up into two concentrated vortex rings in the far wake. The initial size and vortex strength of the vortex rings are defined in this section. The algorithms of the roll-up process is based on the momentum conservation theory [14,19]. The location of the maximum blade bound vortex strength ${\Gamma }_{max}$ is used to define the inner and outer rings. In particular, the near wake trailing vortices from ${\Gamma }_{max}$ to the tip of the blade are roll up into an outer ring and from ${\Gamma }_{max}$ to the root of the blade roll up into an inner ring. The vortex strength of the roll up vortex ring equals to the summation of the vortex strengths of the trailing vortices from which it is formed. The vortex ring strengths ${\overline{\Gamma }}_{in}^b$ and ${\overline{\Gamma }}_{out}^b$ as well as the vortex ring radii ${\overline{r}}_{in}^b$ and ${\overline{r}}_{out}^b$ on blade b can be expressed as

$${\overline{\Gamma }}_{in}^b=-{\int }_{r_1}^{r_2}\frac{d{\Gamma }^b}{dr}dr=\sum _{n=1}^{n_{max}}\Delta {\Gamma }_n^t,$$

(40)

$${\overline{r}}_{in}^b=\frac{1}{{\overline{\Gamma }}_{in}^b}{\int }_{r_1}^{r_2}r\frac{d{\Gamma }^b}{dr}dr=\frac{1}{{\overline{\Gamma }}_{in}^b}\sum _{n=1}^{n_{max}}{r}_n\Delta {\Gamma }_n^t,$$

(41)

$${\overline{\Gamma }}_{out}^b=-{\int }_{r_1}^{r_2}\frac{d{\Gamma }^b}{dr}dr=\sum _{n=n_{max}+1}^{N+1}\Delta {\Gamma }_n^t,$$

(42)

$${\overline{r}}_{out}^b=\frac{1}{{\overline{\Gamma }}_{out}^b}{\int }_{r_1}^{r_2}r\frac{d{\Gamma }^b}{dr}dr=\frac{1}{{\overline{\Gamma }}_{out}^b}\sum _{n=n_{max}+1}^{N+1}{r}_n\Delta {\Gamma }_n^t.$$

(43)

The concentrated vortex strength and radius of the inner and outer ring are defined as an average of the vortex strengths and radii calculated from each blade, that is

$${\overline{\Gamma }}_{out}=\frac{1}{N_b}\sum _{b=1}^{N_b}{\overline{r}}_{out}^b,{\overline{\Gamma }}_{in}=\frac{1}{N_b}\sum _{b=1}^{N_b}{\overline{r}}_{in}^b,$$

(44)

$${\overline{r}}_{out}=\frac{1}{N_b}\sum _{b=1}^{N_b}{\overline{r}}_{out}^b,{\overline{r}}_{in}=\frac{1}{N_b}\sum _{b=1}^{N_b}{\overline{r}}_{in}^b.$$

(45)

3.4. Strength of the Blade Bound Vortex

Based on the discussion above, it can be seen that both the trailing vortex strengths and the far wake vortex ring strengths are determined by the distribution of blade bound vortex strengths. Once the blade bound vortex strengths are known, the near wake induced velocity at any point of the rotor can be determined by Equation (32) and the far wake vortex ring strengths are determined at the roll-up processes and with no change during their convection in the wake. Thus, the blade bound vortex strengths need to be determined as a prerequisite. A method to calculate the blade bound vortex strengths based on the Kutta–Joukowski lift theorem and the blade element theory is discussed below.

3.4.1. Kutta–Joukowski Lift Theorem

The Kutta–Joukowski lift theorem as introduced by Phillips [20] is adopted to solve the blade bound vortex strength. Based on this theory, the lift force $d{\mathbf{F}}_l$ acting on a vortex segment $d\mathbf{l}$ can be calculated by the product of the air density $\rho$, the bound vortex strength ${\Gamma }^b$ on the segment, and the relative wind velocity $\mathbf{V}$, as

$$d{\mathbf{F}}_l=\rho {\Gamma }^b\mathbf{V}\times d\mathbf{l},$$

(46)

where $\mathbf{V}$ includes the contribution from the free stream ${\mathbf{V}}_{\infty }$, the rotor rotation ${\mathbf{V}}^{\Omega }$, the near wake induced velocity ${\mathbf{V}}^{nw}$, and the far wake induced velocity ${\mathbf{V}}^{fw}$, as well as the motion of the platform ${\mathbf{V}}_p$ in the case of a floating offshore wind turbine. $\mathbf{V}$ has three components in the blade coordinate system as $V_{xb}$ in the $x^b$-direction, which is along the pitch axis towards the tip of the blade, and $V_{yb}$ and $V_{zb}$ lie in the blade section and normal to the blade pitch axis, which can be written as

$$\mathbf{V}={\mathbf{V}}_{\infty }+{\mathbf{V}}^{\Omega }+{\mathbf{V}}^{nw}+{\mathbf{V}}^{fw}+{\mathbf{V}}_p=\left\{\begin{array}{c}{V}_{xb}\\ V_{yb}\\ V_{zb}\end{array}\right\},$$

(47)

Substituting Equation (47) into Equation (46), and given that the $x^b$-axis is parallel to $d\mathbf{l}$, it is found that only the $V_{yb}$ and $V_{zb}$ components contribute to the lifting force acting on the blade. Thus, Equation (46) can be rewritten as

$$d{\mathbf{F}}_l=\rho {\Gamma }^b{\mathbf{V}}^n\times d\mathbf{l},$$

(48)

where ${\mathbf{V}}^n$ is the velocity vector that has only the $V_{yb}$ and $V_{zb}$ components, i.e.,

$${\mathbf{V}}^n=\left\{\begin{array}{c}0\\ V_{yb}\\ V_{zb}\end{array}\right\}.$$

(49)

3.4.2. Blade Element Force

The lift force on a blade segment can also be calculated from the airfoil properties of the segment, as indicated in Figure 5. The lift force can be expressed as

$$d{\mathbf{F}}_l=\frac{1}{2}\rho {\left|{\mathbf{V}}^n\right|}^2{C}_{l}cdl,$$

(50)

where ${\mathbf{V}}_n$ is defined in Equation (49), c is the chord length, and $C_l$ is the lift coefficient obtained by means of tabulated aerodynamic parameters (the tabulated data are taken from the NREL certification cases of FAST: https://nwtc.nrel.gov/FAST8) of the airfoil with a given $\alpha$, which is calculated as

$$\alpha =\varphi -({\theta }_t+{\theta }_p),$$

(51)

where ${\theta }_t$ is the twist angle, ${\theta }_p$ is the collective pitch angle, and the angle of relative wind $\varphi$ is the angle between the wind velocity ${\mathbf{V}}_n$ and the plane of the rotor disc, which can be expressed as

$$\varphi =ta{n}^{-1}\frac{V_0(1-a)}{\Omega r(1+a^{\prime })}.$$

(52)

In Equation (52), $\Omega$ is the rotor angular speed, a is the axial induction factor and $a^{\prime }$ is the angular induction factor. Similar to the lift force, the drag force $d{\mathbf{F}}_d$ is calculated from the tabulated $C_d$ and $\alpha$, leading to

$$d{\mathbf{F}}_d=\frac{1}{2}\rho {\left|{\mathbf{V}}^n\right|}^2{C}_{d}cdl.$$

(53)

Based on the discussion above, it is found that both Equations (48) and (50) express the lift force of the blade segment. Thus, by equating them for a blade section i, the following relation is obtained

$${\Gamma }_i^b|{\mathbf{V}}_i^n({\Gamma }_j)\times d{\mathbf{l}}_i|-\frac{1}{2}{\left|{\mathbf{V}}_i^n({\Gamma }_j)\right|}^2{C}_l({\alpha }_i)c_{i}d{l}_i=0.$$

(54)

For each blade section $i=1,\dots ,N_b·N$, an equation similar to Equation (54) can be set up, where j is the number of the vortex ring. For all $N_b·N$ blade sections, a system of $N_b·N$ non-linear equations are obtained and solved for $N_b·N$ unknown values of the bound vortex strength ${\Gamma }_i^b$. Different methods can be used to solve the resulting set of non-linear equations. The trust-region methods are used here, which proved to be more stable and converge more rapidly than Newton–Raphson methods, especially for small and negative angles of attacks. Accordingly, the rotor thrust ${\mathbf{F}}_T$ is expressed as

$${\mathbf{F}}_T=\sum _{b=1}^{N_b}\sum _{j=1}^{N}d{\mathbf{F}}_N=\sum _{b=1}^{N_b}\sum _{j=1}^{N}d{\mathbf{F}}_l\cos \varphi +d{\mathbf{F}}_d\sin \varphi .$$

(55)

4. Simulation Description

This simplified free wake vortex model aims at finding a compromise between physical accuracy and low computational cost. In this section, we compare the computational cost of the present method with other methods available in the literature.

Table 1. Computational time and function of different models.

Method	Computational Time	Modeling Capability
BEM	${\displaystyle C_b{N}_{\theta }{N}_{\Omega }}$	Rotor aerodynamics without wake
AFWVR	${\displaystyle \frac{1}{3}{C}_a{N}_b^2{N}_{\theta }^2{N}_{\Omega }^3}$	Rotor aerodynamics, uncoupled wake aerodynamics
FWVR	${\displaystyle \frac{1}{3}{C}_r{N}_b^2{N}_{\theta }^2{N}_{\Omega }^3+C_{\Gamma }{N}_{\theta }{N}_{\Omega }}$	Coupled rotor aerodynamics and wake aerodynamics
FWVF	${\displaystyle \frac{1}{3}{C}_f{N}_b^2{N}_{\theta }^3{N}_{\Omega }^3+C_{\Gamma }{N}_{\theta }{N}_{\Omega }}$	Coupled rotor aerodynamics and wake aerodynamics
NVLM	${\displaystyle \frac{1}{3}{C}_m{N}_b^2{N}_{\theta }^3{N}_{\Omega }^3+C_{\Gamma }{N}_{\theta }{N}_{\Omega }}$	Coupled rotor aerodynamics and wake aerodynamics

shows the methods under consideration with computational cost and modeling capability: blade element momentum theory (BEM), actuator disc with free wake vortex ring model (AFWRV) of [10,11], the present free wake vortex ring model (FWVR), and two more advanced vortex methods—the free wake vortex filament (FWVF) method of [5] and the nonlinear vortex lattice method (NVLM) method of [6].

In Table 1, $N_b$ is the number of rotor blades, $N_{\theta }$ is the number of time steps per rotor rotation and $N_{\Omega }$ is the number of rotor rotation to be simulated. $N_{\theta }{N}_{\Omega }$ represents the total number of time steps. $C_b$ represents the time cost of evaluating the induced velocity at a point in the field with BEM, $C_a$, $C_r$, $C_f$ and $C_m$ represent the time cost of evaluating the induced velocity of a vortex element (vortex rings or vortex filaments) at a point in the field with AFWVR, FWVR, FWVF and NVLM, respectively. Whether the blades are modeled or not, all four methods are evaluated with the same number of control points on the rotor and with the same number of time steps.

For the BEM theory, the number of control points are constant and each point is considered independently. Thus, the computational time has a linear relationship with the number of time steps, which is the fastest method for aerodynamic simulations. $C_b$ involves the iterative solution of BEM equations and some analytical corrections, such as Prandtl tip-loss, Prandtl hub-loss, and Pitt and Peters skewed-wake corrections when needed.

By contrast, the solution of the other four vortex methods mainly includes two parts: the calculation of the vortex element strength and the calculation of the induced velocity including self-induced velocities and mutually induced velocities. The equations need to be solved increases with the increase of vortex elements (being vortex rings or vortex filaments) in the wake.

The vortex ring strength of AFWVR is directly determined from a prescribed thrust coefficient $C_t$ [10], where the computational time is negligible. By contrast, the FWVR and FWVF methods determine the blade bound vortex strength by iteratively solving the equations based on the 3D vortex lifting law and the blade element theory with tabulated data of lift coefficient $C_l$, and then the vortex strengths are distributed to the new generated vortex rings or vortex filaments in the wake. The computational times to determine the vortex strengths for FWVR, FWVF and NVLM methods at each time step are identical and denoted by $C_{\Gamma }$. Since $C_{\Gamma }$ and $C_b$ are both iterative processes with the same number of equations, they are the same order of magnitude.

The computational times for evaluating the induced velocities of FWVR and FWVF are given in [14,21]. Since the number of vortex filaments is $N_{\theta }$ times the number of vortex rings, the number of induced velocities that need to be calculated at each time step for FWVR is also $N_{\theta }$ times that of FWVF. Since AFWVR has the same number of control points and vortex rings as FWVR, the number of induced velocities that needs to be calculated is the same. $C_a$ involves the analytical solution of a vortex ring induced velocity, which includes the evaluation of the elliptic integrals of the first and second kind. $C_r$ involves the same calculation as $C_a$ as well as the analytical solutions of near wake induced velocities at the rotor points where the time cost is very small. Thus, $C_r$ is approximately equal to $C_a$. $C_f$ involves the solving of the trigonometry functions based on Biot–Savart law, and $C_m$ involves the solving of Biot–Savart kernel functions. Thus, $C_f$ and $C_m$ are considered to be the same order of magnitude and slightly less expensive than $C_r$ [14].

As expected, the above analysis shows that more physically accurate models are also more computationally expensive. BEM is the fastest model among those considered here. AFWVR and FWVR differ in solving the blade bound vortex strength equations, which is approximately the computational time of BEM, while FWVR, FWVF and NVLM differ in a factor $N_{\theta }$ for the computation of the induced velocities, which can be significant. Regarding the modeling capability, BEM evaluates the rotor aerodynamics without considering the wake aerodynamics. AFWVR evaluates both rotor and aerodynamics, but uses a prescribed thrust coefficient to determine the vortex strengths. This prevents further analysis on the influence of the wake on the rotor. FWVR, FWVF and NVLM evaluate the coupled rotor aerodynamics and the wake aerodynamics, thus presenting similar capabilities, despite the fact that FWVR discretizes the wake into vortex rings which partly loss the tangential induction. Thus, Table 1 shows that the present simplified numerical model has low computational cost and relatively high performance compared to other reference models.

5. Results

The results are shown for the NREL 5 MW reference turbine, with rotor blades and tower assumed to be rigid. The basic parameters of the turbine, the airfoils and aerodynamic properties, the rotor speed and blade pitch data corresponding to wind speeds from cut-in to cut-out and other information about NREL 5 MW reference turbine that are used in this model can be found in the NREL report [22]. In this section, the results are shown for three cases: bottom-mounted, floating restricted to one degree of freedom, and floating in three degrees of freedom.

5.1. Bottom-Mounted Wind Turbine

In this section, a bottom-mounted wind turbine with a constant free stream is used to validate the present model. The rotor thrusts $F_T$ for wind speeds varying from the cut-in value of 3 m/s to the cut-out value of 25 m/s are shown in Figure 6. It is compared to different results from the literature, including the results from NREL [22], the result calculated with FAST v8, the results from [16] where the Reynolds-averaged Navier–Stokes (RANS) method was used, and the results from [6] that uses a nonlinear vortex lattice method (NVLM).

Figure 6 shows good agreement of the thrusts with the literature. Form the cut-in wind speed to the rated wind speed, the thrusts from the NVLM model, RANS model and the present FWVR model are higher than the thrusts from FAST v8, while the thrusts at above-rated wind speeds from these three models are lower than the thrusts from FAST v8. However, it should be noted that the wind turbine used by the RANS model is not exactly NREL 5 MW reference turbine but a similar one, which can be a good reference. The thrusts of the NREL 5 MW reference turbine released by NREL [22] are found to be not pure aerodynamic thrusts, but also include the effects of gravity and inertial loads of the rotor (https://wind.nrel.gov/forum/wind/viewtopic.php?f=2&t=917), which are significantly higher than the pure aerodynamic thrusts.

The axial induced velocities of wind speeds 6 m/s, 11.4 m/s and 18 m/s of the bottom-mounted wind turbine are compared with the results from FAST v8 as shown in Figure 7. It can be seen that the two results from 20 m to 60 m along the blade match well with each other. Due to the tip- and root-corrections, the inductions from FAST v8 at the tip and root of the blade are both zero. The induction from the FWVR at the tip of the blade also show a drop for all the three wind speeds, which are mainly contributed by the near wake model and the outer ring of the far wake model. At the root of the blade, the inner ring of the far wake model mainly plays a role of weakening the inductions while the inductions from the near wake model are still strong. With the increase of the wind speed, the effect of the near wake model becomes stronger than the far wake model. Thus, it can be seen that at the wind speed 6 m/s the induction dropped at the root but for the other two wind speeds the inductions are even stronger. The axial induction factors weighted with the swept area of the three wind speeds are 0.300, 0.262 and 0.044 from FWVR and 0.305, 0.242 and 0.035 from FAST v8, respectively.

5.2. Floating Wind Turbine under Single-DoF Motion

Since the results from the present vortex ring method are shown to be reliable for a bottom-mounted wind turbine, the method is applied to a floating offshore wind turbine with a platform motion prescribed in one degree of freedom (DoF).

5.2.1. Thrusts Comparison with Lee

The thrusts are calculated with a wind speed of 8 m/s and a rotor speed of 9.16 rpm. The platform motions are defined as a sine function in Equation (56) with amplitude (A) and frequency (f) given in

Table 2. Amplitude and frequency of prescribed single-DoF motions [6].

Motion	Amplitude A [m or °]	Frequency f [Hz]
Translation (surge, sway, heave)	4	0.1
Rotation (roll, pitch, yaw)	4	0.05
surge	4,8,12	0.03
pitch	2,4,6	0.03

.

$$X(t)=A_1\sin (2\pi f_{1}t).$$

(56)

Figure 8 shows the results of thrusts with translational (surge, sway, and heave, in Figure 8a) and rotational (yaw, pitch, and roll, in Figure 8b) platform motions. The thrust for a fixed rotor is also provided for comparison. For a moving rotor, the thrust forces obtained from both methods have similar amplitude and frequency, which are both consistent with the thrusts shown in Figure 6. Figure 8 shows that the thrust under translational platform motions is most sensitive to surge, and the thrust under rotational platform motions is most sensitive to pitch, for the reason that these motions are in the wind direction, which can significantly influence the relative wind speed normal to the rotor disc. When the rotor has a leeward speed, the relative wind speeds are reduced on the blades and the angle of attacks are smaller, which leads to a lower thrust load than average. When the rotor has a windward speed, the situation is reverse and the thrust load is larger than average. It is also observed that the heave and sway motions of the same amplitude and frequency theoretically have equivalent influence to the rotor thrust. The heave and sway motions can change the relative tangential wind speed on the airfoils, which either increases or decreases the angle of attack according to the azimuthal angle of the blade, thus increasing or decreasing the thrust. However, this influence is relatively small as the change of tangential wind speed is much smaller than the tangential speed coming from the rotor rotation. Figure 8a shows that, with the present FWVR method, the thrust force under heave and sway motions almost overlap and slightly fluctuate around the thrust value of the bottom-mounted turbine, which is deemed to be reasonable.

The effect of the roll motion on the thrust is similar to that of the heave and sway motions, which mainly influence the relative tangential wind speed and the wake induced velocity on the rotor. The yaw motion can reduce the overall wind area of the rotor, hence reducing the thrust load compared to the bottom-mounted turbine, as obtained for both NVLM and FWVR.

The effect of the surge and pitch motions with different amplitudes are further analyzed and compared with the results of [6] in Figure 9. Figure 9a compares the thrust with surge motions with the amplitude of 4 m, 8 m, 12 m and the frequency of 0.03 Hz, while Figure 9b compares the thrust with pitch motions with the amplitude of 2°, 4°, 6° and the frequency of 0.03 Hz. Since the amplitudes of the platform motions are small and the frequency is relatively low, the thrusts change linearly with platform motions for both surge and pitch.

5.2.2. Angle of Attacks Comparison with Sebastian and de Vaal

Additional comparisons are made with the works of Sebastian [5] and de Vaal [14], for different floater types, under the following operating and environmental conditions:

• below-rated: $V_{\infty }$ = 6.0 m/s, $\lambda$ = 9.63, $H_s$ = 1.83 m, $T_p$ = 12.72 s;
• rated: $V_{\infty }$ = 11.4 m/s, $\lambda$ = 7.00, $H_s$ = 2.54 m, $T_p$ = 13.35 s; and
• above-rated: $V_{\infty }$ = 18.0 m/s, $\lambda$ = 4.43, $H_s$ = 4.09 m, $T_p$ = 15.33 s, ${\theta }_p$ = 15°

In these works, the platform motions in time domain are synthesised as sinusoidal series at the first two dominant frequencies asthe chord length

$$X(t)=X_0+A_1\sin (2\pi f_{1}t+{\varphi }_1)+A_2\sin (2\pi f_{2}t+{\varphi }_2),$$

(57)

where the mean value $X_0$, amplitudes $A_i$, frequencies $f_i$ and phase angles ${\varphi }_i$ are given in

Table 3. Parameters of platform motions [5] for the NREL 5 MW turbine with different floater types.

Oper.	Platf.	Mode	${\mathit{X}}_0$	${\mathit{A}}_1$	${\mathit{f}}_1$	${\mathit{\varphi }}_1$	${\mathit{A}}_2$	${\mathit{f}}_2$	${\mathit{\varphi }}_2$
			[m]/[°]	[m]/[°]	[Hz]	[rad]	[m]/[°]	[m]/[°]	[rad]
1	ITI	surge	13.602	0.725	0.007	−1.163	−0.442	0.078	2.609
1	ITI	heave	−0.130	0.318	0.078	1.303	0.254	0.108	2.702
1	ITI	pitch	0.591	1.475	0.078	−0.066	1.630	0.083	1.816
1	OHS	pitch	1.580	−0.084	0.066	1.930	−0.116	0.077	3.113
1	OHS	yaw	−0.021	0.091	0.108	1.983	−0.036	0.120	3.429
1	TLP	surge	1.206	0.436	0.016	−0.831	−0.222	0.077	3.018
2	ITI	pitch	1.722	−0.637	0.065	−0.381	1.677	0.077	1.835
3	ITI	pitch	0.939	1.518	0.066	2.132	2.979	0.078	6.863
3	OHS	pitch	3.324	11.961	0.029	0.420	0.000	0.000	0.000
3	OHS	yaw	−0.222	2.000	0.029	−0.359	3.185	0.058	3.385

.

The mean ${\mu }_{\alpha }$ and standard deviation ${\sigma }_{\alpha }$ of the angle of attack $\alpha$ on the outboard 2/3 of the blade are compared in

Table 4. Mean and standard deviation of $\alpha$ at the outboard 2/3 part of the blade.

Operation	Platform	Mode	WInDS [5]		FWVR1 [14]		FWVR2, tilt = 0°		FWVR2, tilt = 5°
			${\mathit{\mu }}_{\mathit{\alpha }}$ [°]	${\mathit{\sigma }}_{\mathit{\alpha }}$ [°]	${\mathit{\mu }}_{\mathit{\alpha }}$ [°]	${\mathit{\sigma }}_{\mathit{\alpha }}$ [°]	${\mathit{\mu }}_{\mathit{\alpha }}$ [°]	${\mathit{\sigma }}_{\mathit{\alpha }}$ [°]	${\mathit{\mu }}_{\mathit{\alpha }}$ [°]	${\mathit{\sigma }}_{\mathit{\alpha }}$ [°]
1	Monopile	–	3.95	0.23	3.86	0.48	3.82	0.15	3.71	3.23
1	ITI	surge	3.95	0.40	3.87	0.53	3.78	0.23	3.64	3.23
1	ITI	pitch	3.99	2.21	3.90	1.5	3.89	1.92	3.84	3.85
1	OHS	pitch	3.94	0.32	3.84	0.49	3.66	0.25	3.84	3.21
1	TLP	surge	3.95	0.27	3.86	0.49	3.64	0.23	3.70	3.23
2	Monopile	–	6.76	0.37	6.66	0.69	5.84	0.35	5.76	3.35
2	ITI	pitch	6.78	1.67	6.67	1.30	5.82	1.09	5.73	3.47
3	Monopile	–	−0.10	0.80	−0.31	2.24	−0.59	0.05	−0.61	2.91
3	ITI	pitch	−0.08	2.26	−0.28	2.88	−0.59	1.87	−0.61	3.39
3	OHS	pitch	−0.45	3.59	−0.52	3.09	−0.83	2.40	−0.85	3.57

. The current numerical results are presented in the column entitled “FWVR2, tilt = 0°”, without shaft tilt angle. It is shown that ${\mu }_{\alpha }$ and ${\sigma }_{\alpha }$ have the same trend than observed in the literature. In particular, ${\mu }_{\alpha }$ is rather independent of the platform motion, for a given wind speed. The bottom-mounted monopile presents the smallest value of ${\sigma }_{\alpha }$ compared with floating support structures at the same wind speed. This can be explained by the fluctuations of the induced velocities caused by turbine-wake interactions in the floating case. The value of ${\sigma }_{\alpha }$ increases with the increase of platform motion. In addition, the pitch motions of the platform tend to have more effect on the value of ${\sigma }_{\alpha }$ than the surge motions, because the pitch motions generate shear wind velocities on the rotor. The values of ${\mu }_{\alpha }$ are lower than the reference results throughout and the difference are within 1°, which is considered to be acceptable. The column “FWVR2, tilt = 5°” shows the $\alpha$ performance in the design condition of NREL 5 MW with a nonzero tilt angle. It can be seen that the values of ${\mu }_{\alpha }$ are slightly lower than for zero tilt angle, which is mainly because the total wind area is smaller when the rotor is tilted. Moreover, ${\sigma }_{\alpha }$ is significantly influenced by the tilt angle in all cases because the latter can generate fluctuating streams especially on the outboard of the blades.

5.2.3. Thrust Analysis in Time Domain and Frequency Domain

The rotor thrust of monopile in rated condition and five floating support structures in below-rated, rated and above-rated conditions are evaluated in time domain (Figure 10) and frequency domain (Figure 11). Figure 10 shows that the thrust fluctuate around their mean values, after an initial transient. The latter is due to the fact that a certain time is needed for the far-wake vortex rings to develop. Thus, initially, only the near wake mostly contributes to the values of induced velocities. After about 50 s, a steady-state is reached and statistics can be taken. Figure 11 highlights the dominant frequencies in the rotor thrust. For a bottom-mounted turbine on a monopile in rated conditions, the first dominant frequency is at about 0.6 Hz, which corresponds to the vortex shedding frequency. Higher frequencies are multiples of the first one. This is due to a periodic process whereby the thrust force reaches a minimum value when the induced velocity reaches a peak value, and vice versa. For the FOWTs, without exception, the frequencies of the platform motions dominate in the frequency of the thrust force. This is because the velocities induced by the platform motions are larger than the induced velocities from the vortex rings.

Finally, we evaluate the blade bounded vortex strength $\Gamma$ and its relationship with other parameters. As explained above, the value of $\Gamma$ is obtained by equating the lift forces calculated from the 3D vortex lifting law and that from the blade element theory (see Equation (54)). This leads to:

$$\Gamma =\frac{c}{2}{C}_l(\alpha )V^n,$$

(58)

where c is the chord length, $C_l$ is the lift coefficient and $V^n$ is the normal wind speed. The time evolutions of $\Gamma$, $V^n$, $C_l$ and $\alpha$ are shown in Figure 12, for the OHS above-rated condition. In each figure, the continuous line represent the platform pitching motion. The non-dimensional $\Gamma /(\pi R{V}_{\infty })$ (Figure 12a) fluctuates with the platform pitch motion. In addition, the high-frequency oscillation with the period of about 4.96 s is due to the rotor tilt angle. On the outboard of the blade, there are negative values of $\Gamma$, which is consistent with the negative values of $\alpha$ and $C_l$ observed in the same conditions (Figure 12c,d). The normal relative wind speed $V^n$ is less influenced by the platform motion and the rotor tilt, its main fluctuations coming from the constant rotation of the rotor.

5.3. Floating Wind Turbine under Multiple-DoF Motion

To conclude this section, three multiple-DOF operating conditions are conducted, combining the properties listed in Table 3, namely:

• Below-rated: The ITI Energy barge with platform surge, heave, and pitch.
• Below-rated: The OC3-Hywind spar-buoy with platform pitch and yaw.
• Above-rated: The OC3-Hywind spar-buoy with platform pitch and yaw.

5.3.1. The Multiple-DoF Motion Thrusts Comparison with Lee

The multiple-DoF motion thrusts of ITI and OHS with the below-rated operating conditions as well as the fixed-rotor thrusts are compared with the results of [6] in Figure 13. It should be noted that, the platform motions start at a zero azimuth angle in the FWVR code while they start approximately after one rotation of the rotor in the NVLM code. Therefore, to facilitate the comparison, the thrusts from the present FWVR method is shifted with +360° along the azimuth angle-axis. It is found that the fixed-rotor thrust calculated with the present FWVR method is higher than that calculated with the NVLM method, which is consistent with the previous observation. Besides, the thrust amplitudes and frequencies of the two methods for the ITI and the OHS match well.

5.3.2. Induced Velocity Comparison with Sebastian

The wake induction at the rotor as a function of the downstream distance of the vortex rings are computed as the percentage of the total induced velocity at the rotor, and compared with [5] in Figure 14 for the ITI and in Figure 15 for the OHS with the NREL 5 MW wind turbine. The blue dashed lines represent the contribution of the induced velocities from each section along the wake as obtained with the FWVF method, while the red lines are the results of the present FWVR method. The red circles on the line highlight the position of the vortex rings. The black lines and the green dotted lines represent the accumulated percentage of the induced velocities along the wake from the FWVF method and FWVR method, respectively. The induced velocity data output from the present FWVR method was captured at the operating time of 150 s, when a stable condition was reached by the wind turbine. Since the convection of the vortex rings downstream is a dynamic process, the induced velocity captured at different time steps can be slightly different.

In general, the induced velocities are comparable between the two methods. As shown in Figure 14, the accumulated percentage of the induced velocity at distances 1D and 2D behind the rotor differ only by 2%, and this difference decreases to maximum 1.0% as the distance behind the rotor increases. Furthermore, the induced velocities from each section along the wake show similar trend with both methods, the big fluctuations being mainly influenced by the platform motions, particularly from the pitch. For the OHS below-rated condition (Figure 15a), a similar discrepancy is observed between the models, with differences of 2.4% and 2.1% at distances of 1D and 2D behind the rotor and 0.6% further downstream. In addition, the influence of the pitch and yaw platform motions on the induced velocity is relatively small. Above rated power (Figure 15b), slightly larger discrepancies are observed, with differences of 9.8% and 4.8% at distances of 1D and 2D behind the rotor and less than 1% further downstream. Additionally, the present FWVR method predicts significantly larger contribution of the induced velocity from the vortex rings within the 1D distance downstream. This is mainly because the vortex rings convect faster with the increase in wind speed, and the density of rings is sparser at higher wind speed. Because the induced velocity is very sensitive to the distance between the control points and the vortex ring (see Equation (21)) the contribution from the front rings is much more significant than that of the rings further away from the rotor. The small differences in the results can be explained from the differences in wake induction of these two methods. The first difference is that the induced velocities of the wake from the FWVF method are output for each vortex segment in the wake, while the induced velocities from the present FWVR method are output for each vortex ring. As explained in Section 4, there are fewer vortex rings than vortex segments in the wake, thus the induced velocities of the wake from the FWVF method are relatively constant compared to those from the FWVR method. The second difference is that, due to the different modeling properties, the induced velocities from the FWVF method start from a position very close to the rotor, while the induced velocities from the present FWVR method start from the position of the first vortex ring, and the induced velocity at the rotor from the near wake model is not accounted here.

5.3.3. Wake structure discussion

The wake structures of the bottom-fixed wind turbine with the specified below-rated operating condition as well as the FOWTs with the three multiple-DoF platform motions at 150 s of the operating time are shown in Figure 16. The inner and outer rings are distinguished with red and blue colors, respectively, with the vectors indicating the directions where the control points are moving in the next time step. The wake structures of the below-rated operating conditions, as shown in Figure 16a–c, are found to be more compact than the above-rated operating conditions, as shown in Figure 16d, which is due to the different convecting speeds of the vortex rings in the wake in different free stream environments. The wake structure of the bottom-fixed wind turbine has homogeneous vortex ring radii since its blade bound vortex strengths are homogeneous. By contrast, the wake structures of the FOWTs have significant fluctuations because both the blade bound vortex strengths and the rotor position change with time, thus affecting the roll-up process of the wake model. Additionally, the vortex rings in the wake with uneven radii and vortex strengths interacting with each other can cause some fluctuations downstream in the wake, especially below-rated conditions for which the downstream vortex rings are relatively close to one another, which can be seen clearly in Figure 16b,c.

6. Conclusions

In this paper, a simplified free wake vortex ring method for horizontal-axis wind turbine is modified and assessed. A new trailed vortex model with finite length vortex lines is proposed. The NREL 5 MW wind turbine was used to validate the code extensively for the cases of monopile wind turbine, single-Dof platform motions and mutiple-Dof platform motions, respectively, where the thrusts, angle of attacks, and induced velocities are found to have good agreement with the literature. To conclude, the modified free wake vortex ring method proposed here is considered to be effective when solving the aerodynamic load around horizontal-axis wind turbines, on both fixed and floating support structures.