The Stability of Wakes of Floating Wind Turbines

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, influencing the flow reaching downstream wind turbines. In this work, the wake of a FOWT is analysed by means of numerical simulations and 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 high amplitude of oscillations in the streamwise velocities, that may increase fatigue or induce high amplitude motion on downstream turbines.

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, influencing the flow reaching downstream wind turbines. In this work, the wake of a FOWT is analysed by means of numerical simulations and 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 high amplitude of oscillations in the streamwise velocities, that 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 the past few 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 meters), the bottommounted configuration is no longer feasible and Floating Offshore Wind Turbines (FOWT) 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 on the structure of the wake of the turbine, possibly having important performance impacts in 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, Sant, and Micallef 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 modelled by Lagrangianbased 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][16][17][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 to the blade tips, and the goal was to investigate the evolution and growth of instabilities. However, since the focus was not 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 on 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 analyse 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 section II, the numerical method used for the simulations is described. A first-order approximation of the stability properties of the flow is presented in section III. The numerical results are shown and compared to the stability theory in section IV. Finally, the main conclusions are presented in section V.

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 ∞ .
The turbine is positioned at the origin of a cylindrical domain of radius R rad = 5 ( figure 1(a)).
The distances of the inlet and the outlet to the turbine are z in = 5 and z out = 22.10, respectively, for the longer domain used in the simulations reported in section IV A. Around the center of the domain, the elementwise discretization in the streamwise direction is uniform between −0.6 < z < 14.025, with constant spacing ∆z = 0.075. This discretization is consistent with previous studies 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 out = 11.78 and uniform streamwise discretization between −0.6 < z < 8.025 is used for comparison of the growth rate and modes with the stability theory.
As can be seen in figures 1(b) and (c), the discretization is coarser in regions farther from the turbine, in order to reduce computational cost. Dirichlet boundary conditions with constant velocity W ∞ = 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 ∆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 zcomponent 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/R rad ) 2 = 4%. The results were compared to a simulation with a domain with R rad = 10 (the domain and grid described in Ref. 23), with a blockage ratio of 1%. Differences in 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 R rad = 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.

B. Wind Turbine Modelling
The wind turbines are modelled 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 N ACL points (x n ) along the blades. The discrete two-dimensional force vector (F 2D ) 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 , ε ≈ 3.5∆r (where ∆r is the averaged grid spacing) and each blade is discretized with 100 points. The turbine is modelled after the turbine of the Blind Test 31 . The chosen operation condition was the tip-speed ratio λ = ΩR/W ∞ = 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][22][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 Re R = W ∞ R/ν = 50, 000 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 in the order of Re tip c = W ∞ λ c t /ν ≈ 10 5 , 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][20][21][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 figure 2.
In Ref. 28 we considered heave motion only. To assess the amplitude and frequency of heave, Simulations were run considering water depths from 500 m 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 a 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 section 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 (x t , y t , z t ) of the actuator lines: where ω is the angular frequency of motion and the subscripts y and z represent heave and surge, respectively. We define ω * = ω/Ω 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 to 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 ( where φ n = 2πn/N b is the initial azimuthal position of each blade n (n = 0 to N b − 1, where N b = 3 is the number of blades) and A p 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 Sway (translation in the x-direction), yaw (rotation around the y-axis) and roll (rotation around the z-axis) are not simulated. However, the first-order approximation of the perturbations and discussion on the stability of these motions are presented in section III and appendix A.

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

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 (x t0 , y t0 , z t0 ) = (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 (x v0 , y v0 , z v0 ) of the tip vortices emitted by the blades of equation (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 advection of the vortices. The position (x v , y v , z v ) 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 by 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 analysed.
The helical structure of the vortices of equation (6) 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 τ 0 = 0, equation (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 equation (9)  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 equation (10) are the same of equation (9), but the comparison with classical results can be performed without further modification to the reference system by using equation (10). On the other hand, comparison to the simulations should be made considering the different reference systems.

First-order approximation of perturbation imposed by sway, heave and surge
Using the reference system of equation (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 (δ A y sin (ω y τ), A z sin (ω z τ)), in cartesian coordinates. For A x , A y R, first-order approximations of the perturbation in the radial direction and azimuthal directions are The order of magnitude of δ r v and Rδ ϕ 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 1 ⇐⇒ w c /(N b ΩR) 1), so, as an approximation, it is neglected.
Defining θ = Ωτ + φ n = Ωz v0 /w c + φ n , the first-order approximation of the perturbation in the radial and streamwise directions are defined for the helices in function of the azimuthal angle as  This initial perturbation induced by the motion grows in time and in 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 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 equation (14) (pure sinusoidal oscillation for surge and sinusoidal with sin θ and cos θ envelopes for heave and sway), then the analysis performed here would also be applicable.

First-order approximation of perturbation imposed by pitch
From equation (4), applying the change in reference system discussed in section 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 R p R this corresponds to a surge motion with amplitude R p A p . For R p R this corresponds to a motion in z-direction composed of component sin (ω * p (θ − φ n )) sin θ with amplitude RA p , which has the same sinusoidal envelope of the heave component in equation (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 section 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 smallperturbation 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 (ω * (θ − φ n )) sin θ ), (sin (ω * (θ − φ n )) cos θ ) or (sin (ω * (θ − φ n ))). Conceptually, for the purposes of the analysis of the following sections, these components are equivalent to the perturbations of equation (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 2-d Row of Vortices
Several studies made the connection between the stability of helical vortices and a twodimensional 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 figure 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, σ = σ (2h 2 /Γ), has a maximum value of π/2. The growth rate of the vortex pairing mode collapse to π/2 even for non-uniform helices created by turbines under sheared inflow or yawed, if scaled by local properties, as shown by Kleusberg, Benard, and Henningson 21 and Kleusberg 22 . The two modes (δ r m 2d 1 , δ z m 2d 1 ) and (δ r m 2d 2 , δ z m 2d 2 ) predicted by the linear stability theory for each p where i is the complex unit. These modes are illustrated in figure 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 figure 3, which corresponds to the reference system used in this work for a wind turbine, mode (δ r m 2d 1 , δ z m 2d 1 ) 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 θ = θ 0 + 2π j, where j is an integer and θ 0 is chosen to define a cross-section that contains the row of vortices (0 ≤ θ 0 < 2π). This is a slightly different approach from that in Quaranta, Bolnot, and Leweke 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 equation (14), the displacement of the vortices due to surge, heave and sway along the cross-section defined by θ 0 is where the term N b j − n, that can be any integer, has been reinterpreted as the index m of each 2-d vortex ( figure 3).
Projecting the terms of equation (20) onto the modes defined in equation (19), the perturbation 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 (δ r v 2d 1 , δ z v 2d 1 ) tends to dominate. It can be noted that there is a symmetry between the perturbation in the component δ r v 2d and δ z v 2d : a perturbation in δ r v 2d induces a component in δ z v 2d and vice-versa, because the eigenvectors have equal components in r and z.
Qualitatively, the three terms of equation (21)  From the terms of equation (21), it follows directly that the terms iω * θ 0 correspond to a phase shift, and, most importantly, ω * /N b is the dimensionless subharmonic wavenumber p for each motion type. Applying the scaling p = ω * /N b to equation (18), we obtain the theoretical growth rate for a row of vortices shown in figure 4. Equation (18), in principle, is only valid for 0 ≤ p ≤ 1.
For values of ω * /N b that are outside this interval, however, we use p = ω * /N b + j where j is any integer that makes 0 ≤ p ≤ 1 (including for negative values of ω * ).
As a linear approximation, this growth rate is applicable to all motion types, including pitch, yaw and roll, following the discussion of section 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, Benard, and Henningson 21 (14) is composed by two terms and   δ r v2 where (δ r v1 , δ z v1 ) grows and (δ r v2 , δ z v2 ) decays in time and space (for the current choice of reference system).
Using trigonometric identities, equation (23) can be decomposed in and analogously for equation (24). The wavenumber k of each perturbation is given by the term that multiplies θ , while the terms ±φ n are phase shifts between the different vortices n. Thus, a sway or heave motion (x and y-components) with frequency ω * creates perturbations with wavenumbers k = ω * − 1 and k = ω * + 1. The surge motion (z-component), on the other hand, creates perturbations with wavenumber k = ω * . 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 = ω * 21 .
Analogously to the row of vortices, when the real perturbation is projected onto the eigenvectors, a pair of complex conjugate modes is excited, see equations (21) and (22). Therefore, the heave and sway motions excite modes with wavenumbers k = ±(ω * − 1) and k = ±(ω * + 1) and the surge motion excite modes with wavenumbers k = ±ω * . In figure 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 ω * excite eight eigenmodes: • Two unstable eigenmodes with k = ω * − 1 and k = ω * + 1 and their complex conjugates (k = −(ω * − 1) and k = −(ω * + 1)); • Two stable eigenmodes with k = ω * − 1 and k = ω * + 1 and their complex conjugates (k = −(ω * − 1) and k = −(ω * + 1)); while the surge motion excites four eigenmodes: stable and unstable modes with k = ±ω * . As can be seen in equation (25), the sway motion excites the same modes of the heave motion, only with a phase difference of π/2. The example of heave and sway motions for ω * = 1.5 is shown as an example in figure 5. The modes with k = ω * − 1 and k = ω * + 1 are referred as the primary modes and the modes k = −(ω * − 1) and k = −(ω * + 1) as their complex conjugate.  (25), it is expected that each primary mode excited by one frequency belongs to a different branch, because one term has phase difference +φ n and the other −φ n . This is exactly what is observed: the mode with k = ω * − 1 belongs to one branch (related to −φ n ) and the mode with k = ω * +1 belongs to the other branch (related to +φ n ). It is easy to see that the complex conjugate of the primary mode belongs to the opposite branch.
As can be seen in figures 7 and 8, the perturbations created by the heave motion (equation (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, a 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 comparing to the modes of figure 6(a), when considering the modes with k = ±ω * . 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 • There are two eigenvalues for k = ω * − 1 = −0.5, one close toσ /(π/2) ≈ 0.55 and another close toσ /(π/2) ≈ 1.0. The eigenvector that is observed to be identical to the perturbation (see figure 7) is that whose eigenvalue has growth rateσ /(π/2) ≈ 0.55, because its position in relation to the parabola between −1 < k < 2 is the same of ω * = 0. • There are two eigenvalues for k = ω * + 1 = 1.5, both close toσ /(π/2) ≈ 0.55 (it is close to the intersection of parabolas between −1 < k < 2 and between 1 < k < 4). The eigenvector that is observed to be identical to the perturbation (see figure 7) is that whose eigenvalue belongs to parabola between 1 < k < 4, because the position of the eigenvalue in relation to this parabola is the same of ω * = 0.5 in figure 4.
The term "parabola" was freely used here, but we note that the curves do not represent exact parabolas, except for figure 4, that comes from equation (18). The discussion presented here provides an explanation for this rule-of-thumb: the branch that corresponds to −φ n is similar to the branch without phase difference ( figure 6(a)), but displaced one unit to the negative direction in the k-axis (k = ω * − 1), while the branch that corresponds to +φ n is similar to the branch without phase difference, but displaced one unit to the positive direction in the k-axis (k = ω * + 1).
Due to this connection between the excited modes and the stability of 2-d 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 ω * = 1.5 corresponds to the vortex pairing mechanism that has the highest growth rateσ ≈ π/2, which has been discussed in previous studies that excite the vortices using body forces [19][20][21] . Additionally, we note that all motions with an angular frequency ω * = 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 = N b ( 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σ ≈ π/2 are given by ω * = N b ( 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 = ω * (e.g. surge) or the phase shift is ±2π/3 and k = ω * ± 1 (e.g. heave). This means, for example, that a frequency of ω * = 0.5 would never excite the highest growth modeσ ≈ π/2 for k = 0.5, seen in figure 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 θ and cos θ in equation (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, equation (A3), would also create wavenumbers given by k = ±(ω * ± 1). Pitch, equation (17), and yaw, equation (A1), have components with and without terms sin θ and cos θ . Therefore these motions would induce perturbations composed of several wavenumbers: the terms RA p and RA yw induce wavenumbers k = ±(ω * − 1) and k = ±(ω * + 1) and the terms R p A p and R yw A yw induce wavenumber k = ±ω * .
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][51][52][53] are not captured by the model.

A. Direct numerical simulations
Several cases with different types of motion and parameters were simulated, according to table I. The reference case has amplitude A/R = 1% and angular frequency ω = 0.5Ω, according to estimates detailed in section 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/R p = 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][20][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 section IV B, in order to compare the growth rate to analytical stability predictions.  Interestingly, large flow structures, very similar to those observed for ω * = 0.5, can also be noted for the higher frequency, ω * = 2.5. According to the stability theory for a row of vortices, this is expected. Every frequency that can be written as ω * − jN b , with j integer, is indistinguishable from ω * according to the theory presented in section III B. It follows that the configuration of a row of vortices for ω * = 2.5 is equal to ω * = −0.5, which is the mode complex conjugate of ω * = 0.5. Thus, the number of vortices that interact directly is predicted to be the same for 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 ω * = 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 6h, as predicted. The shape of the large flow structures is influenced by the type of motion. For ω * = 0.5,  showing an effect of the (sin θ ) 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 ω * = 0.5 is the same frequency 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 F * z and T * y were defined based on integrals of the streamwise dynamic pressure. F * z is an indicator of the streamwise force that a fictitious disk, with the same radius of the turbine (R = 1) centered at (x c , y c , z c ), 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. T * y 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 − x c < 0 is zero and the velocity on the semicircle of x − x c > 0 (area S/2 and centroid at 4R/(3π)) is constant and equal to W ∞ . and (x c , y c , z c ) = (1, 0, 12). The values of the average in time of F * z and T * y are shown in figure 11, with error bars indicating their root mean square (RMS). In   As can be observed in table III and figure 11, in general, the RMS of the integrals of the dynamic pressure are higher for angular frequencies ω * = 0.5 and ω * = 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 in T * y is observed for a disk partially in the wake ((x c , y c ) = (1, 0)).
Surge and pitch cause great oscillation in F * z 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 low probability of occurring in a real situation. For a sheared inflow, without inflow turbulence, Kleusberg, Benard, and Henningson 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 figures 9 and 10. Without inflow turbulence, the large flow structures are still present at a distance of 6 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 on the initial growth of the instabilities and on the first non-linear effects, since the main mechanisms of vortex interaction are inviscid, as described in section 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 section 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 F * z is lower for the disks closer to the turbine (z c = 6) and inside the wake, while the average value of T * y higher for the disks closer to the turbine and is almost zero inside the wake. Another effect that can be seen for heave in figures 11(b, d) is that F * z inside the wake is higher and T * y partially in the wake is lower for ω * = 0.5 and ω * = 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 Yılmaz and Meyers 54 as a possible flow control strategy of 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 such 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][56][57], where a more uniform wake is obtained.

B. Comparison to the stability theory
In figure 9, it is possible to note that the interaction of vortices occurs further upstream for frequencies closer to ω * = 1.5, while the region of interaction for frequencies ω * = 0.5 and ω * = 2.5 seems to be comparable. This agrees qualitatively with the stability theory (section III), which predicts the highest growth rate for ω * = 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 de- where i = √ −1 and n f is the index (integer) that defines the selected angular frequency, given by ω = n f 2π/∆T . By selecting the frequency of motion for each simulation, the growth of the perturbations can be obtained by analysing the evolution in space of the maximum of the absolute value of Φ 19,21 . The Fourier mode of the streamwise component of the velocity, w, is used to calculate the growth rate σ = d(log(max( w)/W ∞ ))/dz, which is then nondimensionalized and scaled according toσ = σ (2h 2 w c /Γ) 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 ( figure 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 figure 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. 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 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. 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 R p /R, as discussed in section III A 3.
In figure 14, following one vortex, indicated by a gray dashed line, we note that between −π/3 > θ > −2π/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 (a,d), the perturbation (b,e) and in the displacement of the vortex in the non-linear region of the simulation with higher amplitude (c,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 which are practically identical to the perturbation itself. These modes grow exponentially in time and space until non-linear effects become sig- nificant. The flow structures that are formed in the non-linear regime are a consequence of this mechanism and preserves 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. The good agreement of our numerical results with the linear theory suggests that, if there is more than one type of motion or frequencies 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 figures 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 (figure (e)) than for heave ( figure (a)). The linear stability of helical vortices might suggest an explanation. The analysis of section 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 δ z ≈ δ 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 surge. Besides 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 figure 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 ( figure 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 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 sheared inflow conditions or yawed turbines, with the necessary adaptations to the theory of section 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 ω * = 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 growth rateσ ≈ π/2, higher than the growth of disturbances with ω * = 0. 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 behaviour 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.
Appendix A: First-order approximation of perturbation imposed by yaw and roll To derive equation (4), it is assumed that the center of rotation for pitch is located on the yaxis 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 z-direction between the center of rotation and the turbine is in the same order of magnitude of R or the distance in other directions.
Yaw motion (subscript yw ) is conceptually similar to pitch, hence it also generates similar com- 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 yw can be different from zero.
Analogously, for roll motion (subscript r ), where the component with amplitude RA r affects only the azimuthal perturbation while its radial component is cancelled 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 R r A r .