Floating offshore wind turbines (FOWTs) are subjected to platform motion induced by wind and wave loads. The oscillatory movement trigger vortex instabilities, modifying the wake structure and influencing the flow reaching downstream wind turbines. In this work, the wake of a FOWT is analyzed by means of numerical simulations and a comparison with linear stability theory. Two simplified models based on the stability of vortices are developed for all degrees of freedom of turbine motion. In our numerical simulations, the wind turbine blades are modeled as actuator lines and a spectral-element method with low dispersion and dissipation is employed to study the evolution of the perturbations. The turbine motion excites vortex instability modes predicted by the linear stability of helical vortices. The flow structures that are formed in the non-linear regime are a consequence of the growth of these modes and preserve some of the characteristics that can be explained and predicted by the linear theory. The number of vortices that interact and the growth rate of disturbances are well predicted by a simple stability model of a two-dimensional row of vortices. For all types of motion, the highest growth rate is observed when the frequency of motion is one and a half the frequency of rotation of the turbine that induces the out-of-phase vortex pairing mechanism. For lower frequencies of motion, several vortices coalesce to form large flow structures, which cause the high amplitude of oscillations in the streamwise velocities, which may increase fatigue or induce high amplitude motion on downstream turbines.

## I. INTRODUCTION

Offshore wind energy has been receiving an increasing amount of attention in recent years. The reason for this is the several advantages, it has in comparison with onshore wind turbines, such as stronger and more consistent wind conditions, the possibility to install larger turbines, the reduced impact of the noise generated by the turbines, and the mitigated harm it causes to wildlife. In order to make this kind of energy accessible and more affordable, a few technological issues need to be overcome. For example, for deep waters (typically above 60 m), the bottom-mounted configuration is no longer feasible and floating offshore wind turbines (FOWTs) have to be considered, where, inspired by naval architecture solutions, the turbine is usually mounted on a floating platform.

Since this platform is not bottom-mounted, it is allowed (under the constraint of the mooring system) to move under the action of the wind, sea waves, and sea currents. This movement may produce modifications to the structure of the wake of the turbine, possibly having important performance impacts on other downstream turbines. Some studies have been carried out in the past in that direction. For example, the performance of the wind turbine (such as power and thrust) was evaluated using imposed pitch^{1,2} and surge^{3} motion. Later, other studies^{4,5} performed numerical simulations in which the platform motion was not imposed but computed. More recently, the effect of the motion of turbines mounted on spar and submersible platforms on the generation of power was assessed through numerical simulations coupling the blade element momentum to calculate the blade loads, a control model for the turbine, a simplified hydrodynamics model to account for the platform motion, and LES (large eddy simulation) of the free wind.^{6} However, all those studies have paid little attention to the far-wake generated by the turbine, focusing only on near-turbine phenomena affecting the performance of that turbine. The predominant numerical model employed in those works was URANS (unsteady Reynolds-averaged Navier–Stokes) in which the flow close to the blades was resolved.

Vortex instabilities in the near wake of FOWTs were already observed by Sebastian and Lackner^{7} and Farrugia *et al.*^{8} Recently, a few studies^{9,10} have paid more attention to the issue of the far wake when platform movement occurs. To avoid the dissipation of traditional low-order computational fluid dynamics (CFD) methods, the wake was modeled by Lagrangian-based vortical particles,^{9} vortex rings,^{10} or a free vortex wake method.^{7,8} With these vortex-based models, they were able to provide evidence of possible patterns for the interaction of tip-vortices when the platform undergoes motion on each of its six solid degrees of freedom. However, their work focused more on the forces and performance of the turbine than on the wake, so the stability of the wake was not thoroughly explored.

Alternatively, several numerical studies have been carried out focusing on the stability of tip vortices of fixed wind turbines. Those studies are important for understanding how the wakes of several wind turbines interact with each other in a wind farm and how it affects its global performance. Earlier studies^{11,12} developed a theoretical framework to study the stability of helical vortices. The mutual inductance mode predicted by Widnall^{11} has been observed in experiments with models of ship propellers^{13} and wind turbines^{14} as the main mechanism of vortex interaction. The results of the stability theory have been confirmed quantitatively both experimentally^{15–18} and numerically.^{19,20} Later, those instabilities were studied numerically on a flow under different configurations such as shear^{21} and yaw misalignment.^{22} More recently, the interaction of vortices of two in-line turbines was studied.^{23} These recent numerical studies, focused on fixed wind turbines, relied on high-order numerical simulations of the Navier–Stokes equations for their superior resolution of the vortex structures. The perturbations were imposed as body forces on the blade tips, and the goal was to investigate the evolution and growth of instabilities. However, since the focus was not on the flow around the airfoil blades, the blades were not resolved but modeled with the actuator line method (ALM).^{24} The ALM has the benefit of not requiring the detailed mesh necessary to fully represent the blade geometry.^{25} A more faithful representation of the blades, using immersed boundary methods,^{26} observed the same tip vortex instabilities mechanisms. In a follow-up study,^{27} tip vortex instabilities were confirmed as having a major role in the wake recovery of a configuration of two in-line turbines.

Our recent work^{28} connects the study of wakes of heaving turbines and the stability of helical vortices. Using a high-order spectral element method code, we performed numerical simulations of a heaving turbine and qualitatively compared some of the flow features with the predictions of linear stability theory. In the present study, the work of Kleine *et al.*^{28} is extended to include other types of motion. Simulations of pitch and surge motions, in addition to heave, are performed. The details of the derivation of the stability analysis for moving turbines are greatly extended to provide more details and to include all the degrees of freedom: surge, heave, sway, pitch, yaw, and roll of the turbine. The growth rate and the modes induced by turbine motion from the numerical simulations are quantitatively compared to the predictions of the linear stability theory in the present work. A longer domain is used in this study to analyze the evolution in space of the flow structures created by the turbine motion and its impact on the flow in the region where a downstream turbine would be installed.

The present study focus on the stability of the wake of turbines under motion. Hence, the effect of motion on forces and performance is not discussed. The goal is to provide models based on the theory of hydrodynamic stability that could be applicable to several cases and to validate these models. Therefore, we present the model using non-dimensional and normalized parameters. In such a manner, the models are general and could be applicable to most turbines in different conditions. For completeness, some of the developments, descriptions, discussions, and conclusions of Kleine *et al.*^{28} are repeated in this work. This article is organized as follows: In Sec. II, the numerical method used for the simulations is described. A first-order approximation of the stability properties of the flow is presented in Sec. III. The numerical results are shown and compared to the stability theory in Sec. IV. Finally, the main conclusions are presented in Sec. V.

## II. NUMERICAL METHODS

### A. Numerical solver, domain, and boundary conditions

Nek5000, a spectral element method (SEM) code, is used to solve the three-dimensional Navier–Stokes equations in a fixed frame of reference. The spectral element method exhibits low dispersion and dissipation,^{29} which is relevant for stability calculations. In each spectral element, seventh-order Lagrange polynomials on Gauss–Lobatto–Legendre quadrature points are used for spatial discretization and a third-order implicit/explicit scheme is applied for temporal discretization. To stabilize numerical instabilities, filtering is applied.^{22}

All quantities are non-dimensionalized by the turbine radius *R* and the free-stream velocity $W\u221e$. The turbine is positioned at the origin of a cylindrical domain of radius *R _{rad}* = 5 [Fig. 1(a)]. The distances of the inlet and the outlet to the turbine are

*z*= 5 and

_{in}*z*= 22.10, respectively, for the longer domain used in the simulations reported in Sec. IV A. Around the center of the domain, the elementwise discretization in the streamwise direction is uniform between $\u22120.6<z<14.025$, with constant spacing $\Delta z=0.075$. This discretization is consistent with previous studies,

_{out}^{30}which evaluated the minimum grid spacing necessary for accurate simulations using actuator line methods, bearing in mind that in this work we are using seventh-order Lagrange polynomials as basis functions. A shorter domain with

*z*= 11.78 and uniform streamwise discretization between $\u22120.6<z<8.025$ is used for comparison of the growth rate and modes with the stability theory.

_{out}As can be seen in Figs. 1(b) and 1(c), the discretization is coarser in regions farther from the turbine, in order to reduce computational cost. Dirichlet boundary conditions with constant velocity $W\u221e=1$ are imposed on the inflow and lateral boundaries. Sheared and turbulent inflow are not considered in this study. At the outlet, the natural outflow boundary condition is imposed in conjunction with a sponge region, with a width $\Delta z=2.5$, that forces the *x* and *y* components of velocity to zero. Due to the nature of the outflow boundary condition, no forcing of the *z*-component is needed.

There is a confinement effect for this domain, given the relatively low value of *R _{rad}* = 5, which corresponds to a blockage ratio of $(R/Rrad)2=4%$. The results were compared to a simulation with a domain with

*R*= 10 (the domain and grid described in Ref. 23), with a blockage ratio of 1%. Differences on the order of 1%–2% in the streamwise velocity were observed in regions upstream of the turbine and outside the wake. Inside the wake, the differences are smaller. More importantly, the confinement effect did not affect significantly the position of the vortices and the qualitative aspects of the flow. Also, from a theoretical perspective, the influence of the confinement for

_{rad}*R*= 5 on the stability of the vortices is negligible. Therefore, we concluded that the effect of the confinement on the main results and conclusions is negligible.

_{rad}### B. Wind turbine modeling

The wind turbines are modeled using the actuator line method (ALM).^{24} In this method, the blades are represented by body forces calculated from airfoil data and local velocity. Considering the local velocity and the chord and twist distributions, the normal and tangential forces are calculated for $NACL$ points ($xn$) along with the blades. The discrete two-dimensional force vector ($F2D$) is projected on the domain (**x**) using the convolution of the force with a three-dimensional Gaussian kernel:

where *ε* is a smearing parameter. Following parametric studies,^{22} $\epsilon \u22483.5\Delta r$ (where $\Delta r$ is the averaged grid spacing) and each blade is discretized with 100 points.

The turbine is modeled after the turbine of the blind test.^{31} The chosen operation condition was the tip-speed ratio $\lambda =\Omega R/W\u221e=6$, which corresponds to the optimum performance for this turbine (Ω is the angular frequency of the rotor). Hub and tower were not included, since the focus of the work is the stability of the tip vortices. The code was extensively validated for this turbine, showing good agreement between numerical simulations and experimental results^{22,32,33} and was previously used in vortex stability studies^{21–23} with the same turbine and similar discretization and actuator line parameters.

The lift and drag coefficients of the 14% thick NREL S826 airfoil were obtained by Sarmast and Mikkelsen.^{34} For the computational simulation, a Reynolds number based on the radius of $ReR=W\u221eR/\nu =50\u2009000$ is used (where *ν* is the kinematic viscosity). However, since the Reynolds number used to interpolate the airfoil data from the corresponding lookup table can be independent of the CFD Reynolds number, a constant factor of 5.64 is applied to the local Reynolds number used to interpolate the airfoil lift and drag coefficients, so that the tip local Reynolds number (based on the chord at the tip, *c _{t}*) is on the order of $Rectip=W\u221e\lambda ct/\nu \u2248105$, matching the value from the experiments.

^{31}

The ALM has been shown to model very well the near wake, when compared both to experimental^{35} and simpler vortex methods.^{36} The main observed difference is in the vortex core size^{35} (see also Refs. 37 and 38 for an explanation of this difference), which does not have a first-order effect on the stability properties, since the main stability mechanism is inviscid. As can be seen in Refs. 11, 12, and 17, the effect of the core size on the stability and initial dynamics of the vortices is secondary, through the desingularization of the Biot–Savart law. For this reason, previous studies using the ALM had good agreement with theoretical predictions^{19–22} and it is considered an important tool to study the stability of tip vortices.^{25,39}

### C. Turbine motion

When a wind turbine is mounted on a floating platform, it is subjected to motions that result from the loads imposed by the wind on the turbine, the waves and sea currents on the platform, and the reaction and constraints of the mooring lines. Wind and waves are usually aligned, and the yaw control system of the turbine guarantees that the rotor faces the wind directly. The floating platforms that are usually employed for floating turbines, namely spar-buoy, semi-submersible platform, and barge, have a high degree of spatial symmetry with respect to the vertical axis. All these aspects make the movement of the turbine occur mostly at the vertical plane-parallel to the wind. In other words, the most relevant motions of a floating offshore wind turbine are heave (translation in the *y*-direction), surge (translation in the *z*-direction), and pitch (rotation around the *x*-axis), see Fig. 2.

In Ref. 28, we considered heave motion only. To assess the amplitude and frequency of heave, we performed simulations using the software FAST v8^{41} for the NREL 5 MW turbine^{42} mounted on the OC4 semi-submersible platform,^{43} subject to a 10.3 m/s wind speed and waves corresponding to a JONSWAP wave spectrum with significant height (Hs) of 1.5 m peak spectral period (Tp) of 8 s. These environmental conditions are typical of the Southeast Brazilian continental shelf. Simulations were run considering water depths from 500 to 1000 m, with catenary mooring lines designed specifically for each of the depths. We observed heave amplitudes of approximately 1% of the radius (2% peak-to-peak) and a frequency of about half of the rotation of the blades. The distance between the center of the turbine and the axis of rotation of pitch was estimated to be between $1.4R$ and $1.5R$. So, in the current study, we employ these values as starting points for heave, surge, and pitch motions, and investigate how variations in frequency and type of motion affect the wakes produced.

It should be noted that the turbine used to estimate the amplitude and frequency of motion is not the same used in the numerical simulations (described in Sec. II B). Nevertheless, we assume that the order of magnitude of the non-dimensional parameters of the motion is independent of the turbine.

In the current simulations, the motion and the tip speed ratio of the turbine are prescribed. In other words, the turbine motion does not change the tip speed ratio and the non-steady fluid-dynamic effects do not affect the amplitude and frequency of motion. The heave and surge motions were implemented as the movement of the center $(xt,yt,zt)$ of the actuator lines,

where *ω* is the angular frequency of motion and the subscripts _{y} and _{z} represent heave and surge, respectively. We define $\omega *=\omega /\Omega $ as the normalized angular frequency.

Pitch motion, defined as rotation along the *x*-axis, is implemented by changing both the position of the center of the turbine and the angle of the turbine with respect to the freestream. The distance of the center of the turbine to the axis of pitch rotation in the *z*-direction is considered to be negligible compared to the distance in the *y*-direction and the radius of the turbine. The center of pitch rotation is assumed to be at position *R _{p}*, at the

*y*-axis. The equation that defines the position of the tip of the blades $(xb,yb,zb)$ under pitch motion (subscript

_{p}) is

where $\varphi n=2\pi n/Nb$ is the initial azimuthal position of each blade *n* (*n *=* *0 to $Nb\u22121$, where *N _{b}* = 3 is the number of blades) and

*A*is the angular amplitude of motion, measured in radians. According to the estimates performed, the position of the axis of rotation of pitch is rounded to $Rp=\u22121.5R$.

_{p}## III. STABILITY ANALYSIS

In this section, a linear stability analysis of the flow is presented based on a first-order approximation of the perturbations imposed by the motion of the turbine. This was first presented for the heave motion in Ref. 28 and is extended here for the other degrees of freedom. Initially, first-order approximations of the perturbations on the helical vortex system caused by the turbine motion are obtained for all types of motion. The perturbations are then compared to the unstable modes predicted by the classical results of the linear stability of a two-dimensional infinite row of vortices^{44} and perfect helical vortices.^{11,12}

### A. Perturbations imposed by turbine motion

#### 1. Reference system for the stability analysis

For a fixed-bottom turbine, the equation that defines the tip of the blade that rotates in the positive azimuthal direction is

where $(xt0,yt0,zt0)=(0,0,0)$ is assumed without loss of generality.

Previous studies identified that the tip vortices are mainly advected by the flow in the streamwise direction;^{21} hence, a uniform advection velocity *w _{c}* in the streamwise direction is considered. The position $(xv0\u2032,yv0\u2032,zv0\u2032)$ of the tip vortices emitted by the blades of Eq. (5) is

where *τ*_{0} is the time at which the vortex particle was emitted.

Imposing sway, heave, and surge (subscripts _{x}, _{y}, and _{z}, respectively) to the turbine, the position of the tip of the blades is

The turbine motion imposes perturbations to the position of the vortices that grow in time and space. In order to study the evolution of these perturbations, we initially do not consider this growth, assuming only the advection of the vortices. The position $(xv\u2032,yv\u2032,zv\u2032)$ of the tip vortices with turbine motion is then

In the general case, there would be a phase difference between the motions and they could be composed of several frequencies. However, as a linear approximation, these motions do not interact, so the phase difference is not relevant for the present analysis. Therefore, in favor of a more concise notation, the phase is omitted and only one frequency for each motion is analyzed.

The helical structure of the vortices of Eq. (6) changes over time. The vortical particles translate in *z*, making the helical structure translate or rotate, depending on the point of view (for further discussion about this topic, see Okulov and Sørensen^{45}). In order to compare with classical stability results,^{11,12,44} a frame of reference in which the vortices are static should be chosen. There are infinite choices for the pair of streamwise and azimuthal velocities that make the helical structure static, but different reference systems give different frequencies in the stability analysis.^{46} We then choose the reference system that follows the vortical particles, which is the same as that employed by Widnall,^{11} Gupta and Loewy,^{12} and Lamb.^{44} Hence, considering the assumption that the vortices are simply advected by a constant streamwise velocity, the reference system is a translating frame of reference with velocity *w _{c}*. Taking the origin in

*z*of the system as the position where $\tau 0=0$, Eq. (8) becomes

Hence *τ*_{0} can be thought of as the parameter that defines the helical vortices.

It is worth noting that the helical vortices defined by Eq. (9) are left-handed helices, just as the vortices generated in our simulations. However, most of the stability studies consider right-handed helices, including our implementation of the method of Gupta and Loewy.^{12} Also, for positive *w _{c}*, the values of

*z*and

_{v}*τ*

_{0}have different signs, which indicates that the perturbation would have a different sign when written in the function of

*z*: $A\u2009sin\u2009(\omega \tau 0)=\u2212A\u2009sin\u2009(\omega zv0/wc)$. Both of these characteristics are undesirable for the comparison of the perturbations with the stability theory. Hence, without loss of generality, we consider a modification to the vortex system, so the position is given by

which has the desirable properties of being formed by right-handed helices and having perturbations that have a sinusoidal component in *z* with the same sign of the new parameter *τ*, providing greater similarity to the geometry and perturbations of classical stability studies. The stability properties for the geometry of Eq. (10) are the same as Eq. (9), but the comparison with classical results can be performed without further modification to the reference system by using Eq. (10). On the other hand, a comparison to the simulations should be made considering different reference systems.

#### 2. First-order approximation of perturbation imposed by the sway, heave, and surge

Using the reference system of Eq. (10), the position of the tip vortices of right-handed helices without any turbine motion would be

Thus, the perturbation directly induced by the turbine motion is $(\delta xv,\delta yv,\delta zv)=(Ax\u2009sin\u2009(\omega x\tau ),Ay\u2009sin\u2009(\omega y\tau ),Az\u2009sin\u2009(\omega z\tau ))$, in cartesian coordinates. For $Ax,Ay\u226aR$, first-order approximations of the perturbation in the radial direction and azimuthal directions are

The order of magnitude of $\delta rv$ and $R\delta \phi v$ is the same. However, the perturbation in the azimuthal direction is not usually as relevant as the other components for the stability of helical vortices of low pitch ($h/R\u226a1\u21d4$ $wc/(Nb\Omega R)\u226a1$), so, as an approximation, it is neglected. Defining $\theta =\Omega \tau +\varphi n=\Omega zv0/wc+\varphi n$, the first-order approximation of the perturbation in the radial and streamwise directions is defined for the helices in the function of the azimuthal angle as

This initial perturbation induced by the motion grows in time and space, as the vortices move downstream. The growth rate can be estimated using analytical stability models. The applied stability models consider infinite uniform vortices, which is a limitation of the model developed here. Hence, the region of wake expansion is not considered but wake expansion in the region of uniform helices could be partially taken into account by using a different value of *R*.

Another limitation is that only perturbations in position are considered. Other effects may be also relevant for some types of motion. For example, for heave and sway motion, the change in position of the blockage effect of the turbine might impose other perturbations to the wake that are not considered here. The modification of the forces at the blades and, consequently, the circulation is also not accounted. If these effects also have the form of Eq. (14) (pure sinusoidal oscillation for surge and sinusoidal with $sin\u2009\theta $ and $cos\u2009\theta $ envelopes for heave and sway), then the analysis performed here would also be applicable.

#### 3. First-order approximation of perturbation imposed by pitch

From Eq. (4), applying the change in the reference system discussed in Sec. III A 1, the equation that defines the position of the tip vortices under pitch motion becomes

where the approximation comes by assuming that the maximum pitch angle, *A _{p}*, is low. We can directly notice that the motion in the

*y*-direction is of second order. For low amplitudes, this can be neglected and the first-order approximation of the perturbation due to pitch motion becomes

which, written in the radial and streamwise directions, is

Hence, the perturbations imposed by the pitch motion are equivalent to perturbations imposed by the surge motion. For $Rp\u226bR$, this corresponds to a surge motion with amplitude $RpAp$. For $Rp\u226aR$, this corresponds to a motion in the *z*-direction composed of the component $sin\u2009(\omega p*(\theta \u2212\varphi n))\u2009sin\u2009\theta $ with amplitude $RAp$, which has the same sinusoidal envelope as the heave component in Eq. (14).

The first-order approximations of the perturbations imposed by yaw and roll are shown in Appendix A. The change in the direction of the thrust vector in the case of yaw or pitch that causes wake deflection (see Sec. 4.2 of Burton *et al.*^{47}) could potentially also induce perturbations to the vortices, including in the radial direction. These effects are not considered in this small-perturbation model.

As can be observed, all the components of the perturbations induced by pitch, yaw, and roll have the same format of an amplitude multiplied by $(sin\u2009(\omega *(\theta \u2212\varphi n))\u2009sin\u2009\theta ),\u2009(sin\u2009(\omega *(\theta \u2212\varphi n))\u2009cos\u2009\theta )$, or $(sin\u2009(\omega *(\theta \u2212\varphi n)))$. Conceptually, for the purposes of the analysis of Secs. III B and III C, these components are equivalent to the perturbations of Eq. (14). Hence, the detailed analyses are restricted to the case of surge and heave. In the end, a parallel to the case of other kinds of motion is made.

### B. Stability of a 2D row of vortices

Several studies made the connection between the stability of helical vortices and a two-dimensional row of vortices,^{15,17,18,20} showing that the scaled growth rate of long-wave and mutual inductance instabilities collapses to the theoretical value for an infinite row of vortices calculated by Lamb^{44} (see also Refs. 48 and 49). The growth rate *σ* is given by^{48,49}

where *h* is the distance between neighboring vortices, Γ is the circulation and *p* is a dimensionless subharmonic wavenumber (see Fig. 3). The maximum growth rate occurs for *p *=* *1/2, which corresponds to the vortex pairing where neighboring vortices are out-of-phase.

The growth rate when scaled using the distance of neighboring vortices and the circulation, $\sigma \u0303=\sigma (2h2/\Gamma )$, has a maximum value of $\pi /2$. The growth rate of the vortex pairing mode collapses to $\pi /2$ even for non-uniform helices created by turbines under sheared inflow or yawed if scaled by local properties, as shown by Kleusberg *et al.*^{21,22} The two modes $(\delta rm2d1,\delta zm2d1)$ and $(\delta rm2d2,\delta zm2d2)$ predicted by the linear stability theory for each *p* are^{48,49}

and

where *i* is the complex unit. These modes are illustrated in Fig. 3. One of these modes is exponentially damped and the other grows exponentially, depending on the sign of the circulation and the reference system. For the reference system of Fig. 3, which corresponds to the reference system used in this work for a wind turbine, the mode $(\delta rm2d1,\delta zm2d1)$ is unstable.

In order to compare to a row of two-dimensional vortices, we take a cross-section of the helical vortices, taking the values of $\theta =\theta 0+2\pi j$, where *j* is an integer and *θ*_{0} is chosen to define a cross-section that contains the row of vortices ($0\u2264\theta 0<2\pi $). This is a slightly different approach from that in Quaranta *et al.*,^{17} where the minimum distance between neighboring vortices was used. However, for helices of low pitch, the difference between the minimum distance and the distance defined by a cross-section is of second order. From Eq. (14), the displacement of the vortices due to surge, heave, and sway along the cross-section defined by *θ*_{0} is

where the term $Nbj\u2212n$, which can be any integer, has been reinterpreted as the index *m* of each 2D vortex (Fig. 3).

Projecting the terms of Eq. (20) onto the modes defined in Eq. (19), the perturbation $(\delta rv2d,\delta zv2d)$ can be decomposed into the components $(\delta rv2d1,\delta zv2d1)$ and $(\delta rv2d2,\delta zv2d2)$,

and

where the notation $c.c.$ indicates the complex conjugate of the previous term. For the current reference system, mode 1 grows exponentially, and mode 2 decays. Hence, after some distance from the turbine, mode $(\delta rv2d1,\delta zv2d1)$ tends to dominate. It can be noted that there is a symmetry between the perturbation in the component $\delta rv2d$ and $\delta zv2d$: a perturbation in $\delta rv2d$ induces a component in $\delta zv2d$ and vice versa because the eigenvectors have equal components in *r* and *z*.

Qualitatively, the three terms of Eq. (21) are similar. The terms $(Ay\u2009sin\u2009\theta 0),\u2009(Ax\u2009cos\u2009\theta 0)$ and *A _{z}* are the amplitudes of the perturbation induced by each motion. The main difference is that the amplitude for the heave and sway motion is modulated by the angle

*θ*

_{0}, while the amplitude of the surge motion is uniform in the azimuthal direction. This implies that the heave motion has a maximum amplitude at angles $\theta 0=\xb1\pi /2$ and does not induce instability at angles $\theta 0=0$ and

*π*. Due to the three-dimensionality of the problem and non-linear dynamics, the helical vortices created by a heaving turbine would also interact at angles $\theta 0=0$ and

*π*. Nevertheless, the vortices at $\xb1\pi /2$ would have noticeable larger amplitudes.

From the terms of Eq. (21), it follows directly that the terms $i\omega *\theta 0$ correspond to a phase shift, and, most importantly, $\omega */Nb$ is the dimensionless subharmonic wavenumber *p* for each motion type. Applying the scaling $p=\omega */Nb$ to Eq. (18), we obtain the theoretical growth rate for a row of vortices shown in Fig. 4. Equation (18), in principle, is only valid for $0\u2264p\u22641$. For values of $\omega */Nb$ that are outside this interval, however, we use $p=\omega */Nb+j$ where *j* is any integer that makes $0\u2264p\u22641$ (including for negative values of $\omega *$).

As a linear approximation, this growth rate is applicable to all motion types, including pitch, yaw, and roll, following the discussion of Sec. III A 3. Because this is a linear analysis both in the approximation of the perturbation and in the stability theory, the perturbations grow independently of each other. The effect of non-linear interaction of the perturbations imposed by different motions is not considered in this theory.

### C. Stability of helical vortices

Regarding the stability of the three-dimensional helical vortices, the knowledge of the wavenumbers induced by the perturbations is essential to compare to classical results from the analytical linear stability.^{11,12} As discussed in Kleusberg *et al.*^{21} and Kleine *et al.*,^{23} for perturbations imposed as body forces at the tip of the blades, the angular frequency is related to the wavenumber according to the approximate relation $\omega =k\Omega $. However, as it is going to be shown, in the case of perturbations imposed by turbine motion, other wavenumbers can be induced for a certain frequency, depending on the motion of the turbine.

In the two-dimensional case, the perturbation in one direction, when projected onto the eigenvectors, creates a perturbation in the other component with equal magnitude (Sec. III B). This means that the heave and sway motions induce a component in the streamwise direction and the surge motion induces a component in the radial direction, so that $\delta z=\delta r$. In the case of three-dimensional vortices, the most amplified modes predicted by the stability theory also have $\delta z\u2248\delta r$, so that the projection of the perturbations onto the eigenvectors would create perturbations in the other component. Hence, we can assume that the perturbation given by Eq. (14) is composed of two terms

and

where $(\delta rv1,\delta zv1)$ grows and $(\delta rv2,\delta zv2)$ decays in time and space (for the current choice of a reference system).

Using trigonometric identities, Eq. (23) can be decomposed in

and analogously for Eq. (24). The wavenumber *k* of each perturbation is given by the term that multiplies *θ*, while the terms $\xb1\varphi n$ are phase shifts between the different vortices *n*. Thus, a sway or heave motion (*x* and *y*-components) with frequency $\omega *$ creates perturbations with wavenumbers $k=\omega *\u22121$ and $k=\omega *+1$. The surge motion (*z*-component), on the other hand, creates perturbations with wavenumber $k=\omega *$. The surge motion is similar to previous studies where body forces were used as perturbations at the tip, which arrived at the same approximated relation $k=\omega *$.^{21}

Analogously to the row of vortices, when the real perturbation is projected onto the eigenvectors, a pair of complex conjugate modes are excited, see Eqs. (21) and (22). Therefore, the heave and sway motions excite modes with wavenumbers $k=\xb1(\omega *\u22121)$ and $k=\xb1(\omega *+1)$ and the surge motion excite modes with wavenumbers $k=\xb1\omega *$. In Fig. 5, the growth rate predicted from the eigenvalues using the method of Gupta and Loewy^{12} is shown as a function of the wavenumber. All types of modes excite stable and unstable modes (eigenvalues with positive growth *σ* are the unstable modes). A heave motion with frequency $\omega *$ excites eight eigenmodes:

Two unstable eigenmodes with $k=\omega *\u22121$ and $k=\omega *+1$ and their complex conjugates [$k=\u2212(\omega *\u22121)$ and $k=\u2212(\omega *+1)$];

Two stable eigenmodes with $k=\omega *\u22121$ and $k=\omega *+1$ and their complex conjugates [$k=\u2212(\omega *\u22121)$ and $k=\u2212(\omega *+1)$];

while the surge motion excites four eigenmodes: stable and unstable modes with $k=\xb1\omega *$. As can be seen in Eq. (25), the sway motion excites the same modes of the heave motion, only with a phase difference of $\pi /2$. The example of heave and sway motions for $\omega *=1.5$ is shown as an example in Fig. 5. The modes with $k=\omega *\u22121$ and $k=\omega *+1$ are referred to as the primary modes and the modes $k=\u2212(\omega *\u22121)$ and $k=\u2212(\omega *+1)$ as their complex conjugate.

The stability diagram of Figs. 5 and 6 is usually only shown in the literature for positive values of *k*, since it is symmetric. However, since the primary mode $k=\omega *\u22121=\u22120.5$ for $\omega *=0.5$ falls on the negative side, it is useful to look at the two-sided diagram. In the same sense, to say, that a real (without imaginary part) perturbation excites wavenumber $k=\omega *\u22121$ is equivalent to saying that the perturbation excites the conjugate pair with wavenumber $k=\xb1(\omega *\u22121)$.

Even knowing the wavenumber of a perturbation is not sufficient to know which unstable mode is excited. For most wavenumbers, there are three corresponding unstable eigenvalues, as can be seen in Fig. 5 (for higher wavenumbers or higher values of pitch, this may be different, see examples in Ref. 18).

Some of the eigenvalues predicted by the stability theory correspond to modes that have a phase correlation between all vortices at *z *=* *0, while others have a phase difference of $\xb12\pi /Nb$ (see Sarmast *et al.*,^{20} Quaranta *et al.*^{18} for complementary discussion about these modes). This is due to the symmetry of the stability problem regarding the *N _{b}* helical vortices in the azimuthal direction. From the surge motion term of Eq. (25), we note that there is no phase shift between the vortices (at the same streamwise position, the amplitude of the perturbation is the same for all helices). Hence, the eigenmodes that have a phase correlation are those excited by the surge motion, as shown in Fig. 6(a).

The heave and sway motion have terms with a phase shift of exactly $\xb1\varphi n=2\pi n/Nb$ between the vortices, see Eq. (25), thus they excite the modes with a phase difference, shown in Fig. 6(b). In that figure, there are two eigenvalues for each wavenumber. This can be explained by the two possible values of the phase difference between the vortices. One of the branches has a phase difference $2\pi /Nb$ and the other a phase difference of $\u22122\pi /Nb$ [the branch without phase difference was already extracted from the figure and shown in Fig. 6(a)]. Looking at Eq. (25), it is expected that each primary mode excited by one frequency belongs to a different branch because one term has a phase difference $+\varphi n$ and the other $\u2212\varphi n$. This is exactly what is observed: the mode with $k=\omega *\u22121$ belongs to one branch (related to $\u2212\varphi n$) and the mode with $k=\omega *+1$ belongs to the other branch (related to $+\varphi n$). It is easy to see that the complex conjugate of the primary mode belongs to the opposite branch.

As can be seen in Figs. 7 and 8, the perturbations created by the heave motion [Eq. (25)] are almost exactly the same as the eigenmodes predicted by the stability theory of Gupta and Loewy.^{12} In order to obtain real eigenvectors, each complex eigenvector was summed with its complex conjugate, applying the adequate phase. Even though the figures were obtained with completely different methods, they look identical. For all other cases, similar agreement was found. Perturbations for sway motion can also be reconstructed from the same eigenmodes, using appropriate values of phase. For surge motion (not shown here), the same agreement was observed when compared to the modes of Fig. 6(a), when considering the modes with $k=\xb1\omega *$. This indicates that the perturbations imposed by the turbine motion agree almost identically to eigenvectors of the stability theory; thus, only these modes would be excited.

In Ref. 28, we provided a rule-of-thumb to identify which eigenvalue corresponds to each heave frequency when there are multiple eigenvalues for the same wavenumber. The eigenvector that is observed to be identical to the perturbation is that whose primary eigenvalue in the parabolas of Fig. 6(b) has the same relative position in the parabola of a 2D row of vortices (Fig. 4). For example, for $\omega *=0.5$,

There are two eigenvalues for $k=\omega *\u22121=\u22120.5$, one close to $\sigma \u0303/(\pi /2)\u22480.55$ and another close to $\sigma \u0303/(\pi /2)\u22481.0$. The eigenvector that is observed to be identical to the perturbation (see Fig. 7) is that whose eigenvalue has a growth rate $\sigma \u0303/(\pi /2)\u22480.55$, because its position in relation to the parabola between $\u22121<k<2$ is the same of $\omega *=0.5$ in Fig. 4 (diamond-shaped symbol). [The two eigenmodes for the pair $k=\xb10.5$ with phase shift are shown in Figs. 7(d) and 8(d), while the mode without phase shift is equivalent to Fig. 14(e)].

There are two eigenvalues for $k=\omega *+1=1.5$, both close to $\sigma \u0303/(\pi /2)\u22480.55$ (it is close to the intersection of parabolas between $\u22121<k<2$ and between $1<k<4$). The eigenvector that is observed to be identical to the perturbation (see Fig. 7) is that whose eigenvalue belongs to the parabola between $1<k<4$, because the position of the eigenvalue in relation to this parabola is the same of $\omega *=0.5$ in Fig. 4.

The term “parabola” was freely used here, but we note that the curves do not represent exact parabolas, except for Fig. 4, which comes from Eq. (18). The discussion presented here provides an explanation for this rule-of-thumb: the branch that corresponds to $\u2212\varphi n$ is similar to the branch without phase difference [Fig. 6(a)], but displaced one unit to the negative direction in the *k*-axis ($k=\omega *\u22121$), while the branch that corresponds to $+\varphi n$ is similar to the branch without phase difference, but displaced one unit to the positive direction in the *k*-axis ($k=\omega *+1$).

Due to this connection between the excited modes and the stability of the 2D row of vortices, the results of stability of the full three-dimensional helical system might not be needed to estimate the growth rate of perturbations imposed by turbine motion. It is expected that the growth rate can be predicted from the theory of infinite row of two-dimensional vortices, within the restrictions discussed in Refs. 17 and 18.

The surge motion $\omega *=1.5$ corresponds to the vortex pairing mechanism that has the highest growth rate $\sigma \u0303\u2248\pi /2$, which has been discussed in previous studies that excite the vortices using body forces.^{19–21} Additionally, we note that all motions with an angular frequency $\omega *=1.5$ would also excite the vortex pairing mechanism that has the highest growth rate, even if the wavenumbers are *k *=* *0.5 and *k *=* *2.5. Figure 8 shows how the wavenumbers *k *=* *0.5 and *k *=* *2.5 correspond to the vortex pairing, in which neighboring vortices are out-of-phase in relation to each other. The mechanism for which the vortex pairing is excited by wavenumbers different than $k=Nb(j+1/2)$ (where *j* is an integer) has been discussed in Refs. 18, 20, and 46.

We emphasize here, however, that the only frequencies of motion of a turbine that excite the vortex pairing with $\sigma \u0303\u2248\pi /2$ are given by $\omega *=Nb(j+1/2)$, even if other wavenumbers are excited. This is valid not only for surge and heave but for all types of motion where the turbine moves as one. The fact that the turbine moves as one restricts the phase between helices to only two options: either there is no phase shift and $k=\omega *$ (e.g., surge) or the phase shift is $\xb12\pi /3$ and $k=\omega *\xb11$ (e.g., heave). This means, for example, that a frequency of $\omega *=0.5$ would never excite the highest growth mode $\sigma \u0303\u2248\pi /2$ for *k *=* *0.5, as seen in Fig. 6(b). If the excitation is given by other forms, such as actuation on the blades or inflow turbulence, this restriction does not exist and vortex pairing may be induced by other frequencies.

The excitation of wavenumbers different from the normalized angular frequency for heave and sway is due to the terms $sin\u2009\theta $ and $cos\u2009\theta $ in Eq. (14). These terms can be thought of as coming from the decomposition of the oscillation of *y* (for heave) or *z* (for sway) in the radial direction, creating the phase shift between the perturbations on each vortex. Because surge does not have such terms, the wavenumber is equal to the normalized angular frequency. Hence, roll motion, Eq. (A3), would also create wavenumbers given by $k=\xb1(\omega *\xb11)$. Pitch, Eq. (17), and yaw, Eq. (A1), have components with and without terms $sin\u2009\theta $ and $cos\u2009\theta $. Therefore, these motions would induce perturbations composed of several wavenumbers: the terms $RAp$ and $RAyw$ induce wavenumbers $k=\xb1(\omega *\u22121)$ and $k=\xb1(\omega *+1)$ and the terms $RpAp$ and $RywAyw$ induce wavenumber $k=\xb1\omega *$.

A limitation of this model is the approximation that the perturbations are only advected in the streamwise direction. If the wake rotation is relevant, then the perturbations are also advected in the azimuthal direction, which would modify the wavenumbers induced by a certain frequency. Also, the model only considers mutual inductance and long-wave instabilities, which can be calculated from the theory of Gupta and Loewy.^{12} For very large wavenumbers the long-wave assumption is not valid.^{17} Short-wave instabilities^{50–53} are not captured by the model.

## IV. RESULTS

### A. Direct numerical simulations

Several cases with different types of motion and parameters were simulated, according to Table I. The reference case has an amplitude $A/R=1%$ and angular frequency $\omega =0.5\Omega $, according to estimates detailed in Sec. II C. The amplitude of pitch was defined so that the maximum displacement in *z* at the center of the turbine is 1% of the radius, assuming the linear approximation ($1%R/Rp=0.667%$). The effect of the increase in amplitude is to advance the onset of the instability, bringing it closer to the turbine, as discussed in Ref. 28, similar to the increase in the amplitude of body forces or turbulence level shown in previous works.^{19–21} Since the effect of the amplitude is well understood, in this work, we focus on the effect of the frequency and type of motion. A different amplitude, of $A/R=0.1%$, is only used in Sec. IV B, in order to compare the growth rate to analytical stability predictions.

Angular frequency $\omega *=\omega /\Omega $ . | Type of motion . | Amplitude A/R
. |
---|---|---|

0.5 | Heave^{a} and surge^{a} | 1% |

1 | Heave | 1% |

1.5 | Heave^{a} and surge^{a} | 1% |

2.5 | Heave | 1% |

0.5 | Pitch | 0.667% |

Angular frequency $\omega *=\omega /\Omega $ . | Type of motion . | Amplitude A/R
. |
---|---|---|

0.5 | Heave^{a} and surge^{a} | 1% |

1 | Heave | 1% |

1.5 | Heave^{a} and surge^{a} | 1% |

2.5 | Heave | 1% |

0.5 | Pitch | 0.667% |

^{a}

Cases chosen to be quantitatively compared to stability analysis, with an amplitude of $A/R=0.1%$, in Sec. IV B.

The change in the frequency of motion highlights a noteworthy effect. For the lower frequency, $\omega *=0.5$, larger flow structures are created inside the wake, as shown in Figs. 9 and 10. Observing the vortex interaction in Fig. 9, it is possible to notice that these flow structures are created by the merging of the vortices. The formation of similar structures for lower frequencies can also be noted in the simulations reported in Ref. 9. However, in their case, the lower frequency was accompanied by a modification like the motion, from translation (heave, sway, and surge) to rotation (yaw, pitch, and roll), making it difficult to isolate the effects. We show here that the existence of flow structures is independent of the type of motion, occurring for heave, surge, and pitch [Figs. 9(a), 9(e), and 9(g) and 10(a), 10(e), and 10(g)].

Interestingly, large flow structures, very similar to those observed for $\omega *=0.5$, can also be noted for the higher frequency, $\omega *=2.5$. According to the stability theory for a row of vortices, this is expected. Every frequency that can be written as $\omega *\u2212jNb$, with *j* integer, is indistinguishable from $\omega *$ according to the theory presented in Sec. III B. It follows that the configuration of a row of vortices for $\omega *=2.5$ is equal to $\omega *=\u22120.5$, which is the mode complex conjugate of $\omega *=0.5$. Thus, the number of vortices that interact directly is predicted to be the same for $\omega *=2.5$ and $\omega *=0.5$.

The number of vortices that interact with each other can be calculated from the maximum of $1/p$ and $1/(1\u2212p)$ from the relationship $p=\omega */Nb\u2212j$. The mechanism for $\omega *=0.5$ is illustrated in Fig. 3: the vortices $\u22121\u2264m\u22641$ would move to the right and the vortices $2\u2264m\u22644$ would move to the left, forming a group of six vortices that interact directly with each other. The values of Table II agree well with Fig. 9. For example, regardless of the type of motion, in the numerical simulations for $\omega *=0.5$, six vortices interact with each other and coalesce to form the large flow structure discussed in Sec. IV A. For $\omega *=1.5$, two vortices interact with each other and the vortex pairing mechanism is excited, as expected. The type of motion does not influence the number of vortices that interact, as long as the normalized frequency is the same, as predicted by the theory. Even the case of $\omega *=2.5$, which has a more complex dynamics, seems to follow this rule. For this case, the number of vortices that interact is difficult to count; however, the distance between the large flow structures is approximately 6*h*, as predicted.

$\omega *$ . | $p=\omega *Nb$ . | $1p$ . | $11\u2212p$ . | Number of interacting vortices . |
---|---|---|---|---|

0.5 | 1/6 | 6 | 6/5 | 6 |

1.0 | 1/3 | 3 | 3/2 | 3 |

1.5 | 1/2 | 2 | 2 | 2 |

2.5 | 5/6 | 6/5 | 6 | 6 |

$\omega *$ . | $p=\omega *Nb$ . | $1p$ . | $11\u2212p$ . | Number of interacting vortices . |
---|---|---|---|---|

0.5 | 1/6 | 6 | 6/5 | 6 |

1.0 | 1/3 | 3 | 3/2 | 3 |

1.5 | 1/2 | 2 | 2 | 2 |

2.5 | 5/6 | 6/5 | 6 | 6 |

The shape of the large flow structures is influenced by the type of motion. For $\omega *=0.5$, the surge motion creates flow structures that are coherent in the azimuthal direction while the heave creates flow structures that are tilted. The explanation for this difference is discussed in Sec. IV B. The pitch motion is more similar to surge, which is expected since $Rp>R$. Nonetheless, the structures created by pitch are not as coherent in the azimuthal direction as those of surge, showing an effect of the $(sin\u2009\theta )$ term, similar to heave.

The streamwise flow inside the wake experiences higher fluctuations due to these flow structures. These fluctuations might affect downstream turbines, increasing unsteady loads that may worsen the fatigue. In the case of a downstream FOWT, it may increase the amplitude of its motion. It is even possible to imagine an extreme case, where the unsteady flow and the sea waves add to each other with the same frequency, creating very high amplitudes. It is especially worth investigating because the frequency of the flow structures created for $\omega *=0.5$ is the same as the upstream heave motion, caused by sea waves.

In order to estimate the effect of the unsteadiness of the streamwise velocity on a downstream turbine, the parameters $Fz*$ and $Ty*$ were defined based on integrals of the streamwise dynamic pressure. $Fz*$ is an indicator of the streamwise force that a fictitious disk, with the same radius of the turbine (*R *=* *1) centered at $(xc,yc,zc)$, would be subjected to. It is the streamwise dynamic pressure integrated over the area of the disk, normalized by the equivalent value at infinity,

where *S* is the area of the disk and *ρ* is the fluid density. $Ty*$ is an indicator of the torque around the negative *y* direction (torque around the tower, that tends to yaw a turbine), defined as

where the normalization is performed considering an extreme case where the velocity on the semicircle of $x\u2212xc<0$ is zero and the velocity on the semicircle of $x\u2212xc>0$ [area $S/2$ and centroid at $4R/(3\pi )$] is constant and equal to $W\u221e$.

Four positions for the disk were chosen for this calculation. The disks were placed at the streamwise positions *z _{c}* = 6 and

*z*= 12, in positions completely inside the wake, $(xc,yc,zc)=(0,0,6)$ and $(xc,yc,zc)=(0,0,12)$, and positions partially inside the wake, $(xc,yc,zc)=(1,0,6)$ and $(xc,yc,zc)=(1,0,12)$. The values of the average in time of $Fz*$ and $Ty*$ are shown in Fig. 11, with error bars indicating their root mean square (RMS). In Table III, the ratio of RMS to the average of $Fz*$ is shown. The values of Table III deviate from the values found in Ref. 28 because there was a programming error in the code that calculated the RMS of that work.

_{c}Disk position $(xc,yc,zc)$ . | Heave ($\omega /\Omega $) . | Surge ($\omega /\Omega $) . | Pitch ($\omega /\Omega $) . | ||||
---|---|---|---|---|---|---|---|

0.5 . | 1 . | 1.5 . | 2.5 . | 0.5 . | 1.5 . | 0.5 . | |

(0,0,6) | 2.3% | 1.1% | 0.5% | 6.1% | 21.0% | 2.1% | 20.6% |

(1,0,6) | 1.6% | 0.4% | 0.6% | 1.8% | 0.5% | 1.0% | 1.6% |

(0,0,12) | 1.5% | 0.8% | 0.8% | 1.5% | 8.0% | 1.0% | 9.1% |

(1,0,12) | 0.8% | 0.5% | 0.3% | 0.9% | 1.1% | 0.4% | 0.5% |

Disk position $(xc,yc,zc)$ . | Heave ($\omega /\Omega $) . | Surge ($\omega /\Omega $) . | Pitch ($\omega /\Omega $) . | ||||
---|---|---|---|---|---|---|---|

0.5 . | 1 . | 1.5 . | 2.5 . | 0.5 . | 1.5 . | 0.5 . | |

(0,0,6) | 2.3% | 1.1% | 0.5% | 6.1% | 21.0% | 2.1% | 20.6% |

(1,0,6) | 1.6% | 0.4% | 0.6% | 1.8% | 0.5% | 1.0% | 1.6% |

(0,0,12) | 1.5% | 0.8% | 0.8% | 1.5% | 8.0% | 1.0% | 9.1% |

(1,0,12) | 0.8% | 0.5% | 0.3% | 0.9% | 1.1% | 0.4% | 0.5% |

As can be observed in Table III and Fig. 11, in general, the RMS of the integrals of the dynamic pressure is higher for angular frequencies $\omega *=0.5$ and $\omega *=2.5$, both for heave and surge. This indicates that these frequencies might be more dangerous for downstream turbines. As expected, a great oscillation in $Ty*$ is observed for a disk partially in the wake [$(xc,yc)=(1,0)$]. Surge and pitch cause great oscillation in $Fz*$ for a turbine completely inside the wake, because the flow structures are coherent in the azimuthal direction. This would cause high amplitudes in the forces applied to the tower and mooring cables of a downstream turbine. If these effects occur in real-world situations, it might be necessary to employ motion damping, wake steering, change to the rotation frequency, or other control mechanisms to avoid a turbine completely inside the wake of an upstream turbine under surge or pitch with one of the dangerous frequencies. Nevertheless, the effect of a sheared inflow alone has been shown to break the coherence in the azimuthal direction,^{21} which might suggest that the most worrying values of Table III have a low probability of occurring in a real situation. For a sheared inflow, without inflow turbulence, Kleusberg *et al.*^{21} showed that the breakdown of the vortices in the wake occurs first in the positions with lower streamwise velocity, which usually correspond to the region closer to the ground (i.e., with lower *y* values). Hence, even if the disturbance is azimuthally coherent, the flow structures will not be azimuthally coherent if the inflow is not uniform.

In general, the oscillations are reduced with the distance from the turbine, which can be explained by the intrinsic turbulent mechanisms of the wake breaking the structures of concentrated vorticity, as can be seen in Figs. 9 and 10. Without inflow turbulence, the large flow structures are still present at a distance of six diameters from the turbine and cause relevant velocity fluctuations. Hence, the effects of turbine motion can demand higher distances between the turbines. However, before any conclusion is drawn, these phenomena need to be better understood in more realistic inflow conditions, especially with incoming turbulence, and taking into consideration structural properties and the dynamic effects of motion on the operation conditions for the turbine.

The Reynolds number is not believed to have a relevant role in the initial growth of the instabilities and the first non-linear effects, since the main mechanisms of vortex interaction are inviscid, as described in Sec. III. Also, the vortex core size, which is a function of the Reynolds number and of the parameters of the actuator line model (as discussed in Sec. II B, see also Refs. 37 and 38), has only a second-order effect on the stability and dynamics of the vortices in the near wake. However, the Reynolds number and vortex core properties may affect the vortex breakdown and mixing inside the wake, which has an effect on the dynamic pressure, especially for the most downstream position. Hence, there are reasons to believe that the qualitative aspects of Table III are valid for different Reynolds numbers, but more studies are needed to understand its influence.

As expected, the average value of $Fz*$ is lower for the disks closer to the turbine (*z _{c}* = 6) and inside the wake, while the average value of $Ty*$ higher for the disks closer to the turbine and is almost zero inside the wake. Another effect that can be seen for heave in Figs. 11(b) and 11(d) is that $Fz*$ inside the wake is higher and $Ty*$ partially in the wake is lower for $\omega *=0.5$ and $\omega *=2.5$ when compared to the other frequencies. One possible explanation is that the large flow structures provide better mixing with the external flow, bringing streamwise momentum to the wake. This is an indicator that the large flow structures might also have a beneficial effect. The creation of large, coherent flow structures was identified by Yilmaz and Meyers

^{54}as a possible flow control strategy for wind farms. They showed that the creation of coherent vortex rings in the wake of a turbine is an optimum way to increase the power production of a configuration of two in-line wind turbines without inflow turbulence. Turbine motion can be further explored as a possible strategy to create coherent structures that increase power production on wind farms. It should be noted that this strategy is different from the proposal of destabilizing the tip vortices in order to accelerate the wake recovery (e.g., Refs. 55–57), where a more uniform wake is obtained.

### B. Comparison to the stability theory

In Fig. 9, it is possible to note that the interaction of vortices occurs further upstream for frequencies closer to $\omega *=1.5$, while the region of interaction for frequencies $\omega *=0.5$ and $\omega *=2.5$ seems to be comparable. This agrees qualitatively with the stability theory (Sec. III), which predicts the highest growth rate for $\omega *=1.5$ for all types of motion. In order to quantitatively compare the growth rate of perturbations with the stability theory, a lower amplitude of motion of 0.1% is used in the numerical simulations. An amplitude of 1% is too close to the saturation amplitude where non-linear effects start to be relevant, making the linear region too short to allow calculation of the linear growth. The shorter domain is used (*z _{out}* = 11.78), since we are interested in the near-wake dynamics. Qualitatively, the flow behavior is similar to the simulations with an amplitude of 1%, but the transition to turbulence occurs further upstream, as can be seen in Ref. 28.

In order to quantitatively compare the spatial growth rate of the perturbations, the flow is decomposed in its Fourier series in time, as performed in Refs. 19 and 21. During a period $\Delta T$, *N _{S}* = 48 flow-field snapshots $\varphi q,\u2009q=0,\u2026,NS\u22121$, equally spaced in time are used to extract the component $\Phi \u0302$ in each frequency,

where $i=\u22121$ and *n _{f}* is the index (integer) that defines the selected angular frequency, given by $\omega =nf2\pi /\Delta T$. By selecting the frequency of motion for each simulation, the growth of the perturbations can be obtained by analyzing the evolution in space of the maximum of the absolute value of $\Phi \u0302$.

^{19,21}The Fourier mode of the streamwise component of the velocity, $w\u0302$, is used to calculate the growth rate $\sigma =d(log\u2009(max(w\u0302)/W\u221e))/dz$, which is then nondimensionalized and scaled according to $\sigma \u0303=\sigma (2h2wc/\Gamma )$.

^{20,21}Similarly to other studies,

^{18,19,21}after a region of receptivity, the growth is exponential until non-linear effects become relevant and the growth saturates [Fig. 12(a)]. In order to calculate the growth rate, only the region where the growth could be considered exponential was considered.

As can be seen in Fig. 12(b), the growth rate is very close to the values predicted by the stability theory. However, the growth rates are slightly lower, especially for perturbations imposed by heave. The previous works^{21} showed better agreement when using local properties to scale the growth rate of perturbations imposed as body forces to the tip vortices. Since the turbine motion modifies the wake as a whole, more than just perturbing the tip vortices, a possible hypothesis is that these more complex dynamics also influence the stability of the wake. The higher difference for heave seems to agree with this hypothesis, because heave motion also modifies the shape of the wake, while the boundary of the near wake is undisturbed for the surge. More studies are needed to understand this difference. Nevertheless, the good agreement indicates the instabilities predicted by the linear theory are the dominant effect.

Figures 13 and 14 show the real part of the Fourier modes. The perturbations imposed by the heave and surge motions, from Sec. III A 2, are also shown. The perturbations [shown in Figs. 13(b) and 13(e) and 14(b) and 14(e)] were adapted to the left-handed helices, following the reference of Eq. (9) [instead of Eq. (10), used in most of Sec. III].

The similarities between the first-order approximation of the perturbations and the Fourier modes that are present in the wake are clear, which indicates that the modes predicted in Sec. III dominate the flow in the near wake. Also, even for the simulations with an amplitude of motion of 1%, where the linear region is short, the stability theory and the first-order perturbation can explain the flow structures that are seen in the non-linear region of the flow. In Fig. 13, looking at two positions *z* indicated by gray dashed lines, we note that the vortices with positive *y* are displaced in the opposite direction to vortices with negative *y* for heave, while for surge all vortices at a certain position *z* are displaced in the same direction. This explains why the flow structures for heave are tilted. The results for pitch motion seen in Fig. 9(g) show a combination of the effects: the vortices coalesce in one large structure coherent in *z*, but the vortices in the negative *y* position coalesce faster than the vortices with positive *y* since it is formed by the sum of the effects of surge and heave, weighted by the ratio $Rp/R$, as discussed in Sec. III A 3.

In Fig. 14, following one vortex, indicated by a gray dashed line, we note that between $\u2212\pi /3>\theta >\u22122\pi /3$ the heave motion has two lobes while the surge motion only has one lobe. This effect can be noted in the Fourier mode of the simulation with low amplitude [panels (a) and (d)], the perturbation [panels (b) and (e)] and in the displacement of the vortex in the non-linear region of the simulation with higher amplitude [panels (c) and (d)].

Based on the results presented here, we believe that the mechanism that forms the near-wake structures can then be explained as follows. The perturbation on the tip vortices, imposed by the turbine motion, excites vortex instability modes that are practically identical to the perturbation itself. These modes grow exponentially in time and space until non-linear effects become significant. The flow structures that are formed in the non-linear regime are a consequence of this mechanism and preserve some of the characteristics predicted by the linear theory.

Only one type of motion and frequency is considered in each numerical simulation presented here. The linear theory, by assumption, treats different types of motion and frequencies independently. Good agreement of our numerical results with the linear theory suggests that, if there is more than one type of motion or frequency involved, the one that will dominate is that with the higher growth rate according to the linear stability theory. However, further studies are needed to verify if non-linear effects might be relevant for certain cases.

As can be seen in Figs. 9 and 10, the type of motion influences the flow near the center of the turbine. The flow near the hub vortices seems much less disturbed for surge [Figs. 9(e) and 10(e)] than for heave [Figs. 9(a) and 10(a)]. The linear stability of helical vortices might suggest an explanation. The analysis of Sec. III B is not applicable for the hub vortices because the approximations of low pitch are not valid. The radial position of the hub vortices is lower than the distance between neighboring vortices. In this case, a more elaborate model is needed. First, we note that the approximations that $\delta z\u2248\delta r$ and the azimuthal component is not relevant are not valid. For higher pitch, a perturbation in the streamwise direction has lower importance (in the limit of infinite pitch, the vortex is indifferent to a perturbation in the streamwise direction). Hence, the perturbations imposed by heave are more similar to the eigenvectors of the stability theory than perturbations imposed by the surge. In addition to that, Quaranta *et al.*^{18} showed, by transposing the results of Robinson and Saffman^{58} to the helical system, that the growth rate of the modes in each “parabola” of the growth diagram decreases as the wavenumber increases, an effect that is more pronounced for higher pitch (this effect that can also be seen in Fig. 5, to a lesser extent). The “parabola” with the highest growth is that with the peak closer to *k *=* *0, which corresponds to modes with a phase difference [Fig. 6(b)], that are excited by the heave motion. Thus, even if the modes were directly excited, the mode excited by surge would have a lower growth rate. Therefore, both effects indicate that perturbations in the hub vortices by heave motion will reach higher amplitude earlier. In real turbines, the nacelle and its wake will greatly influence the flow in this region, so it is uncertain if such an effect would occur in more realistic configurations. Further analyses of the effect of motion on the hub vortices are left for future studies.

All the analyses performed in this work considered uniform inflow without turbulence. In the case of sheared inflow or yawed turbine, tip vortex instabilities from the linear stability theory were already shown to dominate the near-wake stability when vortices are disturbed by body forces.^{21,22} Hence, we believe that the mechanism described in this study for the formation of the near wake of turbines under motion would also be applicable to shear inflow conditions or yawed turbines, with the necessary adaptations to the theory of Sec. III.

However, the modes described in this work might not dominate the flow for a FOWT in a turbulent atmospheric boundary layer. The turbulence level found in real-world applications might be in the same order of magnitude or higher than the perturbations imposed by the motion. Thus, free-stream turbulence could dominate the near wake, preventing the growth of the perturbations imposed by turbine motion. The large flow structures created by turbine motion with $\omega *=0.5$ only occur due to the phase shift between different helices being restricted because the turbine moves as one. Inflow turbulence does not have this restriction, and, therefore, several frequencies could have a growth rate $\sigma \u0303\u2248\pi /2$, higher than the growth of disturbances with $\omega *=0.5$, as demonstrated by Sarmast *et al.*^{20} Alternatively, turbulence could dissipate or modify flow structures created in the near wake in such a way that they might not be relevant for a downstream turbine. The high turbulence level has been shown to accelerate the disintegration of coherent vortex rings in the simulations of Yilmaz and Meyers^{54} Also, for a fixed-bottom wind turbine with turbulent inflow, De Cillis *et al.*^{59} observed that the modes due to tip vortex instabilities are superseded by modes related to flow fluctuations coming from incoming turbulence. Nevertheless, the motion of floating turbines imposes very specific perturbations, that are not present on fixed-bottom turbines. Hence, further studies are needed to investigate the effects of turbulence on moving turbines.

## V. CONCLUSIONS

The wake of a floating offshore wind turbine was analyzed by means of numerical simulations and comparison with linear stability theory. For the numerical simulations, the use of actuator lines with a high-order SEM code made the analysis of the vortex interactions possible. Two simplified models of vortex stability that are applicable to all types of motion were developed: a model based on the stability of a row of vortices and a model based on the stability of helical vortices.

Using a first-order approximation, it was noted that the perturbations induced by the heave and sway motion have wavenumbers of $k=\omega *\u22121$ and $k=\omega *+1$, in contrast to the perturbations imposed by surge forces, which have wavenumber $k=\omega *$. Modes compatible with these wavenumbers were observed in the numerical simulations in the near-wake and formed flow structures that are observed in the non-linear region of the wake.

Despite the different wavenumbers created by different types of motion, the growth rate depends mostly on the frequency for helices of low pitch and is compatible with a model of an infinite row of two-dimensional vortices calculated by Lamb.^{44} The stability theory of a row of two-dimensional vortices also predicted the number of vortices that interact with each other in the numerical simulations. For $\omega *=1.5$, our simulations showed the out-of-phase mechanism of two vortex pairing observed in other numerical^{19–21} and experimental works.^{13,17,18} For $\omega *=0.5$ and $\omega *=2.5$, larger structures with the interaction of six vortices were observed for all types of motion simulated. These larger flow structures created higher oscillations of the streamwise velocity inside the wake. The RMS of the integrated dynamic pressure was considerably higher for this case, especially for surge and pitch motions, which create coherent structures in the azimuthal direction. This could increase the fatigue loads and the amplitude of the motion of downstream FOWTs.

Our simulations, without inflow turbulence, showed that turbine motion excites the unstable modes of helical vortices. These modes grow in time and space, forming unsteady flow structures that may affect downstream turbines. Free-stream turbulence, not considered here, could dissipate or modify these flow structures. Hence, future studies are necessary to confirm if these effects are sustained under more realistic conditions.

The results suggest that the linear theory can be used to predict which modes imposed by turbine motion will dominate the flow if more than one type of motion or frequency is present. However, further studies are needed to understand the impact of non-linear interactions. This study shows that the general qualitative behavior of the wake can be understood and predicted from relatively simple stability models. A good understanding of the wake effects of turbines under motion is necessary to develop distancing parameters for wind farms of FOWTs. Knowledge about which frequencies and types of motion influence downstream turbines can also guide control laws that avoid dangerous situations or promote favorable flow conditions.

## ACKNOWLEDGMENTS

The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the PDC Center for High-Performance Computing at the Royal Institute of Technology (KTH) and the National Supercomputer Centre at Linköping University. This work was conducted within StandUp for Wind. V.G.K. thanks KTH Engineering Mechanics for partially funding this work. B.S.C. acknowledges the support from FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo), Project No. 2019/01507–8, for this research. B.S.C. also acknowledges the Brazilian National Council for Scientific and Technological Development (CNPq) for financial support in the form of a productivity grant, No. 314221/2021-2.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**V. G. Kleine:** Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Visualization (lead); Writing - original draft (lead); Writing - review and editing (lead). **L. Franceschini:** Conceptualization (equal); Investigation (supporting); Methodology (supporting); Writing - original draft (supporting); Writing - review and editing (supporting). **B. S. Carmo:** Conceptualization (equal); Funding acquisition (lead); Investigation (supporting); Methodology (supporting); Supervision (equal); Writing - original draft (supporting); Writing - review and editing (supporting). **A. Hanifi:** Conceptualization (equal); Funding acquisition (supporting); Investigation (supporting); Methodology (supporting); Project administration (supporting); Resources (supporting); Supervision (equal); Writing - review and editing (supporting). **D. S. Henningson:** Conceptualization (equal); Funding acquisition (lead); Investigation (supporting); Methodology (supporting); Project administration (lead); Resources (lead); Supervision (equal); Writing - review and editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: FIRST-ORDER APPROXIMATION OF PERTURBATION IMPOSED BY YAW AND ROLL

To derive Eq. (4), it is assumed that the center of rotation for pitch is located on the *y*-axis of the turbine. The same assumption is applied here for roll. The center of rotation for yaw is assumed on the *x*-axis of the turbine. For most configurations of FOWT, these are reasonable assumptions. Nevertheless, the formulas presented here would not be applicable for unusual configurations in which the distance in the *z*-direction between the center of rotation and the turbine is in the same order of magnitude as *R* or the distance in other directions.

Yaw motion (subscript _{yw}) is conceptually similar to pitch, hence it also generates similar components

where the position of the center of rotation, *R _{yw}*, would usually be zero, since the axis of rotation is normally the tower of the turbine. For unusual configurations, such as platforms with multiple turbines,

*R*, can be different from zero.

_{yw}Analogously, for roll motion (subscript _{r}),

which becomes

where the component with amplitude $RAr$ affects only the azimuthal perturbation while its radial component is canceled out. Hence, for the purpose of our simplified analysis that neglects the azimuthal component, the roll motion is mainly similar to a sway motion of amplitude $RrAr$.