A spatially nonlinear generalised actuator disk model for the simulation of horizontal axis wind and tidal turbines

Efficient numerical simulation of renewable energy wind and tidal turbines is important for the layout of devices in farms. Computational Fluid Dynamics (CFD) approaches using blade geometry resolved models are computationally expensive. Therefore, most array models use source term representations of rotors, normally actuator disk, actuator line or blade element disk. Unfortunately, these methods rarely capture enough physics to accurately predict power and at the same time correctly characterise the wake velocity field and turbulent structures. This study describes a new Generalised Actuator Disk CFD model (GAD-CFD), that achieves the required accuracy for the simulation of horizontal axis wind and tidal turbines and their wakes. This new method combines a finite volume CFD code with additional source terms representing the rotor, including: correct consideration of losses along the foil by modification of the distribution of downwash; a concise downwash distribution computation; recognition that foil cross section varies along the length; dynamically changing Reynolds numbers and the application of a tip radius correction. Also reported are foil lift and drag coefficients and their variation with thickness, surface roughness and Reynolds number, which is necessary for the proper characterisation the whole rotor. The effectiveness of this approach is investigated and validated against two experiments, and demonstrates improvements over traditional source term methods, in particular the correct CFD approach to tip losses and consequent downstream wake prediction. This study provides confidence in application to both small scale flume studies and large scale array deployments in both the marine and wind


Introduction
The success in recent years of the renewable energy sector in demonstrating significant progress developing offshore wind and marine energy technologies are of particular interest at the current time as the availability of offshore regions for device development and deployment increase.
While wind turbines are migrating from onshore deployments to offshore locations aiming for increased energy yield, the marine energy sector is progressing renewable energy generation technology based on the predictability of tidal flow [1]. The United Kingdom has significant wind [2] and tidal [3] energy density due to its island status and location.
To capture these resources, substantial deployment of large scale horizontal axis tidal (HATT) and wind (HAWT) turbines will be required. The growing demands of these sectors require investments in improved knowledge and technology [4,5]. This can be achieved through the study of: analytical modelling, numerical modelling, practical experiment, and deployment [6e8].
Academic interest in the sector is growing in parallel to the growth of industrial investment. The range of work includes; device design studies [9e12], environmental impact studies [13e15], and the study of the performance of individual devices [16e19] arrays [20e22] and array optimisation [23e25]. The computational efficiency of the approach presented here will make such optimisation procedures more useable for industrial applications.
Computationally efficient CFD models of arrays typically use actuator disk, actuator line or blade element disk approaches as outlined below in the literature review. In all these methods the absence of fully resolved turbine foils in the flow regime requires assumptions regarding the flow characteristics to predict the correct turbine operating conditions. This includes lift and drag characteristics for a full set of Reynolds numbers, changing foil cross section at different radial stations, the capture of a radiused tip (due to design or manufacturing processes), and improved prediction of the distribution of lift and drag along the radial length of the foil for a range of different chord, twist, and section geometries.
The work presented here seeks to validate a new Generalised Actuator Disk CFD model (GAD-CFD) which addresses the issues described above. Each cell in the defined disk utilises the local flow velocity to compute the rotor forces, so can deal with boundary layer shear flows and other non-linearities in the upstream flow. The radial variation in blade design is captured including changes in chord, twist, sectional profile and Reynolds number effects. The model introduces an analytically derived lift distribution method to determine the correct downwash, and thus correct tip loss treatment in the context of a CFD model. This work combines the finite volume CFD code, OpenFOAM [26], with additional energy source terms representing the rotor. The work is correlated with experiments by Mycek et al. [16] (water), and Selig et al. [27] (air).
A review of the literature appropriate to this study is presented in Section 2. This work introduces the detailed model in Section 3. Section 4 is important in its own right, as this discusses the relationship between foil lift and drag and Reynolds number. It is necessary to understand these effects to correctly match the experimental results. A set of case studies are detailed in Section 5. The cases studies simulate experiments from the literature and are chosen to represent a range of scenarios which can be used to validate the model. Section 6 discusses the results of the simulations while comparing and contrasting with the experiments. Conclusions and potential future work are highlighted in Section 7.

Literature review
Ocean energy resource evaluations can be conducted at differing scales. Large oceanographic models simulating turbine energy extraction as sub grid drag coefficients can be used to approximate sites of interest. At this scale resource modelling provides greater coverage than can be achieved with field measurements [28]. Sediment transport assessments [29], and the examination of waves and currents [30], can be studied at this scale.
Within a three dimensional hydrodynamic model to simulate a tidal turbine farm at the Alderney race, [29] implements an additional bed friction term. An extension to this type of model [31] implements rotor source terms in the layers of a 3D coastal model. [25] uses the same approach within an array optimisation. [32] used a two dimensional depth averaged model to predict an array of turbines and structures where the turbines and structures are represented as additional thrust and drag forces. Using the same model, [33] evaluates the effect of array shapes. The authors conclude that denser arrays have more significant, localised effects on water level and suspended sediment concentration.
Energy source models attempt to abstract the rotor geometry into an energy source/sink term, [23] employed an actuator disk approach through a porous jump boundary condition in steady state 3D simulations, while [7] uses multiple actuator disks to represent a set of rotors, investigating the positive effects on performance that blockage can induce with small lateral intervals. [34] abstract this to an actuator fence model, and study the interactions with a simple headland. This work supersedes a more detailed model investigated by [35], who simulated a lateral turbine fence using a 2D CFD model, which showed that headland flow acceleration can be utilised to improve the performance of turbines arrays.
The ability to predict the effect and utility of deployed wind/ tidal farms is important for the viability analysis of a site layout/ location or the design and performance of an individual device. The prediction of this information is essential to the turbine designer. When data from offshore deployments is not available, laboratory scale experiments in wind tunnels, current flumes or wave/current tanks [36e38] can provide adequate data for single or small scale arrays.
A number of experimental studies have been published which can be used for validation and these include works by Selig et al. [27], Bahaj et al. [39], Stallard et al. [40], and Mycek et al. [16]. Selig et al. [27] study the performance of a full scale device at the National Renewable Energy Laboratory (NREL). The device is a three bladed horizontal axis wind turbine generator instrumented to provide pressure data along the aerofoils as well as performance of the full turbine. Bahaj et al. [39] used a laboratory scale tidal stream turbine in a practical study conducted in a cavitation tank and a tow tank to simulate the performance and cavitation effects. Stallard et al. [40] study a group of three-bladed tidal turbines arranged in a flume, noting the effects of wake interaction and recovery. Mycek et al. [16] conducted studies of a well defined rotor geometry situated in a recirculating flume, including a dual rotor layout [21].
Laboratory studies provide valuable insight into the localised performance of wind and marine turbines. The experiments provide detailed data for analysis, and correlation of different modelling techniques. Computational fluid dynamics (CFD) can generate predictions ranging from close approximations to highly accurate simulations of real world scenarios. However with high accuracy comes computational cost. The simulation of entire fields of full scale turbines at a useful level of approximation for studying turbine performance, interaction, and environmental effects is very challenging and costly [41,42]. However, the cost is still significantly lower compared to full offshore deployments. Detailed computational models of individual turbines may include moving meshes, which are modelled to the geometry of the aerofoil. These models resolve aerofoil tip vortices, wake propagation, and turbine performance characteristics. Turbulence is resolved using either Reynolds Averaged Navier Stokes [35] or Large Eddy Simulations (LES) [43].
Addressing the computational cost of transient CFD models, Blade Element Theory (BEMT) approaches were developed from Glauert's propeller theory [44]. BEMT is an analytical approach to solving for forces at the rotor utilising inputs from tabulated rotor geometry and lift/drag data [44,45]. This method is utilised by the wind industry as it is simple, fast and provides good agreement with experiments [46]. This model is also utilised by the marine energy sector [47,48], incorporating corrections for tip and hub losses, yawed and heavily loaded rotors. It also gives a fast assessment of the effects of turbulence [49]. BEMT methods are not able to predict the local velocity field, and therefore the downstream wake. For this a full transient CFD simulation would be required to predict the downstream velocity field, and thus device and environmental interactions. This approach, however, can be prohibitive in terms of computational cost when simulating multiple interacting devices in a deployed scenario.
Blade Element Momentum CFD (BEM-CFD) steady state models [50,51], provide methods of more efficiently simulating HATTs and HAWTs. These methods utilise a radially varying set of blade characteristics, uniformly distributed in an axial direction. Hence, computational cells at the same radius from the rotor centre have the same properties, however, as the flow varies from cell to cell, the resultant forces on the fluid also vary.
A steady state CFD model coupled to a BEMT code [52] was based on the classic Blade Element Theory model [45], to predict the rotor effect on the fluid domain. The authors then examine the power output of rows of turbines where each row is reduced in performance by a constant factor. This type of model has been used to study wake length and inflow characteristics [53], the effect of accelerating flows [54] and the interaction of arrays of turbines [20], validated against [55]. The application of tip loss corrections pertinent to a CFD type model representation takes this approach further [17]. The benefit of a steady state approach, with respect to computational cost, is unequivocal. This technique allows us to move into the realms of array interaction modelling and site design at a more reasonable level of cost.
Actuator line approaches have also been studied where a rotating line of energy sinks represents each of the rotor blades. The method was developed by [56], and further studied by [57] who validated their approach against the NREL 10 m diameter wind turbine experiments [27]. An actuator disk and line approach can be combined with a LES CFD model [58]. These actuator line models do not treat tip loss or downwash effects. The effects of support structures of contra rotating turbines is examined by [22] utilising LES-CFD with an embedded actuator line model.
The extensions discussed by [17] better capture the effects of tip vortices within the steady state model, and is validated with experiments [16]. [59] discuss the idea of correction coefficients to the Shen et al. [60] tip loss functions applied to an actuator line CFD model. These corrections are consistent with the experimental data and are part of an extensive investigation into wind turbine modelling techniques with an important emphasis on tip loss mechanisms [61,62].
The tip loss distribution of [17] follows the classical work of Prandtl and assumes a simple elliptical distribution, which although a good approximation for the geometry specified in [16], can demonstrate significant inaccuracies for other blade geometries, particularly for alternative chord distributions. In addition, the assumption that devices are run in flow conditions independent of the Reynolds number (especially laboratory scale studies), is inaccurate. Lift and drag characteristics significantly change with lower Reynolds numbers and higher turbulence, thus demonstrate varying results when the ratio of thickness to chord length changes along to radial length of the foil.
It can be seen from this extensive body of literature that energy source term approaches to modelling turbines can provide a good estimate of turbine power and thrust. All of them provide an approximation of the downstream wake and this is where the estimates vary. Actuator disks models are the most simple to set up, but do not impart swirl on the wake. Blade element disks correctly predict swirl, but as the effect is rotationally averaged, cannot produce a tip vortex. Actuator lines clearly show a tip vortex and the characteristic corkscrew wake pattern. In all cases, the tip treatment is important for performance and wake extent; performance losses and creation of turbulence in the wake edge must both exist and be numerically consistent. Therefore, in the next section we describe such an approach.

CFD and the governing equations
The OpenFOAM toolbox [26] is utilised for the model implementation. The CFD model requires the solution of the Navier Stokes equations for mass and momentum, where the disc rotor characteristics are incorporated into the momentum equation through an additional source term S i . O'Doherty et al. [63] discuss a range of Reynolds Averaged Navier Stokes (RANS) turbulence closure models. The results of this work show the Reynolds stress model (RSM), Realisable k-ε, and the k-ε RNG variant all demonstrate good fit with experimental data. The effectiveness of the k-ε RNG model [64] for use with large rotating downstream wakes, and ability to effectively capture downstream wake structures, is demonstrated in the work by [17]. The k-ε RNG model is implemented here. Fig. 1 illustrates how a turbine with three hydrofoils is discretised for use with the generalised actuator disk method. The hydrofoil properties are determined at radius r i , and then averaged over 2p radians. This process is repeated for each hydrofoil element over the interval ½r 0 ;r max . Full details of the underlying method are presented in Ref. [17]. The procedure leads to the definition of the following axial and tangential source terms, S a and S t ,which are incorporated into the source term S i in the momentum equation.

The GAD-CFD method
Here C L and C D are the lift and drag coefficients respectively, r is the density, c is the chord length and v R is the resultant velocity.

An extended downwash distribution method
To improve the prediction of the effect of tip losses on the flow field an additional source term, representing the tip vortex induced downwash w (see Fig. 2), is computed. The downwash force is proportional to the force deflecting the flow around the foil (S a and S t ), weighted by a downwash distribution function EðrÞ.
The use of the elliptical downwash distribution EðrÞ2½0; …; 1, where the normalised distribution is 0 at the hub and 1 at the tip, is described in [17]. The authors present an ellipse as a reasonable approximation for a tapered foil design. However with a range of different blade chord distributions which do not conform to an elliptical downwash distribution approximation, a more robust distribution method is proposed. The presented technique utilises the foil's geometry to generate a more accurate representation of the downwash distribution. This is achieved by analytically solving a set of linear equations based on Prandtl's classical lifting line theory, as described in ([65] Section 5.3.2 General lift distribution). The spanwise coordinate transform: is a function of q, where b=2 is the foil radius and 0 q p. The elliptic downwash distribution as a function of q is: This equation can be solved using a Fourier sine series to represent the circulation distribution along a fixed length foil. Thus it can be presumed generally that: where the quantity of N terms defines the required accuracy of the solution, and coefficients A n (where n ¼ 1; …; N) are to be solved.
However, the A n 's must agree with Prandtl's lifting line equation: where this equation is evaluated at a given spanwise location q 0 . Additionally, b, cðq 0 Þ, and a FL¼0 ðq 0 Þ are known values derived from the foil section geometry. Only the A n 's need to be determined.
Thus, defined at a given spanwise location (i.e. a specified q 0 ), this equation is one algebraic equation with N unknowns, A 1 ;A 2 ;…;A n . A set of N independent algebraic equations with N solvable unknowns is evaluated. Once the A n coefficients are determined, the distribution of downwash (GðqÞ) can be calculated by substituting the A n 's in Equation (5). Examples of the normalised computed downwash distribution curves for a number of blade profiles is shown in Fig. 3. It is observed that the distributions vary significantly with respect to the underlying geometry.
The source term S i is appended to the momentum equation and defined thus: where the scalar terms S a and S t are the magnitude of kinetic energy per unit volume in the axial and tangential directions as defined in Equations (1) and (2), S v represents the additional downwash force, and b v a is the axial unit vector, b v t is the tangential unit vector, and c v w is the unit vector normal to a FL¼0 , on the plane S i is an R 3 representing the change in kinetic energy per unit volume in Cartesian space, i.e. dynamic pressure.
Combining the downwash distribution function with the predicted axial and tangential forces, the downwash term can be described as follows: thus: which makes for both a simpler and computationally less expensive implementation. From the above formulation, power and thrust is determined as follows: and where n is the set of cell indices, i, bounded by the bladebox (the domain associated with the rotor), and Vol i is the volume of cell i. Note that Equation (11) is the discretised form of the volume integral of the tangential force multiplied by rotational speed u and is hence the useful power generated.

Additional model refinements
This subsection discusses the additional refinements included in the GAD-CFD model. The purpose of these refinements is to fine tune the predictive capabilities of the model, with the aim of improving model flexibility and simulation robustness.

Tip radius treatment
The tip radius r t referred to here is the radius found at the end  face of a foil which blends between the foil section edge and the end face, see Section C-C in Fig. 4. This is usually provided for the purposes of manufacture (sharp edge removal etc.), and is quite small. However, on some designs the radius in this location can become quite large, and thus effect the tip losses. In this region the lift is effectively zero, while the drag is maintained. To capture this effect within the model the downwash distribution constraint of a FL¼0 is offset inboard by a quantity x which is equal to the specified tip radius, r t , as demonstrated in Fig. 4.

Radially changing foil sections
In a number of foil design cases, the foil section thickness may vary relative to the section chord, i.e. the thickness as a percentage of chord length may vary along the radial length of the foil, see example in Fig. 4. In some cases differing foil sections are applied to the design along the radial length. This can have an effect on the power extraction characteristics, in particular at low rotor speeds. Within this model these section changes are captured through the use of C L and C D lookup tables as a function of radial location linearly mapped between section types. The use of these look-up tables is discussed further in Section 4.

Aerofoil characteristics
This section discusses the relationship between Reynolds number, surface roughness effects and foil performance characteristics. It is necessary to understand these effects to correctly match the experimental results.
Input lift/drag coefficient data, plotted against angle of attack for a full 360 , is taken from a number of sources and tabulated in an ascii text file for use as input to the model. The lift and drag data for given foil sections are generated using a combination of JavaFoil [66] and CFD simulation.
The foil sections used for the simulations range from NACA63418 to NACA63424, and the S809 wind turbine foil section. JavaFoil is utilised for the NACA foil sections, while CFD simulation of the S809 foil section is used as JavaFoil does not include this profile in its standard configuration.
Lift and drag data is required for all these sections for the range of angles of attack 0e360 . The data sources described above provide data approximately in the range of plus or minus 20 . Thus additional information derived from simple flat plate theory is used to complete the range. A simple cubic blending function is used to splice the data into a coherent 360 range.
The majority of experimental lift and drag data is derived at higher Reynolds numbers suitable for aircraft. However, data for the range of Reynolds numbers suitable for flume scale experiments is not available, thus JavaFoil is used to provide this data. JavaFoil demonstrated good correlation for the NACA63418 to NACA63424 range of foils as compared to the experimental work by [67] for identical Reynolds numbers. The possible exception to this is in the stall region, where full stall is predicted to occur at a slightly higher angle of attack and thus reports slightly higher lift values in this region, see Fig. 5. As S809 data from JavaFoil is Fig. 3. Normalised spanwise downwash distribution comparison of the case presented by: Mycek [16], Selig [27], Bahaj [39], and the elliptical distribution from [17]. Fig. 4. Demonstration of section changes along the radial length of the foil, i.e. Section A-A and Section B-B. The treatment of the tip radius is also highlighted, i.e. the offset of the downwash distribution x is equivalent to the radius rt at the foil tip as shown in Section C-C. unavailable, a CFD simulation of the S809 foil section is conducted for the correct Reynolds numbers, however the performance in the stall region is slightly inflated in the same fashion as the JavaFoil predictions discussed above, see Fig. 6.
The lift and drag curves are generated at a set of Reynolds numbers capturing the operating range of the turbines, combined with the NACA standard roughness parameter [67]. The Reynolds numbers are computed based on foil chord length, and consider both the rotational speed of the turbine and the inflow velocity to compute the actual velocity, i.e. jv R j in Fig. 1. This combined with a range of foil sections matching the foil design (at given radial stations) allows us to generate a 3D lookup table accessed as a function of Reynolds number (Re), radial station(r) and angle of attack (a) for both the lift and drag values, e.g. C L ¼ f ða; Re; rÞ.

Reynolds number treatments
Another effect on the performance of the GAD-CFD model is the correct treatment of Reynolds number linked lift and drag characteristics [67]. When simulating deployed cases such that the rotors are not running at Reynolds independent speeds, i.e. the chord based Reynolds number linked lift and drag input data, the simulations will incorrectly predict performance. Not capturing this phenomena within the model can skew the determined power and thrust characteristics. The model incorporates an on the fly calculation of chord based Reynolds number for every cell of interest, and uses this in conjunction with the 3D lookup tables for C L and C D that have been generated using the approach described in Section 4.

Standard surface roughness considerations
The quantity of surface roughness is an important consideration for foil performance as it governs boundary layer thickness, and flow separation characteristics. The NACA standard for surface roughness is described in [67]. These characteristics in turn affect the overall quantity of lift and drag, and also the stall point [65]. Unless the surface of the foil is highly polished, then this performance degrading phenomena should be included in the 3D lift and drag lookup data tables. Of course, when attempting to compare with scale model experiments, the rotors can be smaller than desired and so can be relatively rough in relation to the chord length.

Case studies
The model described in Section 3 is implemented in OpenFOAM [26]. This section describes computational model setup for two case studies including: mesh configuration, boundary conditions, and initial conditions, with respect to the OpenFoam Toolbox.

Mesh configuration
Using the blockMesh utility, the domain is configured as a hexahedral mesh with bounding points defined in Table 1, which captures the domain extents of 18 m Â 4 m x 2 m in x, y, and z. The domain is then subdivided using divisions defined in Table 2. Simple grading is used to refine the mesh at the left, right, and bottom boundaries, see Table 3. With reference to Fig. 7, refinement of the mesh around the turbine and wake region is achieved using the 'snappyHexMesh' utility.
The wake region is defined as a cylinder 0.7 m radius, extending from the rotor to 9 m downstream, i.e. to the domain outflow. The refinement level in this region is specified as level 2 i.e. subdivide the base cell/mesh twice in this region (Note 2 in Fig. 7).
The bladebox region is an annular disk with outer radius equal to the rotor radius, and inner radius equal to the nacelle radius. The bladebox thickness is chosen to be slightly larger than the maximum chord length of the rotor. The rotor assembly and bladebox region are generated as.stl files and imported into the meshing tool. Within 'snappyHexMesh' a refinement level of 5 is used for the rotor assembly and bladebox, i.e. Note 5 in Fig. 7. The region around the rotor assembly (Note 3 in Fig. 7) is set at level 3 up to 0.6 m from the assembly. The final mesh sizes used for the mesh sensitivity study are shown in Table 4. A reasonable level of detail of the nacelle and support is included in the model as shown in Fig. 7.

Initial and boundary conditions
At the inlet a plug flow velocity condition is imposed (see Section 6 for values). A reference pressure of 0 Pa is imposed at the outflow boundary. No-slip boundary conditions are applied to the base, side walls and turbine supporting structure. A log-law wall function is used to compute the boundary layer effect. The kinetic energy initial condition and inlet is set to 0:0019m 2 :s À2 (for the 3% Turbulence Intensity (TI) case), while the dissipation rate initial condition and inlet is set to 0:000097m 2 :s À3 .

Mesh sensitivity study
The study of a mesh with respect to mesh independent results is important for establishing an efficient computational discretisation, while maintaining enough finesse to capture areas of more complex flow behaviour. Within the mesh a region is defined such that it captures the swept volume of the rotor. This refined region supports the computation of time averaged foil forces with respect to the GAD-CFD coupled model. An additional region is defined to capture the downstream wake emanating from the rotor and support structure. The mesh generation process is described in subsection 5.1.1, and the model is solved with parameters set in subsection 5.1.2.
For the mesh independence study an examination of the coefficients of power C P and axial thrust C T is performed across the range of meshes shown in Table 4 from 25,486 cells to 3.7 M cells. Fig. 8 demonstrates that the coefficients tend to converge around Mesh 2 (0.3 M cells), and remain stable thereafter. Fig. 9 reports wake velocities at six locations across the computational domain. The locations are displayed in Table 5. The velocities converge around Mesh 3 (0.9 M cells), remaining mesh independent thereafter. Thus, the results reported in Section 6 are obtained using Mesh 3. This mesh configuration is also utilised for Case Study 2 (subsection 5.2).

Case Study 2: NREL phase III wind turbine
In this subsection the case setup for comparing the GAD-CFD model with the experimental work of [27] for a single turbine is described.

Mesh configuration
The mesh characteristics are defined according to the mesh sensitivity study presented in subsection 5.1.3. Using the block-Mesh utility, the mesh domain is configured as a hexahedral mesh with bounding points defined in Table 6, which captures the domain extents of 18 m Â 4 m x 2 m in x, y, and z. The domain is then subdivided using divisions of 80, 37, and 28 respectively. Simple grading is used to refine the mesh at the left, right, and bottom boundaries, see Table 7. Table 1 Initial block extents, in metres, representing the overall domain.  Table 3 Mesh grading in the x, y, and z directions using the studied configurations.

Start
Middle End  [16]. case using a combination of 'blockMesh' and 'snappyHexMesh' utilities. Note 0 shows the outer/base distribution of cells, while Note 2 shows the level 2 refinements made in the wake region. Note 3 identifies the assembly area refinement, and Note 5 identifies the level 5 assembly region. With reference to Fig. 10, refinement of the mesh around the turbine and wake region is achieved using the 'snappyHexMesh' utility. The wake region is defined as a cylinder 10 m radius, extending from the rotor to 50 m downstream, i.e. to the domain outflow (Fig. 10). The refinement level in this region is specified as level 2 i.e. subdivide the base cell/mesh twice in this region (Note 2 in Fig. 10). The rotor assembly and bladebox is set with a refinement level of 5, i.e. Note 5 in Fig. 10. The region around the rotor assembly (Note 3 in Fig. 10) is set at level 3 up to 10 m from the assembly.

Initial and boundary conditions
Inlet, outlet and wall boundary conditions are imposed as described in Section 5.1.2, with the addition of a no-slip boundary condition on the top wall of the channel. The kinetic energy initial condition and inlet is set to 0:0019m 2 s À2 , and the corresponding dissipation rate values are set to 0:000097m 2 s À3 .

Case runs and results
In the results section comparisons of the GAD-CFD model described in Section 3, and experimental results, are evaluated with respect to performance and wake characteristics. Each subsection reviews the results from each of the individual cases described in section 5. The results are presented with the respective input data for each study, and finishes with a discussion at the end of this section summarising key points from each case.

Case 1: mycek tidal turbine
TSR sweeps from TSR 1 to 7 and at 0:6m:s À1 and 1:2m:s À1 inflow velocities are conducted. The number of iterations to convergence ranges from approximately 300 to 1600, where the convergence criteria is set to 3 Â 10 À3 for all solved equations. Run times range from 10 min to 40 min for each TSR point, on a 10 core PC workstation.
Figs. 11 and 13 demonstrate the performance of the turbine in terms of C p (coefficient of power) and C t (coefficient of thrust) vs. TSR. It can be seen that the distribution of power and thrust correlate well with the experimental results reported in [16] for both 0:6m:s À1 and 1:2m:s À1 flow velocities. The spread of results for the differing flow speeds is due to the considerations of Reynolds numbers when sampling the lift and drag input data during the   Table 5, plotted against mesh number from coarsest to finest. The turbine is operating close to an optimal design Tip Speed Ratio (TSR) of 4.5 based on the local velocity at the rotor. Table 5 Location of samples used for the mesh study.   Table 7 Mesh grading in the x, y, and z directions using the studied configurations for Case Study 2.

Start
Middle End simulations. An over prediction can be seen in the stall region (TSR 2 to 4) of the C p curves.
The GAD-CFD model however only reports thrust acting directly on the rotor, thus a correction needs to be calculated considering the fluid drag acting on the assembly. This issue was examined in [17] and demonstrates good correlation for the combined results of thrust. With this in mind the thrust characteristics correlate well and show a similar trend to the observations of the experimental results.
The importance of modelling the correct assembly structure with respect to the downstream wake characteristics cannot be overstated. In Fig. 12 the velocity plot shows an asymmetric wake recovery in line with observations of the experimental results (i.e. Fig. 9 in [16]). This asymmetry can be observed more clearly in the turbulence intensity plot of Fig. 14. The inflow (left) shows a characteristic asymmetric sideways 'W', and in the centre of the plot there is an increase of turbulence intensity.
The asymmetry and distinct features are the result of the interaction of the assembly structure with the fluid medium and with the swirl imparted by the rotor. Remove the assembly structure and these features disappear (as shown in [17]). This phenomena can be observed more clearly with reference to Fig. 17. The visualisation shows fluid emanating from the assembly structure and propagating down stream in a helical fashion as a result of the wake rotation. The fluid can be seen mixing about the centreline of rotation approximately 5e6 diameters downstream. This fluid flow transports the turbulence emanating from the assembly structure, which pools before dispersing at about the same distance downstream. This pooling of turbulence can also be clearly seen in Fig. 14.
As part of this study a line sample of the velocity and turbulence data is taken behind the rotor to capture the fluid characteristics exiting the rotor. The line source is extended from the centreline of the device to twice the rotor radius, and is repeated downstream at 1.2, 4, 7, and 10 diameters. The results are presented in Fig. 15, which shows the distribution of flow, and Fig. 16 which demonstrates the distribution of turbulence intensity compared to the results from [16].
The results of the velocity and turbulence plots demonstrate good correlation with the experimental results, especially when considering the experimental deviation associated with this type of measurement, see [16] for more information. The distribution and magnitudes of the feature set show good correlation with the experimental results; this is particularly true of the asymmetric flow structures observed. The asymmetry is attributed to the Fig. 11. Comparison of predicted C P with experimental results [16] for a range of TSR's. Flow speed is 1:2m=s and 0:6m=s. Background turbulence intensity is 3%. Fig. 12. Numerical (top) and experimental (bottom) axial velocity slices at hub height. Free stream velocity is 0:8m:s À1 , TSR is 3.67, and background turbulence intensity is 3%. The top image is staged for direct comparison, while the bottom image is Fig. 9a from [16]. Image courtesy of [16] ©Renewable Energy.   Fig. 9c from [16]. Image courtesy of [16] ©Renewable Energy. assembly structure interacting with the rotating wake, see Fig. 17.

Case 2: NREL phase III wind turbine
In contrast to Case 1, lift and drag input data is taken from CFD results which are available at the required Reynolds numbers, aerofoil section thicknesses, and surface roughness (see [27] for more details). These are combined to match the foil design at given radial stations, allowing the generation of the required 3D lookup tables accessed as a function of Reynolds number (Re), radial station(r) and angle of attack (a) for both the lift and drag values, e.g.
Two TSR sweeps are conducted; one with the aerofoil pitch set to 5 , and one at 7 . All cases are run with an inflow velocity of 7:2m:s À1 . The number of iterations to convergence ranges from approximately 250 to 1600, where the convergence criteria is set to 3 Â 10 À3 for all solved equations. Run times range from 10 to 40 min per TSR point, on a 10 core PC workstation.
The results presented in Figs. 18 and 19 show C P and C T , respectively, plotted against TSR for the two different pitch angles reported. The trends of the C P curves demonstrate good correlation with the experimental results presented in [27]. There are no published experimental results to compare C T predictions with. However, the results follow sensible trends with lower coefficients of thrust for higher pitch angles. An over prediction in the C P results is observed in the stall region (approx. TSR 3 to 5) but this improves thereafter, with slightly better agreement for the lower pitch angle than the higher pitch angle. It should be noted that the experimental setup includes a cone angle of 3.4 which is not accounted for in the model, and it is likely that this will have some impact on the model performance predictions. Fig. 20 highlights the difference in wake structure of this wind turbine in contrast to Case 1. This is attributed to the different fluid medium, i.e. air instead of water, and thus the differing viscous properties. The fluid velocity is significantly higher, and is closer to the expected operating conditions of a HAWT.
In Fig. 20 the characteristic asymmetry highlighted in the previous case can be observed although the distribution within the wake differs. Both the velocity and turbulence fields are significantly asymmetric in the near wake region, with reducing asymmetry as mixing occurs further downstream.
The asymmetry is the result of the tower structure interacting with the rotating wake, as can be seen in Fig. 21. This visualisation Fig. 15. Numerical and experimental [16] axial velocity transects at hub height for 1.2, 4, 7 and 10 diameters downstream. Free stream velocity is 0:8m:s À1 , TSR is 3.67, and background turbulence intensity is 3%. The velocity [u*], and width [y*] are normalised by diameter. Fig. 16. Numerical and experimental [16] axial turbulence intensity at hub height for 1.2, 4, 7 and 10 diameters downstream. Free stream velocity is 0:8m:s À1 , TSR is 3.67, and background turbulence intensity is 3%. The velocity [u*], and width [y*] are normalised by diameter, and turbulence intensity [I%] is a percentage. Fig. 17. Visualisation of the simulation showing a turbulence intensity slice taken at hub centre, and a velocity iso surface at 0:32m:s À1 . The iso surface colour is mapped to turbulence intensity. Free stream velocity is 0:8m:s À1 , TSR is 3.67, and background turbulence intensity is set at 3%.
shows the effect of the tower interacting with the fluid medium and propagating downstream with the rotating wake. The wake structure for this case appears to recover further downstream than is available in the simulated domain. Fig. 3 shows how the revised downwash distribution system provides differing results depending on the foil geometry. This combined with the refinements discussed in Section 3.4 and illustrated in Fig. 4, improves the accuracy of the model for predicting power extraction and downstream wake characteristics.

General discussion
In Case 1(Mycek Turbine), the power and thrust (Figs. 11 and 13) show similar results to [17], likely due to the similarity between the analytically computed downwash distribution curve, Fig. 3, and the elliptical downwash distribution presented in [17]. Improvements can be clearly seen with respect to capturing the differing flow speeds, and better capture of the TSR trends influenced by changing foil sections and tip radius inclusion. However, over prediction in the stall region is observed.
Case 2 (Selig Turbine) demonstrates the application of the GAD-CFD model to a wind turbine. The results of the C P curves show good correlation although they over predict in the stall region. The C T results look sensible in the absence of any experimental results. Asymmetry within the wake structure is observed for the wind turbine as a result of the wake interacting with the tower structure. However, the wake distribution in general differs from the marine turbine mainly due to the differing fluid medium and Reynolds number operating range.
The additional refinements introduced, i.e. lift and drag as a function of Reynolds number, tip radius, radially changing foil sections, (see Sections 3.4 and 4), all contribute successfully to an improved representation of power and thrust through the entire range of TSR's studied for each case. The improvements in Fig. 18. Comparison of predicted C P with experimental results [27] for a range of TSR's. Flow speed is 7:2m:s À1 with blade pitch set to 5 , and 7:2m:s À1 with the blade pitch set to 7 . The background turbulence intensity is 3% at the inlet. Fig. 19. Prediction of C T for a range of TSR's. Flow speed is 7:2m:s À1 with blade pitch set to 5 , and 7:2m:s À1 with the blade pitch set to 7 . The background turbulence intensity is 3% at the inlet. Fig. 20. Visualisation of a numerical axial velocity slice (Top), and a numerical turbulence intensity slice (Bottom) taken at hub centre horizontally through the domain. Free stream velocity is 7:2m:s À1 , TSR is 5, and blade pitch angle is 7 . Fluid medium is air. Fig. 21. Visualisation of the simulation showing a turbulence intensity slice taken at hub centre, and a velocity iso surface at 4:05m:s À1 . The iso surface colour is mapped to turbulence intensity. Free stream velocity is 7:2m:s À1 , TSR is 5, and blade pitch angle is 7 . Fluid medium is air.
predicting rotor performance, combined with an accurate representation of the assembly structures, have a positive effect on the downstream wake characterisation. The run times vary between approximately 10 mine40 min for each TSR point. These timings depend on the size of the domain (quantity of cells), and the hardware capabilities (processing speed), and speed of convergence. This is a significant improvement over the use of fully resolved transient models. A rotating reference frame transient model of Case 1 with fully resolved hydrofoil geometry would require a mesh size of at least 15 million elements to maintain an appropriate y þ . This is in contrast to the less than 0.9 million elements used for Case 1. The reason for the additional cell count is due to the amount of refinement required at the hydrofoils to correctly capture the lift and drag characteristics. This combined with the time step size needed to maintain an acceptable Courant number (approximately 1 Â 10 À4 ), with a simulation time of approximately 20 min per time step, would require a minimum simulation time of approximately two days on the same hardware. Therefore, a complete TSR sweep would take approximately twenty days for this reference case on reference hardware.

Conclusions
This paper has described a new computational method for the accurate modelling of the interaction of renewable energy turbines with a fluid. As these devices will provide increasing amounts of global energy in future, this is an important topic. A balanced judgement should be made of the required length scales and accuracy compared to computational cost. The method sits between highly detailed blade resolved models and larger scale oceanographic and atmospheric models. It is important to not only be able to predict turbine performance at peak operating conditions, but throughout the entire TSR range. Additionally, thrust forces are important for structural design and the resulting downstream wake deficit and turbulence describes the interaction with other devices in the farm. Utilising well known loss coefficients can provide good turbine data but inaccurate wakes, and so the approach has a physically correct transfer of energy between the GAD and the fluid. To achieve this, all of the new features are required: a more concise downwash distribution computation, variation of foil section, application of tip radius correction, variation of lift/drag curves with Reynolds number and surface roughness.
The use of analytical methods to successfully and accurately predict the distribution of lift towards the tip of finite wing, i.e. accurately predicting 'roll off' losses for a given foil geometry and chord distribution, are demonstrated to produce excellent results for the accurate prediction of power and thrust.
Allowing for the variation of foil section shape within the model adds to the refinements including the distribution of forces along the foils. This helps produce better characteristics closer to the rotor hub, and also improved prediction in the stall region of the TSR range.
The variation of Reynolds number appropriate to the flow conditions improves the predictive capability both throughout the TSR range and with differing flow speeds. Although less variation is likely at the higher Reynolds numbers seen in full scale deployments, comparison with flume scale experimental studies show significant improvements in predictability.
Tip radius correction produces an effect that is more visible at higher TSRs; this in combination with including surface roughness effects for predicted lift/drag curves reduces over predictions in the overspeed part of the TSR range.
Experimental comparison shows that the downstream wake characteristics are better represented with the revised algorithm, particularly with the inclusion of an accurate representation of the assembly structure. Applied to a selection of cases, this study provides confidence that this approach can be applied to a range of scenarios; both small scale flume studies, and large scale deployments in both the marine and wind environments. The GAD-CFD model by design can take variable inflow conditions thus the prediction of downstream turbines with varying inflow angles and velocity profiles has become computationally viable. Therefore, it is appropriate to consider this model to study turbine arrays and their interaction with respect to local topography and power control. Further possible improvements to the model are related to additional geometric aspects, i.e. capturing a blade coning angle, and by extension blade flexing either by design or environment. This could be coupled with a structural model for design studies. Changing lift and drag characteristics at higher levels of free stream turbulence is an area of further study. As the turbulence levels increase, the quantity of lift and drag changes, as does the stall point relative to angle of attack, and post stall features significantly change. Finally, the model should be combined with turbine control algorithms that consider power capping through stall or pitch control to enable the study of rotors interacting in an array. It is the authors intention to publish the implementation of the model to enable such studies to take place.