Analysis of the aerodynamic performance of the multi-rotor concept

The concept of a large (~20MW) multi-rotor wind turbine intended for offshore installations is analysed with respect to its aerodynamic performance. The effect of closely clustering rotors on a single actuator disk is estimated using two different modelling approaches: a CFD solver in which the rotors are simulated as distinct actuator disks and a vortex based solver in which the blade geometry is exactly considered. In the present work, a system of 7 rotors is simulated with a centre to centre spacing of 1.05D. At nominal conditions (tip speed ratio=9) both models predict an increase in power of ~3% alongside with an increase in thrust of ~1.5%. The analysis of the flow field indicates that in the 7 rotor system the individual wakes merge into one wake at ~2D and that flow recovery starts at approximately the same downstream distance as in the single rotor case. As regards the dynamic implications of the close spacing of the rotors it was found that there is an increase in the loading amplitude ranging from 0.30-2.13% at blade level in rated conditions.


Introduction
The FP7 INNWIND.EU project aims at the investigation and demonstration of new designs for 10-20MW offshore wind turbines and their components and the development of methodologies for asserting innovative subsystem and turbine system designs.
The design of a beyond-state-of-the-art 20MW offshore wind turbine multi-rotor system (MRS), mounted on a common structure support, is amongst the innovative concepts of the FP7 INNWIND.EU Project. This challenging design could lead to the overall reduction of the operation and maintenance costs, as well as the foundation and grid-connection costs in large scale offshore parks. A MRS might have much less weight and cost of blades and subsystems against a single equivalent rotor with the same swept area. Additionally, in a case of a rotor fault the MRS will lose only a few percent of capacity. The main challenges associated with the multi-rotor concept include: the development of a suitable support structure, the aeroelastic and control modelling of the whole system, the exploration of the MRS performance under a set of load cases.
A specific system of 45 rotors closely clustered in tandem on a single actuator plane has been proposed by the University of Strathclyde and structurally designed by CRES [1]. The main characteristics of this rotor are presented in Table 1. For the specific design and prior to its design assessment, a purely aerodynamic performance analysis has been undertaken by CRES and NTUA in view of better understanding the possible effects that could appear due to the close spacing of rotors. This part of the research is reported in the present paper.
The case study refers to a 7 MRS. Flow predictions have been obtained using a RANS flow solver based on the actuator disk modelling of the rotors and an unsteady vortex model that accounts for the blade geometry and in future will be the basis for the aeroelastic evaluation of the design. The main objective of the analysis here presented is to better understand the development of the flow due to excessive blocking caused by the close placement of the 7 rotors and thereby estimate the performance of the system (flow blockage refers to the "triangular" open spaces that form in between the rotors). Since vortex models are usually suppressing viscous effects by definition, it has been decided to also use CFD modelling in order to quantify the eventual viscous effects especially at the tip regions and by that get more confident in the predictions vortex modelling can deliver. In this connection, the overall performance, the corresponding flow field and the loading characteristics are presented and discussed.  [2] using an implicit pressure correction scheme and Wilcox's k-ω turbulence model [3] suitably modified for neutral atmospheric flows. The Reynolds-averaged Navier-Stokes equations for incompressible flows, are solved alongside with the transport equations for k-ω,  [4]. This corresponds to a value of 0.033    . Using the procedure described in [5], the following constants of the k-ω model are derived:

Wind Turbine modelling
It is common practice in CFD models to include wind turbines as actuator disks [6], [7], [8]. This type of modelling adds axial body forces in the corresponding momentum equation for the cells in which the rotor disks are contained. Let, ΔS denote the surface area of the numerical cell with width Δr (Figure 1), ρ the air density, U ax the local axial velocity and C T the local thrust coefficient. Then, the axial force is defined as The local C T is calculated using the blade element theory. Let U ax ,, U cir denote the axial and circumferential flow velocities at a cell containing part of the rotor disk. Since these components also contain the effect of induction, they also define the angle of attack α and through look-up tables the corresponding lift and drag forces L, D. By adding the local twist to α, the relative flow angle φ is obtained and thereby the axial force: φ=atan(U ax /U cir ), F x,n (r,ψ n )=Lcosφ+Dsinφ (6) The above expression refers to a point (r,ψ) of the rotor disk surface and is indexed differently for the three blades (n=1:3, ψ n =ψ+2π(n-1)/3). Rotor loading in the non-rotating frame can be obtained by applying a Fourier coordinate transformation that is also referred to as the Coleman's transformation [9], so that: Then in order to determine the local C T , F bl is interpolated linearly to the yz-plane and so,

Boundary conditions
At the inlet of the computational domain the velocity and k-ω profiles are given by the Monin-Obukhov theory [10] for neutral conditions: where K is the von-Karman constant, z 0 is the roughness length and U ref is the reference velocity at the hub height z hub . Here, z hub is the hub height of the central rotor. Neumann conditions are applied at the outlet boundary as well as at the side boundaries of the domain. The boundary layer height is set equal to 800m (20 D) and the velocity is considered constant above this height, U=U ∞ . As a result, the velocity is set equal to U ∞ at the top boundary and the turbulent kinetic energy is set equal to zero (laminar inflow). On the ground boundary all velocities are zero. Logarithmic wall functions are used to compute the velocity at the first grid point above ground. Figure 2. Numerical grid at the rotors plane x=0.

Computational grid
The computational domain is extended 20D downstream of the rotors level in the axial direction and 11D off the central rotor in the lateral direction. The height of the domain is 25D. In this way, the flow is not restricted by the numerical boundaries, where free-stream Dirichlet conditions are imposed. The inlet boundary is positioned 9D upstream of the rotors, so that a sufficient distance is provided for the flow to develop. The x-and y-axes are selected so that their origin coincides with the hub height of the central rotor.
The grid spacing in the x-direction reduces to its minimum value of 0.05D at the rotors level and increases outwards following a geometrical progression with a ratio of 1.15 until the maximum dimension of the domain is reached. In the vertical direction, the first grid-line is positioned at a height of 0.035D a.g.l., and then the grid spacing follows a successive coarsening and refinement up to the lower tips of the two lower rotors at a height of 4.267D. The grid spacing is kept uniform and equal to 0.04D from the lower tips of the lower rotors up to the higher tips of the two upper rotors and then follows a geometrical progression with a ratio of 1.15 up to the height of the domain (see Figure 2). In the lateral direction the grid spacing is kept uniform and equal to 0.03D all over the area of the 7 MRS, and it follows a geometrical progression with a ratio of 1.2 outwards up to the side boundaries of the domain. The mesh generated in this way comprises 75x137x175 grid points. The same mesh is used for the single rotor case in order to ensure a fair comparison of the predictions.

The GENUVP vortex method
GENUVP is an unsteady potential means of surface source and/or dipole distributions while the wakes are moving vortex blobs [11].The intensity of the surface distributions is determined by the no boundary condition. Along the trailing edge into the wake. The amount of released vorticity is determined by the pressure Kutta condition. The wake evolution is obtained by solving the vorticity transport equation in its The solution process starts at stand still and progresses in time following, sub steps shown in Figure 3. Loads are obtained by integrating the surface pressure distribution given by the solver. corrections can be introduced using is used. Only the blades are included in the simulations surfaces. A surface grid of 26(chordwise) corresponds to an azimuthal step of 7.5 at r=5.30m and therefore through In Figure 4 the blades of the 7 rotor (red dots) indicates the expansion of the wakes and their tendency to mix. Also an overall swirling of the wake under the influence of the mutual interaction of the wakes is noted. The intensity of the surface distributions is determined by the no the trailing edge, there is continuous vortex shedding that feeds vorticity into the wake. The amount of released vorticity is determined by the pressure Kutta condition. The wake evolution is obtained by solving the vorticity transport equation in its Lagr The solution process starts at stand still and progresses in time following, in every . The three sub steps per time step in GENUVP Loads are obtained by integrating the surface pressure distribution given by the solver.
using CL-CD data. In case of dynamic inflow, the ONERA is used. Only the blades are included in the simulations, which are approximated as thin lifting (chordwise)x13(radial) nodes per blade is used, step of 7.5 o . It is noted that the aerodynamic modelling of the through the hub region the flow is expected to accelerate In Figure 4 the blades of the 7 rotors and their wakes are shown. The distribution of vortex (red dots) indicates the expansion of the wakes and their tendency to mix. Also an overall swirling of the wake under the influence of the mutual interaction of the wakes is noted.
Perspective and back view of a 7 MRS. Red dots correspond to vortex blobs. In the Right figure the wake expansion, the tendency of the wake to mix and the overall swirling can be seen. and it follows a geometrical progression with a ratio of 1.2 outwards up to the side boundaries of the domain. The mesh generated in this way comprises 75x137x175 grid points. The same mesh is used flow solver in which the effect of solid boundaries is represented by ed by means of freely The intensity of the surface distributions is determined by the no-penetration , there is continuous vortex shedding that feeds vorticity into the wake. The amount of released vorticity is determined by the pressure Kutta condition. The Lagrangian formulation. in every time step, the three Loads are obtained by integrating the surface pressure distribution given by the solver. Viscous . In case of dynamic inflow, the ONERA model [12] which are approximated as thin lifting while the time step aerodynamic modelling of the blade starts the flow is expected to accelerate. and their wakes are shown. The distribution of vortex blobs (red dots) indicates the expansion of the wakes and their tendency to mix. Also an overall swirling of Red dots correspond to vortex blobs. and the overall swirling can

The test cases
The performance of a 7 MRS is examined, assuming the sheared inflow defined by (9), at nominal operating conditions (TSR=9) corresponding to a Re=3•10 7 for the CFD simulations. The single rotor case is considered in similar operational conditions at the same hub height (H=227m, which is the hub height of the 45 MRS proposed in [1]) and the corresponding results are used as reference values. The rotors are synchronized in speed and azimuth angle which would perhaps not be the case in real operation. It is also noted that neither pitch nor rotational speed control are active.

Multi-rotor system (MRS) performance
In Figure 5 the predicted thrust and power for the 7 rotor system is compared to the performance of the single rotor case. Both models predict an increase of similar level which suggests that the blockage of the flow is not expected to negatively affect the performance. In fact there is an increase of ~3% in power generation alongside with an increase of ~1.5% in thrust. It is noted that the single rotor and the central rotor of the 7 MRS were placed at the same hub height and the calculated loading of the single rotor was multiplied by 7 in order to compare it with the MRS system loading. The single rotor has the same rotor diameter as each of the rotors of the 7 MRS. Choosing the same hub height for both systems means that the single rotor is placed higher as compared to the standard 1D above ground level. Therefore in sheared inflow the single rotor will benefit from higher wind speeds and its performance will be somehow promoted. The individual wake cores of the 7 MRS merge into one wake at a distance which is comparable to the core length of the single rotor. In the single rotor case, flow recovery starts at ~100m, while in the 7 MRS case at ~300m giving a 1:3 ratio which coincides with the diameter ratio of the two systems. In non-dimensional terms the same decrease in the velocity deficit is observed at the same distance. Of course in dimensional terms, the extent of the 7MRS wake will be longer and therefore in a wind-farm design content, the next 7MRS system will be placed at a longer distance as compared to a design using single rotors. Assuming that the rows in a single rotor wind farm design contain 7 rotors, over the same wind farm area more rows could be fitted and more energy could be produced. However if the coverage of the wind farm area in the cross wind direction is also considered, more 7MRS would be placed. A clear comparison between 7MRS and single rotor wind farms would also require an account on the relevant costs and therefore a detailed investigation is needed before conclusions can be drawn in this respect. The effect of the multiple rotors on the wake deficit is also shown in Figure 7, where the axial velocity contours on the vertical plane y=0 are presented. Similar remarks can be made regarding the wake evolution. In these plots the effect of shear is also shown. Note that the lower part of the 7 MRS is closer to the ground level which explains the higher acceleration at z=150m as well as a slightly slower flow recovery. Left: Single rotor, Right: 7 multi rotor system Next, the flow field predictions are compared as obtained by the two models. Since the vortex model is an unsteady solver by definition, mean flow field predictions (averaged over a full rotation) are here used. Figure 8 presents a focused view of the axial velocity contours normalized by U ref at y=0 plane. The lack of viscosity in GENUVP is responsible for the long range conservation of the deficit indicated by the long blue strips in the plots. In the near wake the deficit gradually develops over a long range towards the theoretical U(1-2a) value given by the actuator disk theory. The acceleration seen over the gaps in between the rotors is also conserved as indicated by the red strips in the plots. Instead, in the CFD results, viscous mixing results in a different development of the velocity deficit and a quick damping of the acceleration obtained over the gaps. Similar remarks can be made on the results obtained on the z=227m hub height horizontal plane (Figure 9).

Loading characteristics
The loading characteristics of the multi-rotor system can be obtained from the output time signals of the vortex model. Snapshots of the normalized axial velocity at azimuth angle ψ=0 o are shown in Figure 11 on planes x=0 and x=0.1D. As expected, the axial velocity distribution over each of the system disks follows the rotating motion of the blades. The flow field outside the area occupied by the rotors, as a whole, is only slightly affected. On the contrary the acceleration within the triangular shaped gaps that are formed in between the individual rotors seems to follow the rotation of the blades. In fact this acceleration is triggered by the close passing of the blades. A similar pattern is seen at x=0.1D as regards the distribution of the velocity deficit which also follows the blade passage.  In terms of axial loading, the central rotor (CR_1) receives 2% more average loading and slightly higher amplitude of 0.3%, in comparison to the single rotor case. For the off-set rotors the increase in average loading varies from -0.22 to 2.91% depending on the vertical position of the rotor in the sheared inflow (Figure 12b), while the amplitude in the loading signals increases from 1.49 to 2.13%. These variations could have implications on the blade fatigue loads and should be further examined through aeroelastic simulations.

Conclusions
Two essentially different simulation procedures have been applied in view of better understanding and evaluating the aerodynamic implications of the close spacing in multi rotor systems. On one hand CFD simulations based on a consistent account of the rotor loading using actuator disk considerations have been used in order to obtain a realistic prediction of the performance and the mean flow characteristics of the system. On the other hand a vortex model that realistically accounts for the blade geometry has been used in order to assess the dynamic implications of such a system. By comparing the two sets of results at rated tip speed ratio, a good agreement was found which gives confidence to the predicted behaviour of a 7 MRS. It was found that there is no penalty in power. In fact an increase in power of ~3% was obtained alongside with an increase in thrust of only ~1.5%. This suggests that at least from an aerodynamic point of view, the 7 MRS will perform equally well with 7 stand-alone rotors.
As regards the dynamic implications it was found that at blade level, there is an increase in the loading amplitude ranging from 0.3 to 2.13%. So the next step will be to further evaluate the behaviour and response of multi rotor systems using aeroelastic simulations.