Aerodynamic Analysis of a Two-Bladed Vertical-Axis Wind Turbine Using a Coupled Unsteady RANS and Actuator Line Model

.


Introduction
Climate change mitigation is vitally important for all nations in the world, given that greenhouse-gas (GHG) emissions have increased by over one-quarter since 1995 [1], as reported at the first United Nations (UN) Conference of the Parties (COP).Moreover, energy consumption by developed and developing countries has been projected to increase by 28% from 2015 to 2040 [2].A key approach to replacing fossil fuels as an energy source and limiting carbon release is to invest in renewable energy technology [3].Wind and hydrokinetic energy are particularly attractive options for sustainable electricity generation from low-carbon sources [4], and are likely to become significant contributors to the electricity supply by 2030 [1].Much ongoing research into the development of wind and tidal turbines focuses on horizontal-and vertical-axis turbines [5].Salter [6] compared vertical-axis transverse-flow turbines with horizontal-axis axial-flow turbines in terms of flow impedance, turbulence, blockage ratio, installation, pitch change, and navigation, with tidal flow in the Pentland Firth, Scotland, in mind.Salter found that high blockage (or sweepage), vertical-axis, variable-pitch rotors could lead to substantially higher potential power generation for high impedance flows [6].Such vertical-axis transverse-flow tidal turbines tolerate uneven seabed topography and may attain an even pressure drop by controlling the blade pitch, hence reducing wake turbulence [6].Vertical-axis turbines thus appear to offer a promising near-term technology for tidal energy.Initial study of vertical-axis turbine (VAT) technology began in the 1970s at Sandia National Laboratories where researchers investigated vertical-axis turbine configurations, including Savonius (torque generated from drag) and Darrieus (torque generated from lift) turbines [7,8].The Savonius turbine can accept flow from any direction and is self-starting, with low cut-in speed; however, the Savonius turbine is restricted to fewer applications due to its inefficiency at relatively low tip speed ratios [9].Darrieus turbines have higher cut-in speed than equivalent Savonius turbines, and so rotate faster than the inflow velocity, attaining higher coefficients of performance [9,10], even though their support arms introduce additional aerodynamic drag [11].To solve this problem, Salter and Taylor [12] proposed the innovative vertical-axis rotor system shown in Figure 1.Computational fluid dynamics (CFD) has been widely used in the systematic analyses of vertical-axis turbines [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29].Actuator-type models parameterize the turbine loading and thus reduce computational expense, but do not resolve the fine detail of the blade boundary layers [30].Four approaches have commonly been used to represent turbines in such models, namely: actuator disc with rotation or blade element momentum (BEM) [31][32][33][34][35]; actuator disk without rotation [30,35,36]; actuator surface [37][38][39]; and actuator line [30,40,41].BEM is an analytical method, whereas the actuator disc with rotation model is a combination of blade-element (BE) theory and CFD, which solves the Navier-Stokes equations to satisfy the momentum balance [35].The actuator disc with rotation model is computationally efficient, but does not directly include the influence of vortices shed from blade tips on the induced velocity [31].The uniform actuator disk without rotation model is limited in applicability because of its simplifying assumptions [37], and has proved unsatisfactory as a wake generator method for a crossflow turbine [40].The actuator surface technique accurately predicts the flow structure near blades and in the tip vortex region, but requires a fine mesh passing smoothly over the airfoil surface [38].
Even so, LES is substantially more expensive computationally than RANS, which is why it is used rather sparingly in simulations of turbulent flow past horizontal-axis turbines and vertical-axis turbines.
Typical recent applications of CFD to turbines follow.McLaren [57] reported a numerical and experimental study of the unsteady loading on a small-scale, high-solidity, H-type Darrieus turbine, based on two-dimensional (2D), unsteady Reynolds-averaged Navier-Stokes (URANS) simulations by CFD ANSYS-CFX.The study revealed the dominant effect of dynamic stall on the output power and vibration excitation of the turbine.Nobile et al. [58] later simulated 2D unsteady-flow past a Giromill wind turbine, also using ANSYS-CFX, finding that mesh resolution and choice of turbulence model had a substantial effect on accuracy, with time step having only a slight impact on the numerical results.Biadgo et al. [59] used a stream-tube approach to undertake a numerical and analytical assessment of the performance of a vertical-axis wind turbine comprising a straight-bladed fixed-pitch Darrieus turbine with a NACA 0012 blade profile using ANSYS FLUENT.These numerical predictions were compared with analytical results obtained using a double multiple streamtube (DMST) model, which exhibited inability using both CFD and DMST for the turbine to be self-starting owing to minimum and/or negative torque and performance at very low tip-speed ratios.Bachant et al. [60] developed a validated ALM of a vertical-axis turbine with both high and medium values of solidity, and tested both - RANS and Smagorinsky LES turbulence models in the OpenFOAM CFD framework.Bachant et al. found that RANS models running on coarse grids were able to provide good convergence behaviour in terms of the mean power coefficient.Compared with other 3D blade-resolved RANS simulations [60,61], Bachant et al.'s model achieved approximately four orders of magnitude reduction in computational expense by implementing corrections in sub-models for the effects of dynamic stall, end conditions, added mass, and flow curvature.Given that such models have focused on idealized vertical-axis turbines, further investigation into optimal practical models with fewer correction factors is still required.
Figure 1 shows a group of close-packed contra-rotating vertical-axis rotors, designed by Stephen Salter to maximise the fraction of flow passage swept [12].Blockage is estimated to increase to 80% given the small gaps between the rotors, which are controlled by a hydraulic ram.The rotor diameter should be at least three times the water depth in order to provide stability in pitch and roll of a single rotor, and this should be doubled for a close-packed array.This contributes to a high blockage fraction allowing generation well above the Betz limit for rotors in channels [6].Following Buntine and Pullin [62], the design concept is based on two vortices of opposite-sign cancelling each other out, and thus conditioning the flow though the turbine while lowering the turbulence kinetic energy in the wake.The turbine downstream area will then experience less stream-wise flow variation, reducing mixing loss and therefore enhancing energy extraction.To predict the commercial feasibility of this large-scale marine hydrokinetic application, a numerical model of such devices is required.This paper describes a numerical model of a cross-flow turbine, with the future goal of modelling close-packed tidal rotors comprising many blades.The present model is built upon a previous turbine model, which scales to thousands of cores on a supercomputer [54,56].Although the present focus is on a single rotor, the numerical model can be applied to a large-scale turbine farm in future studies.
Due to a lack of experimental data concerning this type of rotor, the numerical model is first validated against experimental measurements from a two-bladed H-type wind turbine, and then used to predict turbine loading and investigate vorticity distribution in the vicinity of the rotor.
A newly developed, efficient, parallelised, numerical model of vertical-axis turbines, with a fixed tip-speed ratio system and with a torque-controlled system, is presented in the following sections.
This computationally efficient numerical model is coupled with and is developed within the OpenFOAM CFD framework.Unique features of the present model include torque control and active pitch mechanisms.For brevity, only the torque-controlled system is presented in this paper; pitch control mechanisms for solving the dynamic stall problem as well as performance optimization [63,64] will be explored in future work.We believe that the application of the present model to a torquecontrolled vertical-axis turbine gives new insight into the aerodynamic behaviour of vertical-axis wind turbines, in particular the difference in behaviour between an idealised turbine with fixed tipspeed ratio and a more practical turbine with torque control.

Mathematical Model
Flow past a single vertical-axis turbine (VAT) with an arbitrary number of blades is simulated using an adapted version of the Wind and Tidal Turbine Embedded Simulator (WATTES), which is an efficient, parallelised, two-way coupled turbine model of horizontal-axis turbines, scaling to thousands of computing cores [54,56].We denote the newly developed model WATTES-V.A preparatory set-up of the original WATTES model using the OpenFOAM CFD solver was conducted to ensure the codes were correctly coupled [65]; details of the software architecture are provided in Appendix A. This prerequisite ensures that WATTES-V model benefits from the advantages of the original model.One unique feature of the modified WATTES-V model is that it enables torque control; the main benefit of torque-controlled models is their prediction of the dynamic response of the turbine to the flow [52][53][54].The mathematical formulation of WATTES-V is provided below.

Frame of Reference
To calculate the body forces, the coordinates of nodes in the mesh are first translated to the frame of reference of the rotor, in a similar manner to the original WATTES model [54].The centre of the vertical-axis turbine is located at position  (see Figure 2), where  5 6666⃗ = ( 5 ,  5 ,  5 ).The azimuthal angle, which describes the orbital path taken by the first turbine blade, is denoted .In WATTES-V,  starts from the -axis, as indicated in Figure 2. The coordinates of a blade reference frame are denoted  ?,  ?,  ?, with  ?( ?,  ?,  ? ) the origin of the new reference system.In the blade reference frame, the coordinates of a transformed point at position  ?666⃗ = ( ?,  ?,  ? ) are: where Similarly, the localised velocity at a given point is  6⃗ = (, , ), and this is transformed to the rotor's frame of reference as  ?666⃗ = () 6⃗.Once in this frame of reference, the model calculates the momentum source terms, and then a second transformation takes place before passing these back to the CFD solver (cf.Creech et al [54]).To simplify the notation, we denote the transformed coordinates and velocity as  ⃗ and  6⃗ hereafter.

Lift and Drag Calculations
The actuator line method (ALM) [43] creates a distribution of body forces along a set of line segments representing the blades of a turbine.For each turbine rotor, only grid points found within the hollow cylindrical volume  traced out by the rotating blades are considered.
The lift and drag force vectors per unit span on a blade are given by: where  is the fluid density, and  M and  ( are the lift and drag coefficients, which depend on the angle of attack  and the Reynolds number  of the flow over the blade.The magnitude of relative velocity of the fluid over the blade is | ]^N 6666666⃗|, and () is the blade chord length, which can vary along the blade span, but in the present case is constant.As the blades are parallel to z-axis, this is a function of .The unit vectors  M 666⃗ and  ( 6666⃗ are in the direction of lift and drag respectively.Values of  M and  ( are given in tabulated form [54], and as with most models, these are derived from an assumption of two-dimensional flow over the blade.Figure 2   The relative velocity  ]^N 6666666⃗ is calculated for each point within the control volume  at a radial distance  from the rotor center (along -axis) as where  cN 66666⃗ is the blade velocity.For a vertical-axis turbine, the magnitude of  ]^N 6666666⃗ is where  cN =  cN , and  cN is the angular velocity of blade.Note that the spanwise velocity component is neglected here, because the spanwise component of flow velocity is assumed to have minimal impact on the performance of the blade, and so tip-loss effects can be ignored.The azimuthal component of the fluid velocity is given as This is necessary to account for the rotation of the flow, as lift and drag forces act to turn the blades and the generator, resulting in an equal and opposite reaction force acting on the flow, causing it to rotate in the opposite direction to that of the blades [54].
The flow angle relative to that of the fluid is The local angle of attack is then computed from  ]^N as follows: where the local blade angle  is given by The blade pitch angle  U can be actively controlled, as with [54], but for the present validation work it is kept constant at  U = 0.The local blade twist angle  Q is calculated from the blade geometry but we consider straight blades and hence  Q = 0 in the present test cases.
Lift and drag forces per unit span are then calculated using the WATTES-V actuator line representation of each blade, which utilises a two-dimensional Gaussian regularization kernel  O ( O ) [56]: where  cN is the number of blades,  O is the shortest distance between a given point and the  th actuator line.The pointwise lift and drag per unit span,  M 666⃗ O and  ( 666⃗ O , are obtained from Equations ( 3) and ( 4).A two-dimensional Gaussian regularization kernel operates in the blade azimuthal direction and smears the solution in a circle [56], such that: where the distance from the  th vertical actuator line is The value of  was chosen carefully so that it is neither too large (smeared solution) nor too small (extremely high resolution, and correspondingly small time step) [56].Experiments determined that numerical stability was optimal when the Gaussian width was set to twice the local cell length, ∆ , as also by Troldborg [30,43].Other researchers have investigated the effect of the standard deviation (or projection width) on accuracy and stability: Schito and Zasso [30,66] found that the equivalent of the mesh cell width was ideal; Jha et al. [30,67] recommended using an equivalent elliptic planform for its calculation; Martinez-Tossas and Meneveau [30,68] used two-dimensional potential flow analysis to determine the optimal projection width; Tennekes and Lumley [69] recommended the projection width to be of the order of the momentum thickness  †Q [30].Here, the Gaussian width related to mesh size is estimated as ∆ ≈ g n^NN ˆ where  n^NN is the cell volume.
Following Bachant et al. [30], an additional factor  †^T‰ = 2.0 is introduced, and non-unity aspect ratio cells incorporated using  = 2 †^T‰ ∆.This meant that 95.45% ( O ≤ 2) of the Gaussian distribution was captured within the numerical simulation.It should be noted that  is a tuning factor that should be adjusted to the particular circumstances under consideration.
The tangential  Q and normal  S components of body forces acting on the fluid, which are in the opposite directions to the force acting on the blade, are given by Body force components acting on the fluid in  and -axis directions are where  • is also the net thrust component of the fluid to the turbine.Note that  i = 0, as threedimensional flow effects on performance are neglected.
All the calculated force terms are then transformed into body force components, and passed back to OpenFOAM as momentum sources in the Navier-Stokes momentum equation for an incompressible Newtonian fluid given by: in which  6⃗ is velocity field vector,  is fluid density,  is pressure,  is the kinematic viscosity,  is time, and  ⃗ is the body force vector exerted on the fluid.

Power and Torque Calculation
The lift and drag force components acting on the blade exert an equal and opposite reaction on the flow [54].This occurs at each point within the control volume , which is a hollow cylinder of thickness 4 with a radius equal to that of the rotor.This is used to calculate the instantaneous power output of the turbine at time . is blade length and d is span-wise blade element dimension.The total torque acting on the fluid within the hollow cylindrical volume  is The torque on the fluid acts in the opposite direction to the torque that turns the generator to create power  Uož and the torque due to the moment of inertia of the blades  cN , such that  PN = −Ÿ Uož +  cN .Here we have dropped the vector notation for torque, given that the torque vectors are all parallel to the -axis.For a fixed-speed turbine, Using the generator efficiency model from [56] to calculate power, we have where  ]^VN is the actual power,  ` is the drive train efficiency,  a is the generator and power conversion efficiency, and  O`^VN is the instantaneous power output of the turbine.

Torque Control and Thrust
As with the original WATTES, the moment of inertia of the rotor must be defined with torque to accelerate the blades in WATTES-V.Here, it is assumed the majority of each blade's mass is at distance , the rotor radius, from the centre of the rotor, and that each blade is identical to the other.
The moment of inertia for a vertical-axis turbine can then be written as where  cN is the number of blades,  is the mass per unit span, and  is the span length of each blade.We can then use  to define  cN , the torque that accelerates the blades.More details of this, and the time integration scheme used, can be found in [54].
The instantaneous thrust is calculated by integrating the  -direction body forces over the turbine control volume, that is

Turbine Parameterization
Due to the lack of an experimental prototype, the present vertical-axis turbine model is validated against data from wind tunnel experiments involving a two-bladed H-type vertical-axis wind turbine (VAWT) that was equipped with sensors to measure thrust and side loading on the turbine [70] turbulence model [71].We would like further to develop our vertical-axis turbine model by adding solid support struts as a conventional turbine, which would enable our model to be used to represent a wide range of vertical-axis turbines and turbine farms in the future.We thus chose to use a - SST model instead of a - model in this paper.Whilst there would be undoubted merit in exploring the effect of different turbulence models on the results, as undertaken by Barthelmie et al [72], this is beyond the scope of the present work, but is recommended for future study.The goal of the validation test is to check the ability of the newly developed numerical model WATTES-V to determine the thrust and side loading on the turbine for different values of azimuthal angle and tip speed ratio, with future applications to multi-bladed vertical-axis turbines in mind.This also enabled us to investigate the difference in behaviour between an idealised turbine with fixed tipspeed ratio and a more realistic turbine with torque control.(with straight side-walls not allowing the flow to expand) is likely to cause a blockage effect stronger than that in the experiments.The turbine is located 4.5 m downstream of the inlet, at mid elevation of the tunnel.Figure 3 shows a mesh slice in the  - plane, generated using blockMesh and snappyHexMesh utilities in OpenFOAM.The mesh is refined by a factor of 2 using a hexahedral mesh in a rectangular region containing the turbine and near-wake field, following [73].Here, mesh refinement is controlled by the number of cells in the ( • ,  • ,  i ) directions.Simulations were performed using the pimpleFoam solver, a merged PISO-SIMPLE algorithm.It should be noted that the azimuthal angle  used in [70] starts from the -axis, as indicated in the second figure in [70].In accordance with measurements from [70], the azimuth  described in the following sections has been transformed to the experimental coordinate system.

Results and Discussion
Initial and boundary conditions are selected to be approximate those in the physical wind tunnel test section.The inflow velocity is fixed at 4.01 m/s inflow.Lateral, bottom, and top walls of the computational domain are represented numerically by slip-flow conditions.A zero pressure gradient is applied at the inlet, and a fixed pressure prescribed at the outlet with zero gradients for other flow variables.Inlet turbulence intensity is ~10%, with turbulence kinetic energy  of 0.24 m X s X ⁄ and specific dissipation rate  of 1.78 s lW .It should be noted that the computational time for a simulation of ten revolutions was about six core hours for a parallel computation using four computing cores.

Validation and Grid Sensitivity Studies
Sensitivity studies concerning spatial and temporal resolution will be discussed in this section.
We first considered the convergence of turbine mean thrust coefficient for a tip-speed ratio of 3.3, shown in Figure 4. Mesh refinement is conducted by changing the number of cells in the -direction with a fixed cell aspect ratio and mesh topology.The relative error [74] between the results from the two finest meshes is below 0.5%, indicating that mesh convergence had been achieved.The spatial mesh resolution is hitherto set to 150 cells in the stream-wise -direction, with about 18 covering a single blade chord, (where the error between the finest mesh and the mesh employed is about 0.4%), giving a total number of 6.72 × 10 ² cells in the 3D simulation.Details of a mesh sensitivity study of the near-wake vorticity field are provided in the first part of Appendix B.

Two-Bladed H-Type Vertical-Axis Wind Turbine: Fixed Tip Speed Ratio
We now present results obtained for a two-bladed H-type vertical axis wind turbine where the tip speed ratio is set to a fixed value.Figure 5    We next study the effect of tip speed ratio on the mean thrust coefficient, carrying out numerical simulations that reproduce the experimental tip speed ratios of 2.7, 2.9, 3.1, 3. Figure 7 depicts the variation in force coefficients in the -and -directions with azimuthal angle for three selected TSR values.Amplitudes of both the predicted and measured force coefficients increase progressively with TSR.This is because the blade velocity and hence the relative flow velocity experienced by the blades increase as TSR is raised; the increased velocities then augment the blade load.There appears to be satisfactory overall agreement between the numerical predictions and measurements of  ³ and  • for TSR values of 3.3 and 3.7.However, there are more noticeable discrepancies between the predicted and measured values of  ³ and  • for TSR 2.9; this is because the angle of attack exceeds the critical angle for parts of each rotation when TSR is 2.9, causing stall to occur.Force coefficient in y-axis F y TSR-VATm: 3.7 TSR-exp: 3.7 TSR-VATm: 3.3 TSR-exp: 3.3 TSR-VATm: 2.9 TSR-exp: 2.9 Figure 8 illustrates the reduction in the lift coefficient that occurs at TSR = 2.9 as the critical angle of attack of the foil is exceeded at such a low value of TSR.It should be noted that data for cases where TSR < 2.5 were excluded from the experimental analysis because of this kind of poor aerodynamic performance [70].The higher TSR values (i.e.> 2.9) considered in the validation case are sufficiently large to be outside the range in which dynamic stall is likely to occur, and the predicted and measured values of  ³ and  • almost match.However, as the local velocity of the blades increases, so does the local Reynolds number based upon chord length,  n , which in turn affects the dynamic performance of the airfoil.In future work, the dynamic stall problem could be solved for the modelled vertical-axis turbines by controlling the blade pitch to attain an even or higher pressure drop along the whole diameter of a rotor.As shown in Figure 8, the angle of attack on each blade does not exceed 20° at TSR = 2.9 (and higher), and so the dynamic stall problem is not encountered here.The variance is calculated from: where  is a discrete random variable, E is an expectation operator,  S` is the total number of nodes in the region of the blade,  O is the probability mass function,  O is the local angle of attack at point , and  = E[] (or ) is the mean weighted value of angle of attack, given by and A similar method was used in an earlier study [56].The maximum variance occurs at the first rotor blade azimuth of 180 ° and the second blade azimuth of 0 ° or 360 ° (Figure 11).This threedimensionality might be due to shear flow instability, which is similar to that observed for a 2D pitching airfoil when its angle of attack decreases.The variance profiles are asymmetric with azimuthal angle, with large changes occurring after vortex detachment.turbulence kinetic energy profiles for TSR = 2.9, 3.3, and 3.7.U / U 0 - Figure 13 shows the lateral profiles of stream-wise mean velocity (Figure 13(a)-(c)) and turbulence kinetic energy (Figure 13(d)-(f)) at  = 1 − 5 downstream for TSR = 2.9, 3.3, and 3.7 respectively.The near-wake region (roughly   ⁄ < 2) is characterised by a low-momentum area isolated from the ambient flow in the presence of vortices, whereas the transition region (roughly 2 <   ⁄ < 5) is characterised by fast momentum recovery, high levels of turbulence, and expansion of the wake [80].In Figure 13(a)-(c), the asymmetry of mean velocity profiles is more visible closer to the turbine centre in the near-wake region.In Figure 13(d)-(f), the mean turbulence kinetic energy profiles are W-shaped.The two peaks are in accordance with those of the mean velocity profiles; however, the maximum peak of the turbulence kinetic energy is located on the side with negative , not on the side with positive  where the largest velocity deficit is observed.This is presumably due to the (aforementioned) large vortices that shed when the blade motion is in the same direction as the flow velocity in this area.These vortices play a key role, and affect mixing between the ambient flow and the low-velocity wake flow.Comparing the shape of these wake deficits with results from other published models of vertical-axis turbines [30,80], it can be stated that these characteristics of the mean velocity and turbulence kinetic energy profiles agree qualitatively with these previous studies of vertical-axis turbine wakes.For example, the shape of the mean stream-wise velocity profile of the present model corresponds well with those of experimental profiles presented in Figure 9 (left) in [30] and Fig. 5. (a) in [80], where the lowest values of   b ⁄ are both located close to  = 0.35.The shape of the turbulence kinetic energy profile exhibits good agreement with the experimental profile in Figure 9 (right) in [30], especially for areas in the vicinity of both peaks, and is in even better accordance than the University of New Hampshire reference vertical-axis turbine (UNH-RVAT) model used in [30].For the Edinburgh turbine (see Figure 1), the bending stresses at both ends are decreased by a factor of nearly four, with the red rings suppressing tip-vortex losses caused by the adjacent foils at different angles.Although the rings experience drag, the spoked wheel could well be a more efficient load-bearing structure than a tower, which experiences vortex shedding in addition to drag [65].We now present results obtained for a two-bladed H-type vertical axis wind turbine where the rotational speed of the blades is controlled by the torque.Figure 14 shows the limit cycle variation of the turbine angular velocity against azimuthal angle of the first blade, where the rotor is dynamically driven by the incoming wind flow.The predicted mean angular velocity  in the torque-controlled model is 17.87 / .This value is slightly smaller than that of the initial angular velocity (of  The torque coefficient  Í is calculated by using the dynamic generator torque data as

Two-Bladed H-Type Vertical-Axis Wind Turbine: Torque-Controlled Tip Speed Ratio
. Figure 15 shows a comparison of the torque coefficient results obtained for the fixed tipspeed ratio and torque-controlled cases.As shown in Figure 15(a), the predicted  Í for the model with fixed tip-speed ratio keeps changing through one rotor-revolution, whereas the  Í for the torque-controlled model remains almost constant with azimuthal angle.

Conclusions
This paper has presented a newly developed, efficient, parallelised, numerical model that simulates turbulent flow through vertical-axis turbines with a torque-controlled system, as well as with a fixed tip-speed ratio system.This computationally efficient numerical model WATTES-V of a behind close-packed contra-rotating vertical-axis tidal turbines [12]; hence, the support struts and tower shaft for a normal H-type vertical-axis turbine have not been considered herein.The present results show that vortex shedding occurs at the azimuthal position of the first rotor blade, at about 180 °.Vortices detach periodically from the turbine, and the resulting interactions create a complex downstream wake.The angle of attack for each blade did not exceed 20° in the present study, and so dynamic stall could be ignored.However, for future studies based on the present numerical model, either a dynamic stall model could be added as a correction, or a pitch-controlled system could be used to limit the angle of attack to an optimum value.
The wake field predicted by the present vertical-axis turbine model with fixed tip-speed ratio may be divided into two distinct regions.The near-wake region features a low-momentum zone where vortices shed from the turbine have a significant influence on the low-velocity region.The wake deficit in the transitional-wake region exhibits momentum recovery due to entrainment of ambient flow into the wake, and generates asymmetric velocity profiles about the wake centreline.
Analysis of wake turbulence behind a single vertical-axis turbine could facilitate better understanding of key flow features that contribute to wake recovery behind an array of close-packed contra-rotating vertical-axis turbines in future work.The sensitivity study on the turbulence parameters of the inlet flow and the downstream domain length (discussed in the Appendix B) should be useful for future experimental tests and numerical validations.
Dynamic predictions made by the present numerical model with torque-controlled tip-speed ratio are in satisfactory overall agreement with corresponding results from the fixed tip-speed ratio model and experimental data [70] on thrust and lateral loading.In the former case, the rotor is demonstrably driven by the blade-generated lift, which is counteracted by the torque that accelerates the blades and turns the generator.The present model should be useful in the future by enabling predictions of the dynamic response of practical vertical-axis turbines to unsteady flow.

Model Architecture
The Wind and Tidal Turbine Embedded Simulator (WATTES) [56,54] code is an open library source code written in Fortran 95, which employs both the dynamic torque-controlled actuator disc and the actuator line methods with active-pitch correction to simulate the behaviour of multiple wind and tidal horizontal-axis turbines, together with a simplified generator model.Compared with other momentum codes, WATTES predicts the dynamic response of the device to the flow, with lift and drag force components balanced by inertial effects and the resistive torque induced by the generator.
Force components are incorporated within the incompressible Navier-Stokes momentum equations as body force components [54].For computational efficiency, WATTES exploits parallel programming based on multiple instructions multiple data (MIMD) [52]

Appendix B
This Appendix presents results from tests which examine the influence of mesh convergence on the vorticity field in the near wake, the choice of inlet turbulence parameter, and the length of the downstream domain dimension.different meshes, although some slight discrepancies are evident, the relative two-norm errors [56] are 2.80%, 2.44%, 1.94% respectively, which are all under 3% and are within acceptable margins.

Effect of Mesh Convergence on Near-Wake Vorticity Field
We find that a spatial grid resolution of 150 cells in the  -direction, giving a total number of 6.72 × 10 ² cells in a 3D simulation, is sufficient to achieve mesh convergence.

Sensitivity Analysis concerning Inlet Turbulence Parameters
Turbulence intensity (TI) is defined as the ratio of the root-mean-square of flow velocity fluctuations  ?≡ Ö W ® Ÿ • ?X +  • ?X +  i ?X to the mean flow speed  ≡ Ö • X +  • X +  i X [48], and is expressed: where  is turbulence kinetic energy (TKE).A value for TKE at the inlet is thus calculated from Equation (B1) for a given TI [82].The specific dissipation rate  used in the - SST turbulence model in OpenFOAM is calculated using the following formula [83]: where  à is a turbulence model constant equal to 0.09, and  is the turbulence length scale.Figure B2 displays the variation of specific dissipation rate with turbulence kinetic energy, and the thrust coefficient with turbulence intensity for the vertical-axis turbine model.It can be seen that the mean thrust coefficient tends to decrease as TI increases.In particular, as TI varies from 0.1% to 20%, the thrust coefficient decreases by 6.34%.This indicates that the choice of level of turbulence intensity at the inlet can have a substantial effect on the thrust value of a vertical-axis turbine.
However, the change of thrust coefficient is only about 1.72% for a more realistic range of TI between 1% and 10%.

Sensitivity Analysis concerning Downstream Domain Size
To investigate the impact of the limited downstream domain size on the results, we doubled the stream-wise length of the downstream domain, for a case of fixed TSR = 3.3, and compared the thrust and lateral force coefficients obtained using the two domains.Figure B3 shows that very satisfactory agreement is obtained for the values of  ³ and  • on the two domains, for a fixed TSR = 3.3; relative errors between the coefficients obtained using the different domains lie below 0.057%.This confirms that the downstream length utilized in the main paper is adequate.
shows a schematic diagram illustrating a turbine blade with chord, pitch, and path of a single blade.The diagram also indicates the force component vectors that provide loading on the blade.The black dashed circle represents the circular trajectory of a blade.

Figure 2 .
Figure 2. Geometry of and force vectors on a blade of a rotating vertical-axis turbine (VAT).The flow velocity relative to the blades is  ]^N 6666666⃗; the angle of attack  is calculated from the local inflow velocity  6⃗; the freestream velocity  b 6666⃗; and the blade velocity is  cN 66666⃗.The azimuthal blade angle is  with the corrected blade pitch ; and  ]^N is relative angle. M and  ( are lift and drag forces per unit span respectively for the actuator line.
with  O and  O the local coordinates of blade ,  and  are the point coordinates, and the standard deviation  determines the width of the Gaussian kernel.
. The experimental data were collected at the Open Jet Facility at Delft University of Technology[70], which comprised a closed loop open jet air flow of 2.85 m × 2.85 m outlet cross section.The wind tunnel test section was 13 m long.Table1lists the turbine model parameters, derived from[70].The numerical model neglects the rotor shaft and support struts, and utilizes an unsteady Reynolds-averaged Navier-Stokes (URANS) formulation with  - shear stress transport (SST) turbulence closure scheme in OpenFOAM.The URANS approach is an attractive, computationally inexpensive prospect for far-wake simulation[55].The  - SST turbulence model used is the original Menter model[71], which has been used successfully for many different types of flows.The SST (shear stress transport) turbulence model combines the - model in the free shear flow, with the  - model in the near wall boundary regions.It is a robust two-equation eddy-viscosity

Figure 3 .
Figure 3. Computational mesh and boundary conditions, showing plan dimensions of the modelled domain.

Figure 4 (
b) displays timestep resolution test data, evaluated on the 3D grid with 150 cells in the -direction.The relative error is below 0.5%, indicating low sensitivity to temporal resolution.In all these convergence tests, the Courant-Friedrichs-Lewy (CFL) number[75] is below 0.58.In this study, we employed ∆ = 0.03 s, corresponding to 120 time steps per revolution, giving a CFL number of 0.23.Simulations were carried out lasting at least 10 revolutions, with periodic convergence reached after 9 revolutions when the difference in maximum turbine thrust between successive revolutions was 0.06%.

Figure 4 .
Figure 4. Resolution sensitivity of the VAT model: (a) Spatial resolution after 120 time steps per revolution; (b) Temporal resolution on a mesh with 150 cells in -direction.
compares the numerical predictions and measured thrust and lateral force components on the rotor for an incoming flow speed of 4.01 m/s, a fixed pitch angle of 0 °, and a tip-speed ratio (TSR) of 3.7.The measurements were averaged over 22 turbine rotations.It can be seen that the numerical predictions and experimental measurements of the force -and -directions are similar in terms of amplitude and profile, with the maximum thrust loading experienced at the blade azimuth at 90 ° and 270 °.

Figure 5 .
Figure 5.Comparison between predicted and measured [70] thrust and lateral forces on a wind turbine rotor of a two-bladed H-type vertical-axis wind turbine for an incoming flow speed of 4.01 m/s, fixed pitch blade angle of 0 °, and tip speed ratio of 3.7.

Figure 6 .
Figure 6.Comparison between predicted and measured [70] mean thrust coefficient  ³ as a function of tip speed ratio in the range from 2.7 to 3.7, for a wind turbine rotor of a two-bladed H-type vertical-axis wind turbine with incoming flow speed of 4.01 m s ⁄ and fixed pitch blade angle of 0 °.

Figure 8 .
Figure 8. Lift coefficient as a function of local angle of attack at each grid point predicted by the VAT model for TSR = 2.9, compared with measurements for a static airfoil [76]; the red dashed lines show the range of local angle of attack.

Figure 9 .
Figure 9. Flow patterns at eight different phases during a single revolution for TSR = 3.7: (a) Velocity magnitude (m s ⁄ ); (b) -component of vorticity (s lW ); and (c) Turbulence kinetic energy (m X s X ⁄ ) in the central horizontal - plane.

Figure 9
Figure 9 shows plan views of the evolving velocity magnitude, vorticity -component fields, and turbulence kinetic energy contours in the horizontal - plane at eight different phases during one revolution of the 2-bladed VAT operating at TSR = 3.7.The white blades are shown for interpretation only, and the - SST model behaves like a - model even near the blades (as in the free shear flow). is the azimuthal angle of the first blade measured from the experimental turbine in the anti-clockwise direction.The incoming flow passes through an annulus mapped out by the anti-clockwise rotating turbine, with vorticity generated on the surface of the blade and a turbulent wake developing downstream.The rotor interacts with its own wake, especially for azimuthal angles of 90 ° and 270 °, causing the thrust to increase.Vortex shedding starts to occur when the first rotor blade reaches an azimuthal position of about 180 °.Vortices detach periodically from the turbine, and move to the downstream low-pressure wake field.This vortex shedding process drives oscillations in the local flow field affecting the forces on the rotor blades.The highest turbulence kinetic energy is observed at about 90 ° or 270 ° of the azimuthal position.

Figure 10
Figure 10 illustrates an instantaneous three-dimensional vorticity field around the turbine.It can be seen that a smooth, quasi-two-dimensional shear layer, as a consequence of using URANS turbulence modelling, is created behind a blade moving towards the upstream direction.The blade then turns into the downstream direction and sheds large and more three-dimensional (spanwisemodulated) vortices.Strong tip vortices then interact with the shed vortices, and create a complex downstream wake field.

Figure 11 .
Figure 11.Variance of angle of attack as a function of azimuthal angle at TSR = 3.7.

Figure 11
Figure 11  shows the variance in angle of attack experienced by the two blades during a single revolution.This arises due to the blade shedding sheet vortices, which then break up into threedimensional turbulence when the blade moves towards the downstream direction, giving greater

Figure 12 1 .
Figure 12 presents an overview of the downstream wake evolution behind the turbine with the distribution of the mean stream-wise velocity and the turbulence kinetic energy in the near-wake region at   ⁄ = 1.The mean velocity field for TSR = 3.7, shown in Figure 12(a), is obviously asymmetric in the transverse () direction.The mean wake deficit in Figure 12(c) describes the characteristic of the mean velocity as it recovers rapidly on the coarse mesh [77-79] for the three selected tip speed ratios.Minimum values of mean velocities were found predominantly to occur at  ~ 0.35 .In the bypass flow at   ⁄ > 1.5 , the stream-wise velocity component reaches approximately |  b ⁄ | = 1.1, due to the blockage effect.Turbulence kinetic energy profiles in the vicinity of the rotor also exhibit clear asymmetry, with a peak at   ⁄ ~− 0.2.

Figure 14 .
Figure 14.Variation of angular velocity  as a function of azimuthal angle in the torque-controlled system at TSR = 3.3.

Figure 15 .
Figure 15.Variation of torque coefficient  Í as a function of azimuthal angle at TSR = 3.3: (a) Comparison between results from the model with fixed TSR and its torque-controlled counterpart; (b) Enlarged zone of  Í variation in the torque-controlled model.

Figure 15 (Figure 16 .
Figure 16.Comparison between predicted (including the model with fixed TSR and the torquecontrolled model) and measured [70] thrust and lateral force coefficients for TSR= 3.3 as functions of azimuthal angle: (a) thrust coefficient  ³ ; and (b) lateral force coefficient  • .

Figure 16
Figure16displays the thrust and lateral force coefficients as functions of the azimuthal angle obtained from the fixed tip-speed ratio model, the measured data[70], and the torque-controlled single cross-flow turbine was developed within the OpenFOAM CFD framework.The model is based on actuator line theory, and combines classical blade element theory, an unsteady Reynolds-averaged Navier-Stokes flow model, and a - SST turbulence model.This numerical model with fixed tip-speed ratio was validated against experimental data acquired from an H-type 2-bladed vertical-axis wind turbine [70].The model gives numerical predictions in satisfactory overall agreement with experimental data on thrust and lateral loading.It is planned that the present cross-flow turbine model will be employed in future research on wakes through the Message Passaging Interface protocol (MPI).The solution is computed on a number of processors that function asynchronously and independently.The original WATTES model simulated flows using Fluidity, which is an open-source hr-adaptive multiphase computational fluid dynamics (CFD) solver based on an unstructured finite element method and offers anisotropic mesh refinement, developed mainly by researchers at Imperial College London [81].The original WATTES source code was used to represent horizontal-axis turbines within the OpenFOAM [65] CFD framework, and formed the basis of the modified numerical model WATTES-V used herein to simulate flow past a vertical-axis turbine.OpenFOAM is freely available open-source CFD software based on the finite volume method on general unstructured polyhedral meshes, and is written in C++.In order to benefit from the advantages provided by the original WATTES source code, proper coupling of WATTES and OpenFOAM was a necessary prerequisite before the further development of WATTES-V model described in the present study.

Figure
Figure B1 presents horizontal profiles of turbine mean streamwise velocity at a tip-speed ratio of 3.3 in the near-wake region computed on coarse, medium, and fine meshes (with  • = 150,  • = 180, and  • = 375 cells respectively in the x-direction).The figure illustrates model sensitivity to spatial resolution.Satisfactory agreement is generally achieved between the profiles obtained on the Sensitivity of the results to the inlet turbulence parameters is examined by setting different inlet values of  and  calculated from Equations (B1) and (B2) for a range of turbulence intensity values from 0.1% to 20%, with TSR = 3.3.The results are shown in FigureB2.

Figure B3 .
Figure B3.Comparison of predicted force coefficients obtained on meshes of downstream length 6 and 12, for TSR = 3.3 as functions of azimuthal angle: (a) thrust coefficient  ³ ; and (b) lateral force coefficient  • .