Numerical modelling and optimization of vertical axis wind turbine pairs: A scale up approach

The performance augmentation of pairs of vertical axis wind turbines (VAWTs) is known to be dependent on incident wind direction, turbine spacing and direction of rotation. Yet, there is a lack of robust numerical models investigating the impact of these parameters. In this study two-dimensional CFD simulations of an isolated VAWT and of coand counter-rotating pairs of VAWTs were performed with the aim to determine turbine layouts that can increase the power output of VAWT farms. More than 11,500 h of simulations were conducted at a turbine diameter Reynolds number of 1.35 $ 107. A mesh convergence study was conducted, investigating the influence of mesh size, domain size, azimuth increment, number of iterations per time step, and domain cell density. Results showed that mesh size, domain size, and azimuth increment proved to have the biggest impact on the converged results. For the configurations analysed, pairs of VAWTs exhibited a 15% increase in power output compared to operating in isolation, when the second rotor was spaced three turbine diameters downstream and at an angle of 60 to the wind direction. Furthermore, when three turbines were positioned in series, the power output was greater than a pair by an additional 3%. © 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The UK's wind capacity is expected to almost double by 2030 [1]. So far, all large scale wind farms (>40 turbines) are utilising horizontal axis wind turbines (HAWTs), and these are continuously becoming more efficient and larger in size [2] in order to maximise the energy extracted from the given site. Yet, turbulent wakes created by the first row decrease the power output of the turbines behind by up to 40% [3]. Vertical axis wind turbines (VAWTs) could solve this problem since research [4,5] has shown that this type of turbine exhibits the opposite behaviour when composed in wind farms and apparently they enhance each other's performance. Furthermore, maintenance costs are lower due to fewer moving parts, which also makes them easier to install, and, opposite to HAWTs, they can be installed at sites with varying flow conditions (i.e. varying wind direction) [6,7]. Their primary disadvantages are the lower efficiency, reaching 35%e40% in isolation compared to HAWT's that is near to 50% [8], and the low starting torque for some designs e i.e. external power is required to accelerate the turbine at small angular velocities up to its optimal tip speed ratio, l. In summary, harnessing wind power coming from any direction to create energy using VAWT is an attractive option and despite the intensive research in the field, the underlying performance parameters of VAWTs are not well understood. To develop wind farms that can meet future energy demands, performance optimization of VAWT farms is required.
There are studies in literature [2,9e13] investigating optimal VAWT blade designs and geometric properties, with the primary tool being Computational Fluid Dynamics (CFD) and/or wind tunnel experiments. One of the few studies that investigated performance augmentations of VAWT farms was that of Dabiri [5], where experiments in a desert with six 10 m tall times 1.2 m diameter VAWTs were conducted. The experiments investigated the effects of turbine spacing and direction of rotation. It was observed that while HAWTs experienced an overall decrease in power by 20%e 50% when placed in close proximity to each other (1.65 turbine diameter separation), the VAWTs enhanced the overall performance by 5e10%. Furthermore, when the spacing was 4D (D ¼ turbine diameter) for the VAWTs, the downstream rotor exhibited a 5% lower deviation from its isolated performance. This held true for other tip-speed ratios too. Results were also in contrast to HAWTs where a turbine spacing of 15D-20D downwind to fully recover the wake was required [5].
Parneix et al. [14] is an example of another study, where two side-by-side VAWTs were analysed; the CFD model used an inhouse code based on vortex methods to estimate power coefficients. The two-dimensional simulations showed that for a turbine spacing of 1.2D, the power coefficient was 15% higher for both turbines. Additionally, the three-dimensional simulations found an increase of at least 8%.
Finally, Brownstein et al. [4] conducted experiments of a pair of 5-bladed VAWTs in an open circuit, subsonic wind tunnel. The rotor diameter (D) was 0.20 m, and the blades had a NACA 6415 aerofoil shape. The study investigated the rotor performances for varying; array angle (b) between À90 and 90 , direction of rotation, and turbine spacing between 1.25D to 3D. The Reynolds number was 7:3,10 4 for all experiments, and this is more than a factor of 100 lower than what a larger VAWT would experience offshore. Three distinct regions in b were found, and for ba30 , both turbines exhibited performance augmentations between 1.1% and 12.5%. In this regime, a turbine spacing of 1.25D displayed the best improvements.
In this study, the power improvements of a pair of VAWTs configured in two-dimensional CFD simulations for varying; array angle, direction of rotation, and turbine spacing were investigated. To the authors' best knowledge, this is the first attempt at numerically investigating the efficiency augmentations of VAWTs for more than 20 different layouts. The results were compared using a dimensionless parameter and validated against experimental data in literature [4].

Methodology
A frequently applied method for analysing VAWTs is CFD simulations, and in this work, the software package Simcenter STAR-CCMþ 2019.2.1 was used. Three-dimensional CFD simulations are computationally expensive and Bianchini et al. [2] concluded that two-dimensional VAWT studies give accurate results if reasonable mesh, timestep, and geometry settings are applied; thus, it was decided to perform two-dimensional URANS CFD simulations. The CFD simulations were transient, because of the rotational motion, and the implicit method was applied since it was more numerically stable. Balduzzi et al. [9] used a coupled algorithm, however, there were negligible affects when using a segregated solver for small time increments. Coupled flow requires 1.5e2 times more memory, thus, it was decided to use the segregated flow algorithm to develop a computationally efficient model. When wind blows over a VAWT, it occurs at very low Mach numbers and therefore the compressibility effects are minimal to null. As a result, it can be assumed that the air has constant density and is incompressible. k-ε is the most frequently applied turbulence model for analysing VAWTs [9], however, recent papers within the field of VAWTs [9,15] recommended applying the SST model instead of the k-ε model as it gives results that are closer to the experimental data. Hence, in this study the SST (Menter) k-u (SSTKO) turbulence model was applied, and the model differentiated from the standard k-u model (SKO) in the formulation of the production terms P k and P u in the two transport equations (Eqs. (1) and (2) For P k , an additional non-linear production term (Eq. (3)) is added in the SSTKO model, and for P u , a cross-diffusion term (Eq. (4)) is included [16].
The time step of 0.003155 s was derived from a sensitivity analysis during the mesh convergence study, and further information can be found in supplementary material S.4. The optimal setup parameters within the range investigated balancing accuracy and solving time, are stated in Table 1. The mesh convergence study was solved by using the supercomputer at Oxford Advanced Research Computing for a total of 2458 h. The layout simulations were done on the pooled computers at Oxford Brookes University, and the total simulation time was >9000 h.

Rotor geometry
VAWTs can either be comprised of asymmetric or symmetric aerofoils. The classic NACA four-digit symmetric aerofoils have been used by many researchers, because they have been thoroughly tested. Islam [17] performed a numerical study of a wide range of asymmetric aerofoils and compared them to a symmetric NACA0015 profile. The results showed that the asymmetric aerofoils exhibited the highest power coefficients at low tip speed ratios (l < 3), whereas the NACA0015 had a higher efficiency at greater tip speed ratios (l > 3). Brusca et al. [18] concluded that a NACA 0018 was the most efficient aerofoil, and its optimal tip speed ratio occurred around 3.5. Fig. 1 shows the three-bladed design and its geometric dimensions applied in this study.

Domain geometry
The flow must be fully developed over the length L, and the width W must be enough to avoid the boundary effects near the walls affecting the flow in the middle [19]. A domain size that is too small, does not allow the flow to fully develop. Balduzzi et al. [9] concluded that the domain size had to be greater than W ¼ 60D Â L ¼ 90D to replicate open-field-like boundaries (Fig. 2). Furthermore, the first rotor (R1) was positioned L 1 ¼ 20D from the inlet, because the wake was the most important feature and it stretched longer than what occurred in front of the turbine.

Boundary conditions
The rotors were positioned in a rotating region within a larger rectangular domain with stationary walls top and bottom, a velocity inlet, and a zero-gauge pressure outlet (Fig. 3). The rotating regions had a diameter, D RR , of 1.5D. An overset mesh was created between the rotating region and domain using STAR CCMþ, and a slip-condition was applied to the walls to avoid blockage effects. Verified by data from Bockstigen offshore wind farm, Gotland, Sweden the inlet velocity, U 0 , was 10 m/s and the turbulence intensity was 1% [20]. Angular velocities, u, were set at 3.5 rad/s according to Ref. [18]; thus, the turbine diameter Reynolds number was 1.35 $ 10 7 . 25 different layouts were analysed and evaluated with 24 being with two rotors and 1 layout with three rotors. Three variables were investigated; the turbine spacing, dist, the array angle, b, and the direction of rotation of R2 with R1 always rotating in a counter-clockwise direction.

Meshing strategy
A mesh convergence study was conducted using the supercomputer at Oxford Advanced Research Computing for a total of 2458 h. The final mesh predicted moment coefficients that deviated 2.1% from the finest mesh, however, it solved 32 times quicker. Fig. 4 illustrates the final mesh with (a) structured hexahedral cells in the domain including refinement boxes around the wakes derived from the findings and recommendations of [9,21]. Circular rotating regions with triangular cells (b), and tall, dense prism layers near the aerofoils (c) to accurately resolve the boundary layer. The grid convergence study concluded the mesh properties to be 150 prism layers with a first layer height of 5.41 $ 10 À2 mm and prism stretching factor of 1.03. Surface curvature of the blades was set to be 1000 points/circle, and the base size of the rotating region was 0.16 m. The base size of the; domain was 16 m, largest refinement box was 1.6 m, and smallest refinement boxes were 0.2 m. The surface growth rates, 1.05 in the rotating regions and 1.25 in the domain, were kept low to ensure a smooth transition from near the walls to the freestream regions. There were approximately 340,000 cells in each of the rotating regions, and

Convergence
The averaged moment coefficients over one revolution lowered as more revolutions were solved. Previous studies were typically concerned with the percentage change between the moment coefficients of two consecutive revolutions, and the most stringent criterion in literature was the one suggested by Balduzzi et al. [9], who concluded a simulation to be converged when the difference was below 0.1%. All results showed that R2 took the longest to converge, hence why studies regarding convergence criterion for a VAWT in isolation were not strictly applicable to a pair in all instances.
It was found that performing regression analysis using a rational function with degrees between 1 and 15 gave the best results at predicting the moment coefficient. The rational function describing the average moment coefficient, M avg , as a function of number of revolutions solved, was: where N R was the revolution number, m and k indicated the degree of the polynomials in the numerator and denominator respectively, and j and z were coefficients. The best fitting functions occurred for m ¼ k. The mathematical software Maple 2019 was utilised to determine the fitting equations. As initial revolutions have minimum effect on the final convergence, the first 10 revolutions were not taken into account, when using the fitting equations. Simulations were solved for an average of 56 revolutions, and the data points were smoothened for three consecutive revolutions.
The weighted mean square error was applied to evaluate the applicability of each fitting equation. The weighting was derived from the number of revolutions, so that the last revolutions had a greater importance than the initial values. The weighted mean square error, WSME, was calculated as: Where M n was the average moment coefficient for the n'th revolution, d M n was the moment coefficient predicted by the fitting equation at the n'th revolution, and jnj was the cardinality of set n.

Performance indicator
The parameters most frequently used to evaluate the performance characteristics of VAWTs are; the instantaneous moment coefficient (C m ) e indicating the torque generated by the blades e and the coefficient of power (C P ) e indicating the energy efficiency of a turbine. The instantaneous moment coefficient is defined as; where M is the instantaneous moment, r is the air density, A is the cross-sectional area, R is the turbine radius, and U 0 is the freestream velocity. The power coefficient is where P is the power of the turbine, and it is defined as Combining Eqs. (7)e(9), one will obtain [22].
This expression (Eq. (10)) thereby correlates the three most important parameters (tip speed ratio and coefficients of power and moment) in the simplest manner. In this project, all rotors had equal angular speed, and since one was interested in the performance augmentations, the following performance indicator (U) was defined by Eq. (11): where C P_i was the coefficient of power of each rotor, C P_iso was the performance of a single rotor in isolation, and n was the number of turbines in the layout. Therefore, if U > 1, the configuration exhibited a higher power output than if the two turbines were operating on their own.

Rotor orientation
The average performance of co-rotating and counter-rotating pairs for a turbine spacing of 2D, are plotted in Fig. 6 against the array angle b. For both directions of rotation, the efficiency was greatest for 45 jbj 90 , but all compositions except for b ¼ 0 (R2 behind R1) had an efficiency above 1. Therefore, two VAWTs increase each other's efficiencies when positioned in specific layouts controlled by beta. Counter-rotating turbines showed to produce greater power than co-rotating at smaller (jbj 30 ) and the opposite at larger (jbj > 30 ) angles of b. Finally, whether the second rotor was positioned above, b > 0, or below, b < 0, only had a marginal impact on the performance when jbj ! 45 .
Velocity diagrams of ten layouts are illustrated in Fig. 7. When R2 was located in the wake of R1, (b) and (h), R2 did not experience the same kinetic energy of the wind compared to in isolation, hence there was a drop in performance. This occurs for HAWTs in wind farms too [3,5], however, their wake is more persistent (i.e. longer) and therefore the turbines must be spaced further apart. It can also be observed that the blades in (h) were only shadowed for half of their revolution, and as a result, the performance was nearly a factor of 2 greater in comparison to (b), which was fully covered by R1's wake. Fig. 7 (f) and (g) show that if R2 was at the periphery of R1's wake, then R2 bended the wake of R1 causing a greater power output of R2, and a lower power output of R1. This is also seen later in Fig. 8 for the co-rotating plot that the performance of R1 at 30 (R2 at the border of R1's wake) was lower than at 0 (R2 did not lie on R1's wake boundary) before increased performance augmentations again at jbj > 30 (R2 was away from R1).
For larger angles of b; R2 was not in R1's wake e see Fig. 7 (c), (d), (e), and (j); and the pair experienced augmentations. Fig. 7 (c) indicates the distorted flow field around the turbines too. For example, the red regions outside the wake had a higher velocity than the freestream region. Thus, the turbine caused the flow to accelerate.
It was not only the second rotor that was affected; the first rotor also experienced either an increase or decrease in power output. Fig. 8 shows the normalised performance values for the first and second turbines. R1 experienced a greater efficiency as b approached 90 ; yet, R2 reached a peak at jbjz60 . However, as shown previously in Fig. 6, the optimal layout for a pair as a unit occurred at b ¼ 75 . Co-and counter-rotating configurations gave similar results, and the co-rotating configurations exhibited the greatest performance for its R2 at À60 , where the performance augmentation was 1.231. Similarly, the smallest value occurred at 0 , giving a performance of À0.024, i.e. power was required to rotate the turbine. Fig. 8 also shows a spike for R2 at À 90 b À 30 , whereas at 30 b 90 the performance indicator U curve is flatter (it is clearest to see for the co-rotating layout). A hypothesis is that R1 was always rotating in the CCW-direction, and similar to an aerofoil generating lift by redirecting flow momentum downwards and away, then R1 acted in a similar way. This phenomenon is illustrated in Fig. 9. As R1's wake was directed towards the bottom right corner, then it affected the results, because R2 was now closer to the border of R1's wake at À30 than at 30 . Therefore, the performance augmentation U was greater at 30 . Even though, the performance was greater at b ¼ À45 than for 45 , then this is likely to be due to simulation tolerances.

Validation of numerical model
The results from Fig. 8 were compared to Brownstein et al. [4] experimental data as shown in Fig. 10. At this point it should be noted that in Ref. [4] the parameters do not exactly match the conditions of the numerical study such as number of blades and Reynolds number, however this was the closest available to our model study that can be used for validation. The average deviation was 5.1% and 7.8% for respectively R1 and R2, when only considering the results obtained for jbj > 30 . The deviation was greater if all angles were considered (R1: 7.1% and R2: 24.7%), however, these are not of interest when looking for the optimal position due to poor performance. The larger deviation for small angles of b, is hypothesized to be due to three-dimensional effects not being included in the two-dimensional CFD-simulations; e.g. how wind from above may travel to lower altitudes to recover the wake quicker. Furthermore, the asymmetry of the upstream turbine was not as dominant as the one found by Ref. [4], however, the wind tunnel experiment in Ref. [4] with a 5-bladed VAWT design was conducted at a Reynolds number that was a factor of 185 times lower than the one of this work. Nevertheless, two-dimensional CFD simulations are applicable at predicting the performance augmentations experienced by pairs of VAWTs for jbj > 30 . Finally, results were in a very good agreement with [4] where the experimental results showed a spike at bz À 40 , and a flatter curve for b ! 30 .
As described in section 3.1 Rotor orientation, the shift in the performance indicator U was due to the wake being directed downwards. A possible reason to why the numerical results did not show a significant shift was that [4]'s 5-bladed turbine pushed more fluid downwards. Furthermore, the wind speed was nearly a factor 2 (as more representative values monitored by Larsen, G. C. and Hansen, K. S [20]. replicating real working conditions were deployed) lower than this CFD study, and therefore the air had less momentum to 'straighten' the wake. Although, there was a mean 6.5% error in values, then the overall trend of the performance behaviour was captured very well by this methodology. Fig. 11 shows that turbine spacing also proved to influence the performance. The power enhancement increased as R2 was positioned further away from R1. Furthermore, near-all configurations showed an increase in power as b approached 60 . It appears that the enhancement for increasing turbine spacing decayed as b approached 90 . This might indicate that theoretically the optimal distance between rotors was z3D. Nevertheless, the results argue that for larger angles of b, the turbines can be packed closer together without considerable loss in performance. Hence, for a    turbine space of 2D, more turbines can fit within the same area and still generate power close to that of 3D.

Effect of multiple turbines
If the number of rotors in series was increased to three turbines, the results showed that the efficiency for the whole system increased from 1.05 to 1.09, i.e. a further power improvement of 4% points. However, the R1 and R2 had a lower efficiency, yet, it was the high efficiency of R3 that generated a greater overall power output (Fig. 12). In other words, the average performance of R1 and R2 was 0.99, but R3 brought up the average performance to 1.09.

Torque profiles
In order to understand the fluid accelerations caused by R1, the blade torque profiles of the turbines over one revolution were compared (using the results of the last revolution in the simulation). The graphs for two of the layouts are shown in Fig. 13 and Fig. 14. These were particularly chosen because the first (Fig. 13) exhibited a low power output of U ¼ 0.81. R2 experienced higher torque for 0 q 100 compared to operating in isolation, thereafter the moment coefficient was significantly lower up until q z 300 . For R1, the torque was lower throughout the whole revolution. In summary, these lower torque-values caused the layout to perform a power output that was 20% lower compared to isolation.
In Fig. 14, the opposite was true with the turbine experiencing a power greater than isolation conditions, U ¼ 1.10. The performance augmentation was due to the moment coefficient of R2 being greater compared to the isolation case for near all angles of q. For R1, the greater moment coefficient occurred primarily between 120 q 180 and 240 q 360 . Hence, the performance enhancements occurred due to a higher torque during predominately the upwind stroke of the rotation, i.e. 0 q 180 .

Velocity across rotor span
The presence of a second turbine generated changes to the mean flow velocity field. Fig. 15 shows the variation of velocity across the rotor span for varying array angle, b. The velocity was averaged across the vertical (y-direction) rotor diameter span in the middle. Moreover, the velocity magnitudes were averaged over one full revolution, U avg , and then normalised by dividing by the results of a turbine in isolation, U iso . The improvements/deficits in average velocity travelling through R1 and R2 were similar to the performance augmentations previously shown in Fig. 8. For example, the normalised average velocity was greater for R2 than R1 for both cases, and the trend that the incident flow speed on R1 increased as the array angle approached 90 .
The lower output of R1 and R2 was likely to be due to R3 negatively interfering with their wakes. Instead, if the tripleconfiguration utilised a larger turbine spacing and let bz75 (ref. Fig. 6), hypothetically, the wake interaction would be more beneficial, and the fluid would continue to accelerate for each rotor downstream. However, the curve must flatten out at some point, and the wind energy available at the site could potentially be the limiting factor in this case.

Discussion
Mesh sensitivity analyses were carried out for five different parameters (mesh size, time step increment, domain cell density, number of iterations per timestep, and domain size), and these showed that the domain size, mesh size, and time step increment were the values that had the biggest influence on the average moment coefficient within the ranges evaluated. Yet, mesh convergence studies for wind turbines are critical, but they do not imply universal applicability, due to factors such as leading edge erosion and changes in flow conditions.
Results confirmed the potential of VAWT farms, since close-toall layouts experienced performance augmentations within turbine spacings that are not achievable with HAWTs. Interestingly, Fig. 6 indicated that the efficiency improvements in the range of 1.00e1.15 occurred for a broad range of angles, 30 jbj 90 . In other words, this indicates that there should not be a significant decrease in power output if the wind direction changed from north to east. Of the 25 different layouts investigated, a turbine spacing of 3D, array angle of 60 , and both rotors spinning in the same direction exhibited the greatest improvement in efficiency by a 15% increase. This layout was at the limit of the scope, therefore, as evidenced by the other results, which exhibited an optimal angle around 75 , the augmentations are likely to be higher if the 3D layout had an array angle of 75 .
The layout being least prone to changes in wind direction would be~83 , because if there was a slight change in wind direction, the pair would still exhibit performance augmentation. Fig. 16 depicts this argument, and in this diagram b was 90 (a) for illustrating purposes and added simplicity. In configuration (b) the wind direction changes by 15 , and therefore the effective array angle is 105 , which corresponds to the layout analysed with b ¼ À 75 .
According to Fig. 6, this layout (turbine spacing of 2D) resulted in U ¼ 1:10, hence performance augmentations would occur. The same is true if the wind direction changed by À15 . Finally, an array angle of~83 is least prone to changes in wind direction, due to the results indicating that the power enhancements were greater at 30 than À30 . Therefore, the interval in which U > 1.05 happened from 30 b 90 and À 90 b À 45 , and the middle of this region is approximately 83 .
The aerofoil of a VAWT rotated with a higher speed than the fluid surrounding it, thus in the region where the blade was moving downstream, it accelerated the flow around it. On the other hand, when the blade was moving upstream it decelerated the surrounding fluid, and as a result, a distorted velocity field was established. As Fig. 7(c) indicates, there were red regions with accelerated flow and a wake, where the wind speeds were negligible. In the end, these fluid mechanical movements influenced the distribution of flow momentum in the vicinity of the turbines, which suggests that the power enhancements were caused by redirections of momentum near the rotors.
The improvements in incident flow speed lead to a discussion of the CFD setup. This work applied an equal constant angular speed   Torque profile over one revolution for a turbine spacing of 2D, array angle of À75 , and the turbines were co-rotating. The performance was U ¼ 1.10. Fig. 15. Normalised average velocity across rotor span against array angle. U avg was the average velocity over one full revolution for the given rotor, and U iso was the velocity across turbine diameter for a turbine operating in isolation.
to both rotors, and as previously stated, this was an assumption made to simplify the problem. Under real conditions, the angular speed will vary over one revolution, due to varying torque. Furthermore, as evidenced by this study, the varying mean flow field will further increase the variance in torque, and therefore one might question the assumption that rotors should rotate with a constant speed. One may say that the CFD simulation should be setup as an iterative procedure, i.e. the angular velocity of each rotor is changed for each time step to replicate realistic conditions. The tip-speed ratio will also not be equal for R1 and R2, because their incident velocity will be different. From the results, it can be derived that the turbines did not operate at their optimal tip-speed ratio, since a higher freestream across the rotor span of R2 at jbj ! 30 (Fig. 15) caused an improvement in power output (Fig. 8). All these measures increase the complexity of the problem, however, as evidenced by the similarity with the wind tunnel tests, CFD simulations still have validity, if setup correctly.
The CFD simulations demonstrated notable wake interactions. For larger wind farms, the suppressed wakes caused more energy to be available for the subsequent rows, thus maximising the power density. VAWT dynamics are strongly influenced by tip-speed ratio [4], and therefore the efficiency enhancements for a wide range of tip-speed ratios must be further investigated. The results obtained are strictly valid only for the given rotor geometry and boundary conditions simulated, since the quantitative properties are expected to depend on the Reynolds number regime [23].
In summary, it is generally accepted that the primary flow mechanism is proposed to be the flow acceleration around the upstream turbine (R1), which increases the incident wind speed on the downstream turbines (R2, R3). Additionally, the flow in the arrays is reenergised by the turbulent phenomena, and together, these mechanisms are proposed to be the reason why VAWTs exhibit improvements in power compared to operating in isolation [24].

Conclusions
In this paper 25 different layouts were investigated. Results show that VAWTs increase each other's performance by up to 15%, and this optimal layout was for a turbine spacing of three turbine diameters, an array angle, b, of 60 , and when the rotors were corotating. Yet, this layout was at the limit of the scope, and the other results indicated an optimal angle around 75 , hence the augmentations are likely to be higher for layouts in this region. Key findings of the study were: a. As jbj approached 90 , the performance of R1 increased. b. R2 peaked in power augmentation at jbjz75 c. The total efficiency increased as the turbine spacing increased d. Increasing the number of turbines further increased the overall efficiency. e. Two-dimensional CFD simulations produced accurate results when compared to wind tunnel tests, since the values were within 6.5% of experimental data for the augmented layouts. f. Greater performance of pairs was predominately due to a distorted flow field that was established in the vicinity of the VAWTs, and in these regions, fluid travelled with a greater speed than the freestream velocity.
In the future, the study has scope of expansion on bigger infrastructure by adding more turbines or going three-dimensional. Moreover, the rotor geometry was derived from numerical studies for boundary conditions and Reynolds numbers that were not exactly equal. Accordingly, there is great potential in conducting a design optimization study with the application of wind farms in mind to achieve further advances in performance. One might find that the optimal design for a turbine in isolation is not identical to the one optimised for wind farm configurations.
The potential applications for VAWTs are endless, because the turbines are cheaper and easier to manufacture and maintain. For example, remote villages and islands that primarily rely on electricity from diesel generators or off the coast of UK. The common factor is that VAWTs farm are likely to not be limited by the efficiency of the turbines, but by the wind energy available at the given site.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work. Prof. I. Tzanakis is thankful for the support received from the Engineering and Physical Sciences Research Council (EPSRC) through the UltraMelt2 (grant number: EP/R011044/1) project.