Numerical Investigation of Aerodynamic Performance and Loads of a Novel Dual Rotor Wind Turbine

Abstract: The objective of this paper is to numerically investigate the effects of the atmospheric boundary layer on the aerodynamic performance and loads of a novel dual-rotor wind turbine (DRWT). Large eddy simulations are carried out with the turbines operating in the atmospheric boundary layer (ABL) and in a uniform inflow. Two stability conditions corresponding to neutral and slightly stable atmospheres are investigated. The turbines are modeled using the actuator line method where the rotor blades are modeled as body forces. Comparisons are drawn between the DRWT and a comparable conventional single-rotor wind turbine (SRWT) to assess changes in aerodynamic efficiency and loads, as well as wake mixing and momentum and kinetic energy entrainment into the turbine wake layer. The results show that the DRWT improves isolated turbine aerodynamic performance by about 5%–6%. The DRWT also enhances turbulent axial momentum entrainment by about 3.3%. The highest entrainment is observed in the neutral stability case when the turbulence in the ABL is moderately high. Aerodynamic loads for the DRWT, measured as out-of-plane blade root bending moment, are marginally reduced. Spectral analyses of ABL cases show peaks in unsteady loads at the rotor passing frequency and its harmonics for both rotors of the DRWT.


Introduction
Modern utility-scale horizontal axis wind turbine (HAWT) rotor blades are aerodynamically optimized in the outboard region, whereas the blade root region is designed primarily to withstand structural loads.Therefore, very high thickness-to-chord ratio airfoils, which are aerodynamically poor, are used in the blade root region to provide structural integrity.Up to 5% loss in wind energy extraction capability is estimated to occur per turbine due to this compromise [1].This "root loss" occurs even in turbines that operate in isolation, i.e., with no other turbine nearby.Most utility-scale turbines are deployed in clusters, with multiple turbines operating in proximity of each other.Array interference (wake) losses resulting from aerodynamic interaction between turbines in wind farms have been measured to range between 8% and 40% depending on wind farm location, farm layout/wind direction and atmospheric stability condition [2].
Flatback airfoils [3] and vortex generators [4] have been used in the blade root region to mitigate blade root losses.Improvements in wind farm efficiency have been sought by optimizing wind farm layout so as to minimize wake interference between turbines [5][6][7].Ideas for wind farm efficiency improvement include the use of counter rotating vertical axis wind turbines (VAWTs) for large-area wind farms to increase power production per unit land area [8].Other ideas have been pursued for horizontal axis wind turbines in wind farms with set turbine layout.Redirecting the wakes of upstream turbines through yaw misalignment [9,10] is one such concept.By yawing an upstream turbine, the force exerted by the turbine on the flow (reaction to the thrust force) is turned slightly in the cross-flow direction.A component of this force then acts in a direction perpendicular to the flow velocity, serving as a centripetal force to curve the mean flow and divert the flow/wake away from the turbines immediately downstream.Another idea [11] aims to reduce wake loss through manipulation of the turbulence in the turbine wake by changing the induction factor for the turbine rotor.This can be achieved by various means, such as altering the pitch of the blades, the RPM of the rotor or the yaw of the nacelle.
Wind turbines with two rotors (typically arranged in tandem so that the incoming flow stream area is unchanged) are known as dual rotor wind turbines (DRWTs).Newman [12] developed a multi-rotor actuator disk theory and demonstrated that a turbine with two equally-sized rotors could capture up to 8% more energy than a corresponding single-rotor wind turbine (SRWT).Previous research on DRWTs has been focused on increasing wind energy capture by harvesting the kinetic energy left in the wind after it passes through the turbine rotor.Jung et al. [13] explored a 30-kW counter-rotating dual-rotor wind turbine.It featured an eleven-meter diameter main rotor with a 5.5-m auxiliary rotor located upwind of the main rotor.This DRWT uses a bevel gear to couple the counter-rotating shafts.The authors used quasi-steady strip theory and a wake model to predict the performance of several DRWT configurations.They predicted a 9% increase in the turbine power coefficient (C P ) when compared to a single-rotor configuration.Other studies have led to patents including Kanemoto and Galal [14,15] who propose a DRWT with two differently-sized upwind rotors driving a generator with two rotating armatures.
The DRWT technology by Rosenberg et al. [16,17] (see Figure 1) takes a different approach; it aims at reducing losses (blade root and wake losses) in wind turbines and wind farms.It utilizes a secondary, smaller, co-axial rotor to mitigate blade root losses and to enhance mixing of the turbine wake.Rosenberg et al. [16] and Selvaraj [18] introduced this turbine technology and presented preliminary aerodynamic analyses of a DRWT design.The conceptual 5-MW offshore turbine by the National Renewable Energy Laboratory (NREL) [19] was used as the baseline single-rotor design and also as the main rotor for the DRWT.The secondary rotor was designed using an inverse design approach based on the blade element momentum theory.The design and optimization approach used Reynolds averaged Navier-Stokes (RANS) computational fluid dynamics (CFD) simulations with an actuator disk representation [20] of turbine rotors.RANS CFD analyses showed an increase in C P of around 7% with the DRWT.
In this paper, we extend the analyses of [16][17][18] by including the effect of the atmospheric boundary layer and investigate turbine aerodynamic performance and loads, as well as wake mixing.We present comparative (between DRWT and SRWT) isolated turbine aerodynamic analyses for uniform inflow with no incoming turbulence and two atmospheric stability conditions: neutral and stable.An improvement of about 6% in C P through root loss reduction is demonstrated.The analysis of turbine wake shows increased momentum and kinetic energy entrainment in the wake layer of the DRWT.One of the concerns with the DRWT technology is the potential of increased unsteady loads on the primary rotor due to its proximity with the secondary rotor.These loads are computed numerically using large eddy simulations (LES) and reported as power spectral densities of out-of-plane blade root moment.No significant increase in loads is observed for the DRWT.
The remainder of the paper is laid out as follows.A summary of the numerical method utilized in this study and its validation against experimental data are presented first.Code validation results are also presented in this section.Section 3 summarizes the computational setup, grids and simulation details, including the assumptions made in the present calculations.Aerodynamic performance results are described in Section 4, wherein comparisons are drawn between an SRWT and a DRWT operating in uniform inflow and in neutral and stable atmospheric boundary layer (ABL) flow.Section 4 also investigates the aerodynamic loads on the two rotors of the DRWT for the different inflow conditions.The conclusions are presented in the final section.

Figure 1.
A cartoon of the dual-rotor wind turbine (DRWT) technology proposed by Rosenberg et al. [16].

Numerical Method
A wide range of methods can be used to model wind turbine and wind farm aerodynamics.Analytical models [21,22] and semi-analytical models, such as blade element momentum theory, vortex lattice methods, etc. [23], have been extensively used for the design and analysis of wind turbines operating in isolation.
Wind turbine wake dynamics and wind farm aerodynamics have been investigated using parabolic methods [24] and computational fluid dynamics (CFD) methods.CFD methods can range from time-averaged Reynolds averaged Navier-Stokes (RANS) simulations [25] to large eddy simulations (LES) [26] that resolve energy containing turbulence in the atmosphere and turbine wakes.Vermeer [27] provides an overview of wind turbine aerodynamics, as well as wind farm aerodynamics through a survey of existing numerical, as well as experimental work.Sanderse [28] reviews different numerical methods currently being used for aerodynamic analysis of wind turbines and wind farms.
Recent research on numerical modeling of wind turbine and wind farm aerodynamics has largely focused on using LES (see, e.g., [29][30][31][32][33][34][35]).Jimenez et al. [29,30] used incompressible LES to model the aerodynamics of a single wind turbine (modeled as an actuator disk) in an atmospheric boundary layer.Troldborg et al. [33] investigated the aerodynamic interaction between two turbines using an actuator line model coupled with an LES flow solver.Aerodynamic interaction between the turbines was simulated for varying atmospheric turbulence intensity, distance between the turbines and partial and full-wake operation of the downstream turbine.Porté-Agel et al. [34] investigated wake losses in an offshore wind farm with varying wind direction using LES.Stevens et al. [35] investigated the effects of the alignment of turbines in a wind farm and identified optimal staggering angles to use for micrositing (wind farm layout).
The Simulator for Wind Farm Application (SOWFA) [36,37] software is the LES flow solver used in this work.SOWFA solves spatially-filtered, incompressible forms of continuity and Navier-Stokes equations (see Equation ( 1)).The grid-filter width, computed as the cube-root of the cell volume ∆ = (∆ x ∆ y ∆ z ) 1/3 , is used as the spatial filter width.Unresolved, sub-filter (or subgrid) scale stresses introduced by spatial filtering are modeled using a subgrid model.Turbine rotor blades are parameterized using the actuator line model (ALM); blade geometry is not resolved.The actuator line model uses lookup tables for airfoil polars to compute sectional lift and drag forces and applies them as body forces.The governing equations are written in spatially-filtered quantities (denoted by overhead (˜)) as: In the above equation, θ is the potential temperature, α is the thermal diffusivity of the fluid and f i is the the force exerted by turbine rotor blades computed using lookup tables for airfoil lift and drag polars; p * = p/ρ 0 + u j u j /2 is modified kinematic pressure; τ ij = u i u j − u i u j is the subgrid scale (SGS) stress tensor; q j = u j θ − u j θ is the SGS heat flux; ρ 0 and θ 0 are constants based on the Boussinesq approximation.For simulating the atmospheric boundary layer, the flow is driven by a pressure gradient, δ i1 F P ; the coordinate system is such that index "1" corresponds to the streamwise direction, "3" points up and normal to the ground and "2" is determined by the right-hand rule.The DRWT is modeled in SOWFA by simulating the two rotors of the DRWT as two single-rotor turbines operating in tandem.
The deviatoric part of the SGS stress tensor (τ ij ) is modeled using an eddy-viscosity model, and the SGS heat flux with an eddy-diffusivity model: q j = u j θ − u j θ = −(ν sgs /Pr sgs )∂ θ/∂x j , where S ij = 1/2 ∂ u i /∂x j + ∂ u j /∂x i is the resolved strain-rate tensor and Pr sgs is the SGS Prandtl number.The mixing length model by Smagorinsky [38] is used to model eddy viscosity as ν sgs = (C S ∆) 2 | S|.In the original model, C S was assumed to be a constant, but dynamic calculation of this coefficient has been used in recent years [39,40].Improved, tuning-free, scale-dependent SGS models have also been developed (see e.g., [31]) and used for atmospheric flow and wind farm simulations.The standard Smagorinsky model with C S = 0.135 is used here.SOWFA uses a finite volume formulation, and the discretization is second order accurate in space (central) and time (backward).Details about the SOWFA software can be found in [41].A two-step solution procedure is used.In the first (precursor) step, the turbines are removed, and turbulent flow in the domain (the ABL) is simulated using periodic boundary conditions in the streamwise and cross-stream directions; the flow is driven by an adjustable pressure gradient.After the solution reaches a statistically-stationary state, time-accurate data are sampled at every time step on the inlet plane(s) of the computational domain and stored.These data are specified as a boundary condition for the subsequent wind farm calculations.
Viscous effects are negligible everywhere except near surfaces (ground, in the present case) due to the high Reynolds number in ABL flows.Energy containing eddies near the ground can become very small, and resolving such small scales can lead to exorbitant grid sizes.Surface flux models for stress and heat (e.g., Moeng [42]) are therefore usually used in wind farm computations.Moeng's models require as input the surface roughness height, h 0 , the horizontally-averaged surface heat flux, q s , and a measure of the horizontally-averaged shear stress specified as friction velocity, u * .While h 0 and q s are directly specified (from estimates in the literature for different surfaces, sea, grasslands, forest, etc.), u * is approximated using the Monin-Obhukhov similarity theory [43].
Lee et al. [44] coupled the LES solver for wind farm aerodynamics and SOWFA with the structural dynamics solver in the FAST (Fatigue, Aerodynamics, Structures and Turbulence) code [45] to enable calculation of fatigue loads due to atmospheric and wake turbulence.Through this coupling, the simplified aerodynamics module (including the turbulent inflow model) in FAST is replaced by LES, which provides much higher fidelity in resolving the flow.Since the focus of this paper is limited to aerodynamic performance and loads, and not aeroelasiticity, the turbine is idealized by assuming the rotor blades to be infinitely stiff.Structural dynamics can have a significant effect on blade/turbine/tower loads, and hence, the present loads analysis is only a preliminary investigation.In order to avoid confounding effects from the controller, the turbine is assumed to operate at a fixed, user-specified rotation rate.The turbine RPM and incoming mean wind speed are set to achieve the design tip speed ratio.The instantaneous tip speed ratio fluctuates with the incoming turbulent wind.

Validation
The SOWFA solver is first validated against experimental data for the three-bladed, stall-controlled, 100-kW Tellus turbine (measurement data from [46]).This turbine is referred to as the Risø turbine in [46].The turbine rotor diameter is 19 m, and its blade chord and twist distributions are shown in Figure 2a. Figure 2b compares LES predictions of the power curve against measured data, as well as against blade element momentum (BEM) theory predictions.Stall-controlled turbines run at a fixed rotational velocity.Above a certain wind speed corresponding to rated power, they begin to stall to reduce power generation.This turbine stalls for wind speeds greater than 12 m/s.It is well known that spanwise flow alleviates stall in 3D blades and allows the blade to operate at higher angles of attack than a 2D blade would.This alleviation of stall cannot be simulated with the actuator line method.Stall-corrected 2D polars can be used to partly address this weakness of the model, but it was not pursued here, as the focus is on evaluating DRWT performance at the design condition.Therefore, the comparison in Figure 2b is limited to pre-stall operation.Comparisons of the radial distributions of sectional torque (tangential) force coefficient, c τ F , between the predictions made by LES and BEM theory methods for two rotor tip speed ratios, λ = Ω r tip /u ∞ = 6.0 and 7.0 (close to the design λ), are presented in Figure 3. c τ F is the non-dimensional aerodynamic force on a blade section that generates torque and is defined as: where F τ is the component of the net aerodynamic force (per unit length) in the plane of rotor rotation, ρ is the fluid density, u rel is blade relative flow velocity, c is section blade chord, c l and c d are section lift and drag coefficients, respectively, and φ is the angle that the blade relative velocity vector makes with the plane of rotor rotation.The agreement between the two solvers is very good except very near the point where the blade cross-section abruptly transitions from an airfoil shape to a cylinder (around 0.3× the tip radius).The differences between the two models emanate from the different blade discretization and interpolation used to calculate sectional airfoil properties (chord, twist and polars) and are magnified at the geometric transition point.

Computational Setup
A two-step approach is used for the numerical predictions.The atmospheric boundary layer (ABL) flow is computed first on one grid, and then the flow around the wind turbine is computed on another, more refined grid.Investigations are conducted for uniform inflow and two atmospheric stability conditions: (1) neutral; and (2) stable.In the simulations presented here, the wind blows from the southwest; the wind direction is 240 degrees (measured clockwise) from the north, such that the west and the south boundaries of the computational domains are the inlet, while the east and the north boundaries are the outlet (see Figure 4).
Three coordinate systems are utilized in this study.A frame of reference attached to the ground, described by unit vectors ê x, ê ỹ, êz , is utilized for the CFD calculations.The freestream wind blows at an angle φ = 30 o to ê x.A new coordinate system with its x-axis aligned with the flow direction is therefore defined by the unit vectors êx , êy , êz , where ê x. êx = cos(30 o ).Lastly, a cylindrical coordinate system, with its axis aligned with êx , is used to compute momentum and energy entrainment into the turbine wake layer (see Section 4.2.3).This coordinate system is defined by the unit vectors êr , êθ , êx .The details are provided in Appendix A.  The atmospheric boundary layer (ABL) is developed in a computational domain by performing "precursor" simulations.A precursor calculation simulates an infinitely long domain in the horizontal directions (the Earth's surface) by using cyclic boundary conditions.The intent here is to create statistically steady ABLs under different stability conditions and not to capture the transient effect caused by the diurnal or seasonal variation of the ABL.Wall models by Moeng [42] are used to estimate the surface stresses (viscous and SGS) and temperature flux at the bottom boundary.The aerodynamic surface roughness is h 0 = 0.1 m, which is a typical value for a terrain with low crops and occasional large obstacles [47].The top boundary is at 1 km, which is many diameters away from the turbine.The velocity normal to the boundary is set to zero.The pressure boundary condition is zero-gradient, and the temperature gradient is specified to be 0.003 K/m at the top boundary.Temperature inversion is applied in the domain with the width of 100 m.The temperature at the bottom of the inversion layer is 300 K and at its top is 308 K. Above the inversion layer, the temperature gradient is 0.003 K/m.The inlet is on the south and the west boundaries and a zero pressure gradient boundary condition to simulate outlet conditions is applied on the north and the east boundaries.
Since the objective of the precursor simulations is to establish a fully-developed atmospheric boundary layer flow, wind turbines are excluded from these simulations.The flow is driven by a pressure gradient, which is adjusted to achieve the desired flow speed at the turbine hub height.Random perturbations are applied to the flow initially to trip the boundary layer.Precursor simulations are carried out for a long enough time (5 h of simulation time) to achieve statistical stationarity.With the mean flow speed at the turbine hub height set at 8 m/s, it takes about nine minutes of simulation time for the flow to cross the computational domain from the inlet to the exit.In 5 h of simulation time, the flow cycles the domain approximately 34 times, which is sufficient to achieve a fully-developed ABL.
Time-accurate flow data sampling is initiated at the inlet boundaries after reaching a statistical stationary state (at t = 5.0 h).Sampling is performed for 1000 s of simulation time, which is the total simulation time for the wind turbine calculations.These time-accurate flow data are prescribed as an inflow boundary condition in the wind turbine simulation.This is one approach to prescribe time-accurate inflow boundary conditions.Another approach is to use synthesized turbulence, where the time-accurate flow information can be constructed using analytical models of turbulence spectra (e.g., Kaimal, von Karman, etc.) [48].Once statistical stationarity is reached (at t = 5.0 h), the entire flow field from the precursor simulation is also stored and used to initialize the flow in the wind turbine simulations.
Each wind turbine simulation starts at t = 5.0 h, and a total of 1000 s of flow is simulated.The initial and inlet boundary conditions are prescribed using the precursor simulation data as described above.The computational domains for precursor and wind turbine simulations need not be identical.On the contrary, it is desirable to make the domain of the precursor calculation considerably larger than the domain for the wind turbine simulation to account for the length scale disparity between the physical phenomena of interest in these simulations.In the precursor simulations, the energy containing turbulent eddy size can be of the order of a kilometer (planetary boundary layer height), while in the wind turbine simulation, the flow physics of interest is turbine wake for which the relevant length scale is of the order of the turbine diameter (∼ 100 m for utility-scale turbines).
A shorter domain is used for wind turbine simulations in comparison with precursor simulations to allow for higher spatial resolution of the turbine wakes.As seen in Figure 4, the precursor runs are performed on a computational domain with dimensions 3 km × 3 km × 1 km.The domain is discretized into 288 × 288 × 100 hexahedral cells.The precursor simulations are performed in parallel on 128 cores, and each run takes about 50 h of wall-clock time.The wind turbine runs are performed on a domain of size 2.2 km × 1.5 km × 1 km that is discretized into 220 × 150 × 100 hexahedral cells.Two levels of mesh refinement are used in the vicinity of and downstream of the wind turbine.Each grid cell in a refinement block is split into half along each direction (i.e., it is divided into eight cells).In the refined block, the cell size is 2.5 m in each direction.The final mesh has a total of about 14 M cells.The blockMesh utility in OpenFoam is used to generate the computational meshes.It takes approximately 120 h of wall-clock time to simulate 1000 s of flow in each wind turbine simulation using 128 cores in parallel.Figure 4b shows the turbine location in the computational domain, as well as the topology of the refinement zones.
LES simulations are carried out using the optimum DRWT identified in Rosenberg et al. [16].The non-dimensional chord and twist distributions of the main and secondary rotors of the DRWT are shown in Figure 5.The blade chord (c) and radius (r) are nondimensionalized by the respective rotor tip radii.To enable direct comparisons, simulations are also performed for a conventional single rotor wind turbine (SRWT), which is the conceptual NREL 5-MW offshore reference turbine [19].The NREL 5-MW turbine rotor is used as the primary rotor of the DRWT.

Results and Discussion
The objectives of the DRWT technology by Rosenberg et al. [16] are two-fold: (1) minimize root losses; and (2) increase entrainment from the upper atmosphere into the turbine wake.An increase in C P between 5% and 7% through root loss mitigation was already demonstrated for uniform inflow calculations in Rosenberg et al. [16].Here, we analyze the effect of the atmospheric boundary layer (ABL) on blade root loss, wake mixing and aerodynamic loads.

Atmospheric Boundary Layer
The precursor (ABL) simulations are analyzed first.The simulated flow data are averaged in time and space (along the horizontal directions) to obtain mean velocity profiles for both the neutral and stable ABL conditions.A pressure gradient along the flow direction is imposed and is continually adjusted until the desired mean flow speed (u ∞,h = 8 m/s for the simulations here) is achieved at the turbine hub height.The atmospheric stability is varied between different simulations by changing the heat flux through the bottom boundary.Zero net heat flux is prescribed for the neutral ABL simulation, while a value of −0.4 K•m/s is used to simulate the stable ABL condition.A negative value implies heat flux out of the computational domain.
Figure 6a plots the normalized mean velocity profiles for two ABL simulations; the freestream mean wind speed at the turbine hub height ( ū∞,h ) is used for the normalization.The uniform flow case does not require a precursor simulation; a uniform wind speed of 8 m/s with zero inflow turbulence is used for that case.As expected, the vertical shear in the mean wind speed near the ground in the stable case is much higher than in the neutral case.Figure 6b compares the streamwise turbulence intensity, defined as 2 / ū∞,h , where u x is the streamwise component of the wind velocity.Equation ( 4) is used to calculate streamwise turbulence intensity from the computed Reynolds stress tensor.In both ABL cases, the streamwise turbulence intensity decreases with height across the turbine rotor.At the hub height, the streamwise turbulence intensity (σ u x /u ∞,h ) is about 6% for the stable ABL case and about 10% for the neutral ABL case.Stable stratification suppresses atmospheric turbulence, and hence, the streamwise turbulence intensity is less for the stable ABL case compared to the neutral case.Atmospheric boundary layer stability is characterized by the Richardson number, which can be written as: where z is height from the ground, θ is potential temperature, u is flow speed, L is the Obukhov Length and φ m and φ h are nondimensional temperature and wind profiles, respectively.For z/L ≥ 0, φ h = 0.74 + 4.7z/L and φ m = 1 + 4.7z/L.z is set as 50 m for the ABL calculation of the Richardson number [49,50].For neutral stability, L → ∞ and Ri → 0. For the stable condition simulated, Ri = 0.125, which falls in the "slightly stable" category according to the classification of atmospheric stability by Sedefian and Bennett [49].They collected meteorological data over a year at a site in Staten Island, NY, USA.They categorized all of the possible ABL stability conditions into seven groups, ranging from strongly unstable to strongly stable ABL.As they reported, neutral and slightly stable atmospheres (as characterized by Ri using Equation (3)) were observed for 31.1% and 12.3% of the time, respectively.Due to their high probability of occurrence in the field, these two stability conditions were selected for the simulations.

Aerodynamic Performance
The aerodynamic performance of the DRWT design is analyzed and compared against the corresponding SRWT.Analysis is presented for the three inflow conditions described in Section 3: (1) uniform inflow; (2) slightly stable ABL; and (3) neutral ABL.The emphasis is on evaluating the ability of the DRWT to mitigate blade root loss and wake loss.
Figure 7 visualizes the flow field downstream of the SRWT turbine operating in a neutral ABL.The figure is an iso-surface plot, showing surfaces with a constant value of Q-criterion.In a vortex, Q = 0.5(|Ω| 2 − |S| 2 ) is greater than zero, where Ω and S are vorticity and rate-of-strain tensors, respectively, and |T| denotes the Euclidean norm of the tensor T. The tip and hub vortices can be seen in the near-wake region, but quickly disintegrate in one to two rotor revolutions due to the high atmospheric turbulence.Some structures can be seen underneath the turbine, which correspond to the ABL turbulence near the ground.

Root Loss Mitigation
Airfoils in the root region (approximately inner 25%) of conventional turbine blades have very high thickness-to-chord ratios and, thus, have poor aerodynamic performance.The smaller secondary rotor in the DRWT proposed by Rosenberg et al. [16] aims to mitigate the losses due to the aerodynamically non-optimum root region of the larger primary rotor blades.Since the secondary rotor is much smaller in size compared to the primary rotor, it has significantly lower loads; loads increase in proportion with the cube of rotor diameter.The secondary rotor can therefore be designed using relatively thin, aerodynamically-optimized airfoils that are efficient at extracting wind energy passing through the main rotor blade root region.
The C P of the DRWT is computed as the sum of the powers generated by the two rotors, normalized by the power in the wind going through the main rotor disk area.RANS calculations by Rosenberg et al. [16] showed an increase in C P of about 7% with the DRWT.The results from the LES simulations for the three inflow conditions, conducted as a part of this study, corroborate the findings of Rosenberg et al. [16].The DRWT outperforms the conventional SRWT for all inflow conditions in terms of time-averaged power produced, demonstrating the "root loss mitigation" ability of the DRWT.There is also a marginal reduction in the time-averaged out-of-plane blade root bending moment (M OOP ) of the primary rotor of the DRWT.A summary of these results is presented in Table 1.The percentage differences in the table are calculated as 100×(DRWT-SRWT)/SRWT.The largest increase in power is observed for the neutral ABL condition and smallest for the uniform inflow case.Comparative wake mixing analyses are performed for the DRWT and SRWT to investigate the ability of the DRWT to reduce wake loss.The setup and inflow boundary condition for the neutral and stable cases are described in Section 3.For the uniform inflow simulation, the precursor run is not required, and a uniform wind velocity and zero turbulence intensity are specified at the inlet boundaries.The flow angle is kept the same (30 degrees w.r.t.ê x) as for the two ABL cases.

Mean Flow
Time-averaged velocities in turbine wakes are compared first.Averaging is performed over 110 revolutions of the main rotor.Large scale eddies in atmospheric turbulence contain significant energy at frequencies much lower than rotor passing frequency; hence, averaging over a long time (in comparison with the rotor rotation period) is required.Contours of mean streamwise velocity are plotted on an x-z plane in the streamwise direction (see Subplots "a" and "b" in Figures 8-10) and on y − z planes in the cross-stream direction at four downstream locations: 2D, 4D, 6D and 8D, where D is the main rotor diameter (see Subplots "c" and "d" in Figures 8-10).The cross-stream contour plots are drawn over disks of a diameter 1.4-times the turbine rotor diameter, such that the top of each disk corresponds to the 12 o'clock position of the rotor and the bottom to the six o'clock position; the view is looking from downstream (but at an angle, so the disks look oval and not circular), and the main turbine rotor is rotating counter-clockwise in the figures.The path traced by the rotor tip is marked with the dashed curves.For a more quantitative comparison, streamwise velocity is azimuthally (circumferentially) averaged; averaging denoted by angle brackets as ūx .Radial profiles of ūx , normalized by ū∞,h , are compared between the SRWT and the DRWT in the wake region (see, e.g., Figure 8e).
Figure 8 compares the SRWT and DRWT designs for the uniform inflow case.Large differences in mean velocities near the center of the disks are evident up to 8D downstream of the rotor.The SRWT is not effective at extracting energy from the wind in the blade root region, and hence, the wind flows through there unharvested.In contrast, the secondary rotor of the DRWT efficiently extracts the energy from this streamtube and, hence, leaves a larger momentum (velocity) deficit in the wake.Reduced values of ūx /u ∞,h are therefore observed in the wake of the DRWT as compared to SRWT for r/r tip < 0.4.This deficit reduces with downstream distance due to enhanced wake mixing in the DRWT.The circumferentially averaged plots show that by 8D downstream, the difference in the velocity (and hence, kinetic energy) deficit between the DRWT and SRWT is considerably reduced.Also noted in the uniform inflow case is the slight flow acceleration near the ground under the turbines, which occurs due to the contraction of the streamtube caused by the slip wall boundary condition imposed on the ground.As a result, the wake deficit shifts up and away from the ground as it convects downstream (see Figure 8).
The addition of the atmospheric boundary layer (ABL) considerably changes the turbine performance and wake dynamics.The wall shear reduces the velocity near the ground.Hence, the acceleration effect observed in the uniform inflow case, due to streamtube contraction, is not observed in the ABL cases.The effect of wake rotation is evident in the ABL cases as seen by the azimuthal locations of highest mean velocity in Subplots (c) and (d) of Figures 9 and 10.The main rotors of the turbines are rotating in the counter-clockwise direction as seen from downstream, and the turbine wake rotates in the opposite (clockwise) direction.As the wake rotates, it pulls the higher momentum fluid from above (12 o'clock position), down and to the right in the figure (clockwise).Wake rotation also pulls the low momentum fluid in the boundary layer near the ground, up and to the left in the figures.The primary mechanism of wake mixing is turbulent momentum transport across the turbine wake layer, which is proportional to the turbulence level in the incoming wind.Higher turbulence intensity in incoming flow therefore leads to a higher wake mixing rate.The turbulence in the ABL enhances wake mixing, and hence, the velocity deficits in the wakes are reduced for the two ABL cases in comparison with the uniform inflow case.The wake deficit for either turbine is highest for the uniform flow case and smallest for the neutral ABL case.This behavior of the turbine wake mixing rate increasing with inflow turbulence intensity has been previously observed both experimentally [51] and numerically [52] for conventional, single-rotor turbines.
The interesting observation in the present study is that the difference in the velocity deficits between the SRWT and DRWT also reduces faster for the ABL cases when compared to the uniform flow case; the difference is negligible by 8D downstream (compare Subplot "e" in Figures 9 and 10 with Figure 8).This suggests that the presence of atmospheric turbulence enhances the wake mixing rate of the DRWT more than that of the SRWT.The increase in wake mixing rate is due to the ABL turbulence augmenting the interaction between the trailing wake/vortex systems of the two rotors of the DRWT; this interaction leads to increased turbulent momentum transport into the turbine wake layer, as shown later in Section 4.2.3.The largest increase in wake mixing rate is observed for the neutral ABL case (Figure 10).
It should be noted that the secondary rotor of the DRWT used here is not specifically designed or operated to target wake mixing.The results obtained here, of enhanced wake mixing rate with a DRWT, indicate that the technology has potential to reduce wake losses if its design and operation are optimized for that purpose.Furthermore, it is observed that even with the de facto design/operation of the DRWT, the additional wake deficit due to the secondary rotor of the DRWT is sufficiently mixed out such that it will not adversely impact the performance of the immediately downstream turbine if it is placed at least 8D away.Therefore, the gains achieved by mitigating root losses for isolated turbines will be realized in array configurations (wind farms), as well.

Turbulence Intensity
Mechanical turbulence, generated due to mean velocity shear in turbine wake and tip vortices, is an important contributor to wake mixing.Streamwise turbulence intensities in the wakes are therefore compared between the DRWT and SRWT configurations for the three inflow conditions in Figures 11-13.In all three cases, streamwise turbulence intensity near the turbine hub height is higher for the SRWT until about 4D downstream of the turbine.This is because of two reasons: (1) due to the "leakage" of high-momentum fluid through the blade root region, a higher velocity shear is present in the wake of the SRWT than that of the DRWT; compare, e.g., Subplots (a) and (b) in Figure 8. Mean velocity shear is a source of turbulence; hence, there is higher shear-generated turbulence in the wake of the SRWT near the turbine hub height up to 4D downstream of the turbine; and (2) the secondary rotor of the DRWT dampens the large scale eddies in the incoming ABL flow.This phenomenon of a turbine rotor acting as a damper/high-pass filter has been reported by Chamorro [53] for conventional single rotor wind turbines.
While the SRWT shows higher turbulence intensity in the wake up to about 4D, the turbulence intensity for the DRWT is higher further downstream for all inflow conditions.This increase in turbulence production in the DRWT for x > 4D is due to the interaction between the trailing wake/vortex systems of the primary and the secondary rotors of the DRWT.The turbulence production mechanism is velocity shear, which is increased by this interaction.Higher turbulence intensity in the wake explains the higher wake mixing rate for the DRWT noted in the previous section.The difference in turbulence intensities between the DRWT and the SRWT is highest for the neutral ABL case (which has the highest ABL turbulence of all of the cases simulated), re-emphasizing that atmospheric turbulence promotes the interaction of the wake/vortex systems of the two rotors.
While increased turbulence intensity enhances the wake mixing rate and therefore boosts wind farm efficiency, it does have a potential to increase fatigue loading on downstream turbines.This potential issue with the DRWT technology will be addressed in the future.Aerodynamic loads for a DRWT operating in "clean" flow, however, are analyzed in Section 4.3.

Momentum Entrainment
Utility-scale turbines in wind farms are typically installed in systematic arrangements (arrays).Wake losses in a wind farm are most severe when the wind direction is aligned with the turbine rows.In such extreme cases, it is meaningful to look into the entrainment of high-momentum fluid into the turbine wake layer.For a very large, closely-packed turbine array, the turbine wake layer would be horizontal, stretching vertically from H − r tip to H + r tip , where H is the turbine hub height and r tip is the rotor tip radius.For an isolated turbine analysis, however, the turbine wake layer is cylindrical, co-axial with the rotor, and has the same diameter as that of the turbine rotor (see Figure 14).
Turbulent transport of momentum from outside the turbine layer has been identified to be the dominant mechanism for re-energizing the flow in large wind turbine arrays [54].Turbulent momentum flux through the cylindrical turbine wake layer (see Figure 14) is analyzed to compute turbulent momentum and energy transport and investigate wake recharging.The present analysis is performed for a turbine wake layer that extends from the turbine rotor location to 8D downstream.Time averaged velocity and the Reynolds stress tensor (u i u j ) are interpolated onto this cylindrical surface using the Tecplot 360 software (www.tecplot.com).
Turbine wake layer Entrainment of high momentum fluid into this cylinder is induced by turbulent stresses, particularly the stress term u r u x , where the subscript "r" denotes the radial component, the prime denotes a perturbation quantity and the overline denotes a time-averaged quantity.A cylindrical coordinate system ( êr , êθ , êx ) with its axis aligned with the freestream flow direction êx is used (see Appendix A).Equation ( 6) relates u r u x to the Reynolds stress tensor computed in the ground frame of reference ( ê x, ê ỹ, êz ).Since êr points radially outward, a negative value of u r u x implies transport into the cylinder (turbine wake layer).
Figure 15a-f plots the Reynolds stress term u r u x normalized by u 2 ∞,h on the unwrapped cylindrical surface that defines the turbine layer.The color map in the figure is reversed to make intuitive sense: negative values indicate net flux into the cylinder.In the figure, θ = 0 • corresponds to the 12 o'clock position and ±180 • correspond to the six o'clock position.It can be noticed that the highest value of −u r u x is not at the 12 o'clock position, but skewed about it.This behavior has been noted previously by Wu and Porté Agel [55], although for the stress term u x u z .A clear trend of increasing momentum entrainment with increasing ABL turbulence is seen in the figure for both SRWT and DRWT.The x-location where peak entrainment occurs varies both with atmospheric turbulence intensity and azimuthal location.The peak values occur closer to the turbine as inflow turbulence intensity increases (compare Subplot "a" with "b" and "c" of Figure 15).This is because the incoming turbulence disintegrates the tip vortex system quickly.For ABL cases, peak −u r u x occurs around one rotor diameter downstream of the turbine near the ground (θ = ±180 • ), while near the 12 o'clock position, the peak values occur further downstream.While there is a significant difference in the spatial distribution and magnitudes of u r u x between the three inflow conditions considered, the differences between the SRWT and DRWT are subtle.The DRWT shows a higher level of momentum entrainment, especially for x > 5D.To calculate net entrainment from all around the cylinder, u r u x is further averaged azimuthally.The azimuthal averaging operation is denoted by angle brackets; thus, the quantity u r u x represents the temporallyand spatially-averaged value of u r u x .u r u x multiplied by the circumference of the cylinder (2πr tip ) is the turbulent momentum entrainment per unit length into the cylinder.Subplots (g-i) in Figure 15 show the variation of u r u x (normalized by u 2 ∞,h ) with distance from the turbine rotor location.
The y−axis is reversed in the plots as negative values imply positive entrainment.Variation is plotted on the same scale to contrast the mixing rates between the different inflow conditions.The neutral ABL case shows the highest entrainment, while the uniform inflow case shows the lowest entrainment.The plots for the turbulent flux of kinetic energy (u x × u r u x ) are very similar to those for momentum flux and, hence, are not shown.
The overall percentage changes in turbulent flux of axial momentum and kinetic energy due to the DRWT are provided in Table 2.The increase in entrainment for the DRWT over the SRWT is highest for the neutral ABL case and lowest (indeed there is a slight reduction) for the uniform inflow case.A modest 3.29% increase in momentum entrainment is observed for the neutral case.The reader is reminded that in these simulations, the secondary rotor is operated at the tip speed ratio that gives the maximum aerodynamic performance for isolated turbine operation in uniform inflow conditions.No attempt is made here to optimize the secondary rotor design/operation to enhance wake mixing.Notwithstanding the sub-optimal design/operation, the DRWT still shows higher levels of entrainment than the SRWT for the two ABL cases.These results demonstrate that in a wind farm, the turbines operating in wake flow can benefit from enhanced entrainment provided by the DRWT.Additionally, the fact that the turbine operation and the geometry of the secondary rotor can be optimized to actively target wake mixing leaves potential for further improvement in wind farm efficiency from the DRWT technology.

Aerodynamic Loads
Aerodynamic loads in terms of blade root bending moments are analyzed in this section.The approximations made to carry out this analysis are: (1) the rotor blades are assumed to be infinitely rigid (i.e., no deformation in turbine geometry is permitted); (2) the turbine is operated at a fixed rotation rate irrespective of the incoming wind speed (hence, the tip speed ratio is fluctuating due to inflow turbulence); and (3) the blade geometry is not resolved in the simulations, and hence, the potential interaction between the rotors due to blade thickness is not captured.
Figures 16 and 17 compare the turbine power and out-of-plane blade root bending moment (M OOP ) in the time and frequency domains for the two ABL cases simulated.The dynamic loads in the uniform flow case are very small compared to the ABL cases, and hence, those results are not presented.For the DRWT, turbine power is the sum of powers from the two rotors; blade root moment however is only compared for the main rotor unless otherwise stated.Figure 16 compares the DRWT and SRWT power and loads for the stable ABL condition.The secondary rotor in the DRWT efficiently extracts power near the main rotor blade root region, and hence, the DRWT produces higher net power than the SRWT (see Subplot "a" in Figures 16 and 17).In the plots, power is normalized by 1/2ρ ū3 ∞,h A d and blade root bending moments by 1/2ρ ū2 ∞,h A d r t , where A d = πr 2 t and r t is the main rotor tip radius.
Fluctuations in power output are caused by spatial and temporal variations of the incoming wind.While the temporal variations are only due to turbulence, spatial variations also occur in the mean wind speed in the ABL; the mean wind speed varies monotonically with height from the ground in the area swept by the turbine rotor (see Figure 6a).Each blade of the turbine rotor(s) experiences a one per revolution variation/excitation because of this spatial variation of the mean wind speed -highest wind at the 12 o'clock position and lowest at the six o'clock position.This periodic (since the rotor rotation speed is fixed) excitation results in deterministic fluctuations in turbine power and blade loads.The torque contribution from each blade of a turbine rotor adds linearly (scalar addition) to give the net rotor torque (and power); blade root bending moments however are about different axes for each blade, so they need to be added as vectors in order to compute the moment on the turbine shaft.Here, we analyze blade root moments.Since the turbine blades are assumed to be identical, one-per-revolution fluctuation in the torque/power for each blade causes a one-per-BPF variation in torque/power for the rotor, where BPF stands for blade passing frequency.Blade root bending moments however are for each rotor blade, which therefore fluctuate at the fundamental frequency of one per turbine revolution, or 1/rev.Deterministic power fluctuations at the BPF and its harmonics are observed in the power spectral density (PSD) plots (see Subplot "c" in Figures 16 and 17) due to the variation of mean wind speed with height as discussed above.In these figures, the main rotor passing frequency is used to normalize frequency.Power fluctuations at secondary rotor passing frequency and its harmonics are also present, but they are much smaller in magnitude and are therefore inconspicuous in the figure.While the increase in mean power from the DRWT is measurable (between 5% and 6%), the change (reduction) in power fluctuations observed is insignificant.Power fluctuations are not desirable, and hence, no increase in such fluctuations with the DRWT is beneficial.The difference in out-of-plane blade root bending moment (M OOP ) is also insignificant; the DRWT showing a very small reduction at high harmonics of rotor passing frequency (see Figure 17).The reduction is expected to be small since the contribution to the root bending moment from the blade root region, where the blade relative flow velocity and the moment arm are both small, is insignificant.The results suggest that the aerodynamic interaction between the main rotor and the secondary rotor has little impact on unsteady aerodynamic loads experienced by the main rotor.
Figure 18 compares the M OOP for the SRWT operating in the stable and neutral ABL conditions.The PSD of the M OOP is higher for the stable case at the fundamental rotor passing frequency because the vertical shear in the mean wind speed is higher for the stable case (see Figure 6a).Vertical shear in mean wind results in a 1/rev variation in the angle-of-attack, and hence, loads, on the rotor blades.At frequencies greater than the fundamental rotor passing frequency, the neutral case shows higher values of M OOP than the stable case.This behavior can be explained by the difference in the size of the turbulent eddies in the two ABL cases.A wind turbine rotor can "chop" through a large, slow-moving turbulent eddy multiple times, resulting in blade loads at the rotor passing frequency and its harmonics.The integral length scale represents the size of the largest eddies in a turbulent flow.In the present problem, the integral length scale is computed at the hub height using a two-point correlation of streamwise wind velocity, defined as R uu (r) = u(x) u(x + r êx ) , where the angle brackets denote spatial averaging over the horizontal plane at the turbine hub height (x.êz = h) and the overline denotes time averaging.The two-point correlations are computed using the precursor simulations for the two ABL cases and are compared in Figure 18b.The integral length scale (L) is computed as L = R max 0 R uu (r) dr, where the upper limit of the integral R max = 3 km.The integral lengths are found to be approximately 24 m and 130 m for the stable and neutral ABL cases, respectively.The higher values of M OOP at frequencies greater than the fundamental rotor passing frequency can be attributed to the larger eddy sizes (higher integral length scale) in the neutral ABL case.

Secondary Rotor
Figure 19 presents the power spectral densities of the secondary rotor M OOP for neutral and stable conditions.Contrasting the magnitudes in Figure 19 with those in Figures 16 and 17 elucidates that the mean and the fluctuating loads on the secondary rotor are in fact orders of magnitude smaller than those on the main rotor.Therefore, the secondary rotor can be designed with relatively thin airfoils.The abscissa in Figure 19 is normalized by the secondary rotor blade passing frequency.While the decay in power spectral density with frequency is monotonic for the main rotor, which is not the case with the secondary rotor.The M OOP at the fourth harmonic of the rotor passing frequency is higher than the second and third harmonics.This is likely an effect of the aerodynamic interaction between the main and the secondary rotors.

Conclusions
A numerical study is conducted using large eddy simulations to investigate aerodynamic performance and loads for the dual rotor wind turbine (DRWT) technology proposed by Rosenberg et al. [16].The LES solver is first validated against experimental data, as well as blade element momentum theory results for a conventional, single-rotor turbine.A two-step procedure is used to first simulate the atmospheric boundary layer and then wind turbine aerodynamics.The DRWT is analyzed for three different inflow conditions: uniform inflow, stable and neutral atmospheric boundary layer.The following conclusions are drawn from the results of this study.
1.The DRWT operating in isolation shows aerodynamic performance (C P ) improvement of about 5%-6% for all inflow conditions.The performance benefit is obtained due to efficient extraction of energy (using the smaller secondary rotor) from the streamtube going through the blade root region.2. The DRWT enhances wake mixing and entrainment of higher momentum fluid from outside the wake layer when the atmospheric (freestream) turbulence is moderately high (as in the neutral stability case simulated here).A modest (≈ 3.2%) increase in momentum entrainment is observed with the de facto DRWT design.3. The enhancement in wake mixing is associated with the increased turbulence intensity in the wake.This could potentially increase fatigue loads on the downstream turbines.4. Spectral analysis of aerodynamic loads (measured as rotor power and out-of-plane blade root moment) shows negligible reduction for the main rotor in the DRWT.Unsteady fluctuations in rotor power are observed at blade passing frequency, while fluctuations in blade root moments are at the rotor passing frequency and its harmonics.These fluctuations occur because of the azimuthal variation (due to the ABL) in the incoming mean wind, as well as turbulence in the wind.
The enhanced wake mixing rate observed with the DRWT is promising.It proves that the technology has the potential to improve wind farm efficiency if the secondary rotor is optimally designed and/or operated for that purpose.

Figure 2 .
Figure 2. Verification of LES predicted results for the Tellus (Risø) turbine.(a) Non-dimensional rotor blade chord and twist variation; and (b) power curve variation with wind speed compared against data and blade element momentum (BEM).
(a) CFD domain for precursor simulations (b) CFD domain for wind turbine simulations

Figure 4 .
Figure 4.A schematic showing the computational domains for the atmospheric boundary layer (precursor) and the main (wind turbine) simulations.The entire box in (a) is the domain for precursor calculations; the shaded area shows the smaller domain for wind turbine calculations.In (b), Points A-D are lateral midpoints of the rectangular refinement zones.

Figure 5 .
Figure 5. Radial distributions of the (a) primary and (b) secondary rotor blade chord and twist of the dual-rotor wind turbine (DRWT) configuration analyzed here.

Figure 6 .
Figure 6.Time averaged streamwise wind speed and streamwise turbulence intensity profiles for the two atmospheric boundary layer (ABL) conditions simulated.The mean wind speed is ū∞,h = 8 m/s.The turbulence intensity is zero for the uniform inflow case.

Figure 7 .
Figure 7. Iso-surfaces of the Q-criterion of the SRWT operating in a neutral ABL.The contours are colored by streamwise turbulence intensity: red and blue showing high and low turbulence intensity levels, respectively.

( a )Figure 8 .
Figure 8.Comparison between SRWT and DRWT of normalized mean streamwise velocity, u x /u ∞,h : (a,b) on the x-z plane passing through the rotor hub; (c,d) on cross-stream (y-z) planes; and (e) circumferentially averaged values at four downstream locations (x/D = 2, 4, 6, and 8) for the uniform inflow condition.

( a )Figure 9 .
Figure 9.Comparison between SRWT and DRWT of the normalized mean streamwise velocity, u x /u ∞,h : (a,b) on the x-z plane passing through the rotor hub; (c,d) on cross-stream (y-z) planes; and (e) circumferentially averaged values at four downstream locations (x/D = 2, 4, 6, and 8) for the stable ABL condition.

( a )Figure 10 .
Figure 10.Comparison between SRWT and DRWT of normalized mean streamwise velocity, u x /u ∞,h : (a,b) on the x-z plane passing through the rotor hub; (c,d) on cross-stream (y-z) planes; and (e) circumferentially averaged values at four downstream locations (x/D = 2, 4, 6, and 8) for the neutral ABL condition.

Figure 14 .
Figure 14.Investigation of turbulent momentum transport into the turbine wake layer: (a) a schematic; and (b) cylindrical surface showing the turbine wake layer through which turbulent momentum flux is computed to quantify entrainment in the turbine wake.

Figure 15 .
Figure 15.Radial transport of streamwise momentum into the turbine wake layer.Color contours show u r u x /u 2 ∞,h on the cylindrical surface of Figure 14 that has been cut at θ = ±180 • and unwrapped.

Figure 16 .
Figure 16.Stable ABL simulation results.Turbine power and out-of-plane blade root bending moment (M OOP ) compared in the time and frequency domains.

Figure 17 .
Figure 17.Neutral ABL simulation results.Turbine power and out-of-plane blade root bending moment (M OOP ) compared in the time and frequency domains.

Figure 18 .
Figure 18.Autocorrelation of the velocity fluctuation as a function of distance.

Table 1 .
Percentage difference in time-averaged turbine C P and M OOP .

Table 2 .
Percent change (DRWT-SRWT) in cumulative streamwise momentum and kinetic energy for the three inflow conditions.Inflow Condition %∆(u r u x ) %∆(u x × u r u x )