Vertical-axis wind-turbine farm design: Impact of rotor setting and relative arrangement on aerodynamic performance of double rotor arrays

The impact of rotor setting and relative arrangement on the individual and overall power performance and aerodynamics of double rotor vertical axis wind turbine (VAWT) arrays is investigated. Eight rotor settings are considered: two relative rotational directions (co-rotating, CO, and counter-rotating, CN), two relative positionings (downstream turbine positioned in the leeward, LW, and windward, WW, of the upstream rotor), and two phase lags ( ∆ θ = 0 ◦ and 180 ◦ ). For each of the eight rotor settings, 63 different relative arrangements are considered resulting in 504 unique cases. The arrangements are considered within 1.25 d ≤ R ≤ 10 d ( d = rotor diameter) and 0 ◦ ≤ Φ ≤ 90 ◦ , where R and Φ are relative distance and angle of the rotors, respectively. Unsteady Reynolds-Averaged Navier–Stokes (URANS) CFD simulations, validated with experimental data, are employed. The results show that the power performance of the array is significantly influenced by the relative rotational direction and positioning, ∼ 8% in power coefficient ( C P ), while it is marginally dependent on relative phase lag. The different performance of the studied arrays is because of different parts of the downstream turbine revolution being affected by the wake of the upstream turbine and dissimilar strength/width of the shear layer created in the two rotors’ wake overlap. The preferred rotational direction for WW arrays is co-rotating while for LW arrays counter-rotating is favored. For the same arrangement, counter-rotating turbines with LW relative positioning have the highest C P due to their downstream turbine blade moving along the flow direction in the wake overlap region resulting in little energy dissipation and weak shear layer. In contrast, counter-rotating arrays with WW relative positioning have the lowest C P , because the downstream turbine blade moves against the flow in the wake overlap region, resulting in extensive velocity deficit and a thick, strong shear layer.


State-of-the-art and research gaps
Wind farm rotor setting and arrangement are the two key factors that need to be optimized to minimize the wake losses within a wind farm (González-Longatt et al., 2012;Bartl et al., 2012;Cazzaro and Pisinger, 2022) and consequently maximize the overall power output of the farm (Shakoor et al., 2016;Dabiri, 2011;Ti et al., 2021). Wind farm arrangement is defined as the relative distance, angle, and positioning of the rotors. Rotor settings include, but are not limited to, rotational speed and direction and could be individually adjusted for each rotor so that each rotor would operate optimally based on its own local conditions. The local flow field of a rotor is heavily affected by the adjacent rotors. Therefore, the rotor settings of adjacent turbines are mutually dependent.
The effect of wind farm rotor setting (Montoya et al., 2014;Yang et al., 2019;González et al., 2015) and arrangement (Chowdhury et al., 2012;Kusiak and Song, 2010;Chen et al., 2013) on the performance of horizontal axis wind turbines (HAWTs) has been extensively investigated. According to these studies, the larger the relative distance between adjacent rotors within a farm, the higher the overall power performance of the farm will be (Vermeer et al., 2003;Sanderse, 2009;Azlan et al., 2021). A minimum distance of 5−7d (d = rotor diameter) is reported to result in the optimal performance of the HAWTs without negatively affecting the adjacent rotors (Choi et al., 2013;Bai et al., 2013;Behnood et al., 2014). However, this is not always the case for vertical axis wind turbine (VAWT) farms where placing VAWTs in certain arrangements in close proximity can even increase the overall  (Dabiri, 2011;Sahebzadeh et al., 2020;Zanforlin and Nishino, 2016;Peng et al., 2021).
The research on HAWT wind farm control has also made considerable progress, where it was shown that the use of a wide range of strategies such as power de-rating, yaw-based wake steering and turbine repositioning could significantly enhance the overall power production in the farm (Kheirabadi and Nagamune, 2019;Andersson et al., 2021;Johnson and Thomas, 2009). The same direction could potentially be pursued in VAWT farm design research, with even more promising prospects considering the recent findings on closely-packed VAWT arrays (Dabiri, 2011;Sahebzadeh et al., 2020;Hansen et al., 2021). Earlier studies have suggested that closely-packed VAWTs, set in certain optimal layouts, can benefit from mutual synergic interactions and operate as efficiently as an isolated solo VAWT or possibly even slightly better (Dabiri, 2011;Sahebzadeh et al., 2020;Whittlesey et al., 2021). Furthermore, it has been shown that the relative rotational direction of adjacent rotors can further contribute to their power efficiency. For instance, some studies have reported higher power performance for counter-rotating rotor pairs, compared to co-rotating ones (Dabiri, 2011;Zanforlin and Nishino, 2016;Hassanpour and Azadani, 2021). The aforementioned findings indicate a great potential for compact VAWT farms with higher power generation per area of land (power density), compared to HAWT farms; a potential that is achievable through optimal individual rotor setting and relative arrangement in wind farms.
This potential for high power density VAWT farms has led to an increasing interest in the characterization and optimization of VAWT farms (Dabiri, 2011;Sahebzadeh et al., 2020;Shaheen and Abdallah, 2017;Shaaban et al., 2018). Table 1 presents an overview of the studies on the optimal layout design of VAWT farms. The table focuses on studies in which both the rotor setting and relative arrangement (farm layout) of the turbines have been investigated. The table details the method of each study, the number of studied rotors, arrangements (layout), rotor settings, the studied parameters, relative rotational direction(s), and the focus of the study. According to the table, while some experimental studies have been conducted on VAWT farm design, 2D URANS simulations have been widely used. Relative rotational direction, phase lag and distance have been studied in different arrangements with two focus points of power performance and wake aerodynamics. The table also shows that the number of investigated rotor settings and relative arrangements has been rather limited, not sufficient to provide a comprehensive understanding of the impact of these two key factors on the performance of the studied arrays.
The effect of relative arrangement (i.e., layout) of co-rotating double rotor VAWT arrays on the individual and overall power performance of the turbines has been studied in detail (Dabiri, 2011;Zanforlin and Nishino, 2016;Chen et al., 2017). Sahebzadeh et al. (2020) developed a performance map for optimal layout design of such arrays and identified the flow acceleration between the rotors as the main contributing factor to the enhanced power performance. This finding is in line with the numerical results of Alexander and Santhanakrishnan (2019) and the wind-tunnel measurements of Ahmadi-Baloutaki et al. (2016).
Further investigation of the literature points to a potential in relative rotational direction (RD), phase lag (∆θ ) and positioning (RP) (i.e., the position of the downstream turbine in respect to the upstream rotor in the arrangement) of VAWTs for increased power performance. In this regard, rotors can be converging, i.e., blade(s) pushing the flow downstream, or diverging, i.e., blade(s) pushing the flow upstream. However, the conclusions of previous studies regarding the impact of relative rotational direction, phase lag and positioning are not always consistent. For example, Shaheen and Abdallah (2017) conducted a study on double rotor VAWT arrays and used their findings to design a triple rotor VAWT array which showed a 30% power performance increase compared to an isolated rotor. This is while Shaaban et al. (2018) stated that the average power coefficient of closely spaced VAWTs cannot surpass that of an isolated solo  CN Peng et al. (2022) CFD (3D LES) 2 R, RD, σ , PA, AF 4 (4) C P , Wake CO, CN FM = field measurement, CFD = computational fluid dynamics, URANS = unsteady Reynolds-Averaged Navier-Stokes, WT = wind-tunnel measurement, DES = detached eddy simulation, LES = large eddy simulation, R = relative distance, Φ = relative angle, RD = rotational direction, WD = wind direction, λ = tip speed ratio, RP = relative positioning, ∆ θ = phase lag, AF = airfoil shape, PA = pitch angle, σ = solidity, C P = power coefficient, C T = thrust coefficient, PD = power density, CO = co-rotating, CN = counter-rotating.
rotor. According to Chen et al. (2017), the effect of the relative rotational direction on the performance of VAWT arrays is more significant than that of the relative distance. This was confirmed by the field measurements of Dabiri (2011), wind-tunnel measurements of Vergaerde et al. (2020a,b,c) and numerical studies of Zanforlin and Nishino (2016), Bangga et al. (2018), Jiang et al. (2020a) and Jiang et al. (2020b) according to which converging counter-rotating VAWTs result in less energy dissipation between the rotors compared to co-rotating arrays, resulting in a flow with higher mean energy moving downstream. However, these findings were contradicted by Peng et al. (2020) and Peng et al. (2022), who concluded that rotational direction has no significant impact on the power performance of VAWT arrays. Zanforlin and Nishino (2016) and Chen et al. (2017) studied the impact of wind direction on the individual and overall power performance of the turbines in double rotor VAWT arrays for only a few wind directions. This literature review indicates the need for further research to clarify the potential of relative rotational direction, phase lag and positioning in increasing the power performance of VAWT arrays and to comprehensively understand the underlying physics. This is due to the following gaps: 1. The majority of studies have considered a limited number of relative arrangements, mainly focusing on side-by-side rotors. This leaves our knowledge of the generalized dependency of the arrays' performance on relative distance and angle deficient. As a result, the impact of relative rotational direction, phase lag, and positioning on the power performance of VAWT arrays is not comprehensively quantified; 2. The impacts of physical mechanisms associated with relative rotational direction, phase lag and positioning on the power performance of VAWT arrays are not yet fully understood; 3. There is little consensus in the literature on the optimal rotor setting in different arrangements and in some cases, the findings of different studies are contradictory. As a result, an overall map to indicate the optimal rotor setting as a function of relative arrangement cannot be developed.

Objectives and novelties
This study seeks to address the above-mentioned gaps by investigating the power performance and aerodynamics of double rotor VAWT arrays, as the smallest generating cell in a wind farm, in different arrangements with different rotor settings. The objectives of the work are (i) to understand the impact of relative rotational direction, phase lag, and positioning of the turbines on the power performance of the individual rotors in a double rotor VAWT array as well as the overall array; (ii) to identify the associated underlying aerodynamic mechanisms; and (iii) to develop high-resolution power performance maps for different rotor settings and arrangements, which is needed to identify the optimal cases.
To realize the aforementioned objectives, extensive highfidelity URANS simulations, validated with experimental data, are carried out to study the impact of relative rotational direction (RD), phase lag (∆θ ) and positioning (RP) on the individual and overall power performance and aerodynamics of double rotor Darrieus H-type vertical axis wind turbine arrays in a wide range of relative distances (R) and angles (Φ). In total, 504 unique cases Tip speed ratio (based on U ∞ ), λ [-] 4 Chord-based Reynolds number, Re c [-] 1.57 × 10 5 are investigated in 8 different rotor settings and relative rotations and 63 different arrangements.
The important impact of relative arrangement on double-rotor VAWTs was shown in Sahebzadeh et al. (2020). However, to the best knowledge of the authors, the interlinked impact of rotor setting and relative arrangement on the aerodynamic performance of double-rotor VAWT arrays has not yet been systematically investigated. Therefore, this study investigates the impact of rotor settings and the concurrent impact of rotor setting and relative arrangement on the aerodynamic performance of VAWT arrays by considering a wide range of rotor settings (RD, ∆θ and RP) in a wide range of relative arrangements (R and Φ). The findings will better clarify the blade, near-wake, and far-wake aerodynamic mechanisms leading to increased/decreased power performance of double-rotor VAWT arrays. Furthermore, this study will develop a map of the optimal rotor setting corresponding to the highest power performance for different double-rotor VAWT arrangements. To the best knowledge of the authors, this map does not yet exist in the literature and can be a useful tool in designing optimal wind farms.

Outline
The rest of the paper is organized as follows: Sections 2.1-2.3 detail the geometrical and operational characteristics of the turbines, computational domain and grid, boundary conditions and other computational settings, respectively. This is followed by solution verification and validation in Section 2.4. Test cases are introduced in Section 3. The power performance of the upstream turbine, downstream turbine and the overall array is investigated in Section 4. Blade aerodynamics and turbine wake are analyzed in Section 5. Discussion and conclusions are provided in Sections 6 and 7. Table 2 details the geometrical and operational characteristics of the studied turbines. These characteristics are selected with respect to the validation study, which will be presented in Section 2.4. The rotors are considered with only one blade in order to reduce the computational cost of a large number of transient simulations. It should be noted that according to the findings of Rezaeiha et al. (2018c), the number of blades has a marginal effect on wake characteristics and aerodynamic performance of low solidity VAWTs operating within their optimal regime. Less-aerodynamic bodies, i.e., the shaft and connecting rods, are excluded from the rotors. This is in line with the findings of Rezaeiha et al. (2017c) and Tummala et al. (2016), which reported a systematic drop in turbine power performance due to ignoring these bodies. Therefore, the aforementioned simplifications are expected to have a marginal effect on the conclusions of this study.

Geometrical and operational characteristics of the turbines
The rotation direction of the upstream rotor (TurbI) in all cases is counter-clockwise, while for the downstream rotor (TurbII) it can be clockwise or counter-clockwise, depending on the relative rotational direction (RD) of the case.
A solo rotor with the same operational and geometrical conditions as the turbines in the double-rotor array (Table 2) is also simulated to be used as a comparison baseline.

Computational domain and grid
A two-dimensional computational domain is considered to reduce the computational costs of the large number of transient simulations. The computational domain consists of a fixed surrounding domain and two rotating cores. Note that the comparison of the results of 2D and 2.5D computational domains for URANS simulations of VAWTs show a systematic difference (e.g., Rezaeiha et al., 2017aRezaeiha et al., , 2018b. Following the best-practice guidelines for VAWT CFD simulations (Rezaeiha et al., 2017a(Rezaeiha et al., , 2018b, a 35d × 40d (width × length) computational domain is developed resulting in a 2D blockage ratio of ( 2d W ) ≤ 5% for the double rotor array (Fig. 1). In all simulations, the minimum distance between the upstream turbine center and the domain inlet is considered to be 15d. This minimizes the effect of the domain boundaries on the turbines' upstream induction field. The minimum distance between the downstream turbine center and the domain outlet is considered to be 10d allowing the wake to fully develop in the stream-wise direction before reaching the outlet. For all simulations, an equal distance (d s ) is considered between the centers of the two rotors and the lateral boundaries. In order to allow the wake to fully develop in the lateral direction, this distance is always equal or larger than 10d (d s ≥ 10d).
High-quality unstructured computational grids are developed, containing 0.7 to 1.4 million cells, depending on the arrangement of the array. A body of influence is considered in the vicinity of the rotors, extending 15d downstream of each rotor, in which the computational grid is further refined to better capture the nearand far-wake interactions. Each blade is discretized by 800 cells along its circumference. For all the studied arrays, the boundary layer of the airfoils is discretized by 20 cell layers with a growth rate of 1.1. This results in maximum and average y + values of 4.2 and 1.8, respectively. Fig. 2 depicts the computational grid for a sample case.
A domain with the same 5% blockage ratio as the double rotor arrangements is developed for the solo rotor. In this case, about 0.5 million cells with the same topology as the double rotor cases are used. This results in the average and maximum y + values of 1.7 and 4.1, respectively.

Boundary conditions and other computational settings
Incompressible unsteady Reynolds-averaged Navier-Stokes (URANS) equations are solved employing the commercial CFD code ANSYS Fluent v19.1 (ANSYS, 2020). The boundary conditions and the rest of the computational settings are presented in Table 3. The sliding grid technique (Steijl and Barakos, 2008) is used for the interface between the two rotating cores and the fixed domain. The computational settings are set according to the best-practice guidelines for accurate CFD simulation of vertical axis wind turbines (Rezaeiha et al., 2017a(Rezaeiha et al., , 2018b(Rezaeiha et al., , 2019b. The results of steady RANS simulations (with the blade of both rotors fixed at θ TurbI = 0 • ) are used to initialize the transient simulations.
The selected 4-equation transition SST (γ -Re θ ) turbulence model (Menter et al., 2004) accounts for laminar-to-turbulent transition by considering two extra transport equations for momentum-thickness Reynolds number (Re θ ) and intermittency (γ ), in addition to the two transport equations of SST k-ω turbulence model (D.C., 1998). The choice of this turbulence model is based on the results of the validation studies (Section 2.4) and the best practice guideline for the accuracy of turbulence models for high-fidelity CFD simulations of VAWTs according to which the performance of this model is found to be superior to other eddy-viscosity turbulence models (Rezaeiha et al., 2019b) when compared against the more complex scale-resolving VAWT simulations (Rezaeiha et al., 2019a) and the experimental measurements. In order to limit the turbulence production in the stagnation regions, the production limiters developed by Menter (1994) and Kato (1993) are employed. The curvature correction developed by Smirnov and Menter (2009) is used as a modification to the turbulence production term to account for the effects of system rotation and streamline curvature on the turbulence model.
All simulations are performed on 24-core CPU nodes (Xeon E5-2690v3, 2.6 GHz) with 64 GB of memory per node, accounting for an average of ≈2800 core hours per simulation.

Solution verification and validation
A comprehensive grid-sensitivity study is carried out in which three uniformly-refined fine, medium and coarse grids of a sample double rotor array are investigated to ensure the grid independency of the results. The number of cells for the fine, medium and coarse grids are 1,381,353, 745,422 and 414,096, respectively. For the grids, the maximum y + values are 2.7, 3.8 and 5.7, respectively. The grid convergence index (GCI) (Roache, 1997) for the three grids is calculated using the overall power coefficient of the array with a safety factor of F s = 1.25 following the best practice guidelines. For the medium-fine grid pairs, the GCI coarse and GCI fine values are 3.2 × 10 −3 and 2.4 × 10 −3 , respectively. These values correspond to 1.6% and 1.2% of the Richardson extrapolated overall power coefficient. Furthermore, the differences between the stream-wise and lateral velocity profiles in the wake of the array in the three grids are negligible with a maximum and average deviation of 1.5% and 0.2% between the coarse and medium grids for stream-wise velocity, respectively. This is 0.8% and 0.1% between the medium and fine grids. The maximum and average deviation of lateral velocity between the coarse and medium grids are 0.15% and 0.07%. This is 0.05% and 0.03% between the medium and fine grids. Accordingly, the medium grid is selected for the remainder of the research. More detailed information about the grid-sensitivity analysis is extensively presented in Sahebzadeh et al. (2020).
The CFD model is validated using three separate validation studies in which the CFD results are compared against three separate sets of wind-tunnel measurements on one-bladed twobladed and three-bladed Darrieus H-type vertical axis wind turbines conducted by Simão Ferreira et al.   (Menter et al., 2004) Temporal discretization order 4-stage Runge-Kutta scheme (Gottlieb and Shu, 1998) Spatial discretization order Second-order upwind scheme (Barth and Jespersen, 1989) Pressure-velocity coupling scheme SIMPLE (Patankar and Spalding, 1972;Patankar, 1980) Azimuthal increment  the CFD results and the experimental data for azimuthal angles of θ = 108 • , 133 • , and 223 • are 4.5%, 7.4%, and 11.6%, respectively. The highest deviation is observed at θ = 90 • . The CFD results also show a shift towards smaller azimuthal angles. This may be attributed to the different methods available for calculating circulation strength in CFD and experiments where differences in the size of the integration windows could be a source of deviation. ii. Validation study for two-bladed turbine: In the wind-tunnel measurement by Tescione et al. (2014), the stream-wise and lateral velocity components in the near wake of a two-bladed turbine operating in an optimal regime were measured. For this validation study, CFD simulations are performed for the two-bladed turbine used in the experiment. Fig. 3b compares the CFD results and experimental data for stream-wise velocity along the lateral direction. For these two parameters, the deviation between CFD and wind tunnel is 16% and 2.8%, respectively (Fig. 3b).
It should be noted that these deviations are asymmetric, where more significant deviations are observed on the windward side of the rotor. This is believed to be partly due to the higher possibility of flow separation in this region than the leeward side, which introduces more difficulties in modeling by 2D URANS. The higher possibility of flow separation is because, in this range of y/d, the blade moves against the flow as opposed to the leeward side in which the blade moves along the flow. This also applies to the less aerodynamic connecting struts when moving in the windward region, resulting in larger separations compared to the leeward side. Since these connecting struts are present in the experiment and excluded from the CFD simulations (see Section 2.1), this is also believed to have contributed to the deviations. iii. Validation study for three-bladed turbine: CFD results of the power coefficient of a three-bladed turbine are compared against that of a three-bladed turbine in the windtunnel measurements by Raciti Castelli et al. (2011). The deviation between the CFD results and experimental data for λ = 2.51, 2.04, and 2.64 are 0.5%, 3.4%, and 6.8%, respectively (Fig. 3c). The largest deviation is observed at λ = 3.08. This is in line with the findings of Rezaeiha et al. (2017c), according to which the impact of the lessaerodynamic structural components on the C P and the considerable drag that they can generate is more pronounced in higher tip speed ratios. While these components are part of the experiment, they are excluded from the CFD simulation (see Section 2.1). Furthermore, measurement uncertainties are not provided in Raciti Castelli et al. (2011).
Considering the large blockage ratio of ≈10% in this experiment, this can lead to large deviations between CFD and experiment.
In addition to the items mentioned above, the following reasons could also contribute to the deviations between CFD and

Test cases
The studied parameters are as follows: • Relative rotational direction (RD): Co-rotating in which both turbines rotate in the counter-clockwise direction (CO) or counter-rotating in which TurbI rotates in the counterclockwise direction while TurbII rotates in the clockwise direction (CN).
• Relative phase lag (∆θ ): The difference between the two turbines' azimuthal angles (θ TurbIθ TurbII ). • Relative angle (Φ): The angle between the line connecting the two turbines' centers and the X -axis.
The analysis is performed on an array of two Darrieus H-type VAWTs in 504 different test cases. This consists of 8 unique rotor settings each in 63 unique relative arrangements (see Table 4).
According to previous findings, the effect of adjacent VAWTs on one another in R > 10d distances is marginal (Sahebzadeh et al., 2020;Rezaeiha et al., 2018b). Therefore, R = 10d is selected as the largest relative distance to be studied. Relative phase lags of ∆θ = 0 • and 180 • are selected to represent the smallest and largest possible phase lags between the two rotors.
For convenience, a naming scheme is defined where each name consists of the following three parts in order: (i) the relative rotational direction of the two rotors, (ii) the relative positioning, and (iii) the phase lag. For example, CN-LW-180 indicates the cases with a counter-rotating relative rotational direction in which TurbII is positioned in the leeward side of TurbI and there is a 180 • phase lag between the two rotors.
The combination of 3 studied rotor setting parameters, each with 2 distinct values (Table 4), results in 8 unique rotor settings (2RD × 2∆θ × 2RP = 8 rotor settings). Fig. 4 shows these 8 unique rotor settings, i.e., combinations of relative rotational directions, phase lags and positionings. Fig. 5 shows a schematic of the studied relative arrangement parameters in cases with WW and LW relative positionings. Fig. 6 illustrates a schematic of all the studied arrangements in WW and LW relative positionings.

Impact of rotor settings and relative arrangement on power performance
In this section, the power coefficient of the turbines (C P ) is employed to evaluate the impact of rotor setting and relative arrangement on the power performance of each individual rotor as well as the overall array. In Section 5, a comprehensive aerodynamic analysis is carried out to understand the underlying physics causing this impact.
The results in the following sections are normalized by their counterpart value for the solo rotor.
The individual power coefficient of the upstream turbine, C TurbI P , and downstream turbine, C TurbII P , are calculated using Eq. (1): where M, Ω, q, U ∞ and A are moment, rotational speed, dynamic pressure, freestream velocity and turbine swept area, respectively.
The overall power coefficient of the double rotor array, C Overall P , is calculated as the arithmetic mean of C TurbI P and C TurbII P using Eq. (2). (2) For a comparative analysis, all test cases are compared against a reference counterpart arrangement with two co-rotating rotors in WW relative positioning and ∆θ = 0 • phase lag (i.e., CO-WW-0), in the same relative arrangement (i.e., relative distance and angle). A new variable termed as power coefficient index is defined as the difference between the C P of the respective case and its counterpart reference case (with CO-WW-0 rotor setting and the same relative arrangement), normalized by the C P of the isolated solo rotor (C Solo P = 0.287), see Eqs.
(3)-(5).   Fig. 7a and b show that C TurbI P is influenced by the induction field of TurbII. According to the results, the extent of this influence depends heavily on the rotor setting and the relative arrangement of the two turbines, which are elaborated as follows.

Impact of relative rotational direction (RD):
• As shown in Fig. 7a and b, in a WW relative positioning, co-rotating turbines result in higher C TurbI P values while in an LW relative positioning, counter-rotating turbines are preferred.
• The highest difference between co-and counter-rotating turbines in WW arrays is obtained for the R/d = 1.25 and Φ = 30 • arrangement in which the counter-rotating array (CN-WW-0) has 3.6% lower ∆C TurbI P compared to the co-rotating array (CO-WW-0).
• The highest difference between co-and counter-rotating turbines in LW arrays corresponds to the same R/d = 1.25 and Φ = 30 • arrangement, in which the ∆C TurbI P of the counter-rotating array (CN-LW-0) is 4.1% higher than that of the co-rotating array (CO-LW-0) ( Fig. 7a and b, Fig. 8a-d and Fig. 9a-d).
• These differences are because for 10 • < Φ < 30 • , the largest portion of TurbI revolution is located in TurbII induction field. For Φ > 30 • or R/d > 1.25, the ∆C TurbI P difference between the 8 studied rotor settings (Fig. 4) decreases considerably since a smaller part of TurbI revolution is located in TurbII induction field, while for Φ < 10 • , the part of TurbI revolution located in TurbII induction field is not substantially different for different rotor settings. • The effect of relative phase lag is negligible for co-rotating turbines, while slightly more apparent for the counterrotating cases. The maximum difference of ∆C TurbI P between turbines with ∆θ = 0 • and ∆θ = 180 • is observed for the CN-WW-0 and CN-WW-180 arrays in R/d = 1.25 and Φ = 15 • . In this arrangement, the ∆C TurbI P of the array with ∆θ = 180 • is about 2% higher than the array with ∆θ = 0 • . For all other rotor settings and relative arrangements, the impact of relative phase lag on ∆C TurbI P is less than 1%.
• For the counter-rotating turbines, the impact of relative phase lag on ∆C TurbI P highly depends on Φ (see Fig. 8a-d and Fig. 9a-d). For 0 • < Φ ≤ 45 • , the CN-LW-0 arrays have a higher C TurbI P than the CN-LW-180 arrays, and the CN-WW-180 arrays have a higher C TurbI P than the CN-WW-0 arrays. This trend is reversed for 45 • < Φ ≤ 90 • .

Impact of relative positioning (RP):
• For co-rotating turbines, a WW relative positioning results in C TurbI P enhancement, while the opposite applies to the counter-rotating turbines for which an LW relative positioning results in C TurbI P increase.
• For co-rotating turbines, the highest difference between LW and WW arrays is reported for R/d = 1.25 and Φ = 45 • , where the ∆C TurbI P for the CO-LW-180 array is about 2.5% less than the CO-WW-180 array. • For counter-rotating turbines, the highest difference between LW and WW arrays occurs for R/d = 1.25 and Φ = 30 • where the ∆C TurbI P of the CN-LW-0 array is 5.7% higher than that of the CN-WWW-0 array ( Fig. 7a and b, Fig. 8a-d and Fig. 9a-d).

Impact of relative arrangement (R and Φ):
• Increasing the values of R and Φ decreases the impact of relative rotational direction and positioning on ∆C TurbI P . For relative distances of R ≥ 3d or Φ ≥ 75 • , the ∆C TurbI P is below 1%.
• The impact of relative phase lag on ∆C TurbI P remains marginal regardless of the relative arrangement.   Fig. 4) and relative arrangements, respectively.

Impact of relative rotational direction (RD):
• For a WW relative positioning, selecting co-rotating turbines generally leads to higher C TurbII P values while in an LW relative positioning, counter-rotating turbines are generally superior. Note that the same trend is observed for TurbI. • The highest impact of relative rotational direction on LW arrays is observed for R/d = 1.5 and Φ = 15 • arrangement, where counter-rotating turbines (CN-LW-180) have a 6.8% higher ∆C TurbII P compared to co-rotating turbines (CO-LW-180) ( Fig. 7c and d, Fig. 8e-h and Fig. 9e-h).
• For a given R/d, the highest impact of relative rotational direction on ∆C TurbII P is observed at Φ = 15 • . This is because this relative angle results in the largest portion of TurbII revolution being located in the accelerated flow region between the two rotors, where depending on the relative rotational direction, C TurbII P may increase or decrease substantially. This is further investigated in Section 5.    for counter-rotating turbines (CN-WW-0) is about 4.1% less than co-rotating turbines (CO-WW-0).
• The highest impact of relative rotational direction on LW arrays occurs at R/d = 1.5 and Φ = 15 • arrangement, in which the ∆C TurbI P of counter-rotating turbines (CN-LW-0) is 3.8% higher than that of the co-rotating turbines (CO-LW-0) ( Fig. 7e and f, Fig. 8i-l and Fig. 9i-l).
• These differences are because a Φ = 15 • relative angle results in the largest portion of TurbII revolution being located in the accelerated flow region between the two rotors, where depending on the relative rotational direction, C Overall P may increase or decrease substantially. This is further investigated in Section 5.

Impact of relative phase lag (∆θ ):
• The aforementioned marginal impact of relative phase lag on the power performance of both the upstream and downstream turbines translates to its marginal impact on the overall power performance of the array. Therefore, the ∆C Overall P contour plots in Fig. 7e-f are only presented for ∆θ = 0 • .

Impact of relative positioning (RP):
• The impact of relative positioning on ∆C Overall P is similar to that on ∆C TurbII P , meaning that an LW relative positioning increases the overall power performance of the array compared to a WW relative positioning. This increase is more extensive for counter-rotating cases.

Impact of relative arrangement (R and Φ):
• Similar to the observations made for TurbI and TurbII, increasing R and Φ decreases the impact of relative rotational  (Fig. 4), three distinct regions, i.e., set of relative arrangements, can be recognized in the R − Φ  • Fig. 10 shows that there exists a common wake as well as a common optimal region that are shared between all the studied rotor settings. Furthermore, the extents of the wake, optimal and minimal interaction regions are shown to be dependent on the rotor setting.
• Accordingly, Φ = 75 • and 1.25 ≤ R/d ≤ 2.25 range is the shared optimal region between all studied rotor settings, resulting in overall power performance increase regardless of the relative rotational direction, phase lag, or positioning.
• The optimal region extends farther than the aforementioned common R − Φ range depending on the rotor setting. The • The maximum overall power performance in all the studied cases is C Overall P / C Solo P = 1.021, reported at R/d = 1.25 and Φ = 75 • for a CO-WW-180 array.
The aforementioned observations are consistent with the wind-tunnel measurements of Brownstein et al. (2019) and Müller et al. (2021), who have also reported the impact of rotor settings and relative arrangement on the individual and overall power performance of the array. Fig. 11 presents the optimal relative rotational direction and positioning, resulting in the highest C Overall P for all the studied arrangements. It can be seen that if TurbII is placed on the leeward side of TurbI, counter-rotating turbines are preferred. In contrast, TurbII with a co-rotating relative rotational direction on the windward side of TurbI results in higher overall power performance. The figure shows the individual and interlinked impact of relative rotational direction and positioning on the power performance of the farm. It should be noted that relative phase lag has a negligible effect on the power performance of the turbines.

Blade aerodynamics and turbine wake
A comprehensive aerodynamic analysis is carried out to understand the underlying physics corresponding to the impact of relative rotational direction and positioning on the power performance of the double rotor VAWT arrays. The aerodynamic analysis is performed for (i) blade, (ii) near-wake and (iii) far-wake. The instantaneous momentum coefficient (C m ), the contribution  of the four revolution quartiles to the generated power compared to that of the solo rotor, the contribution of each quartile to the generated power, and the instantaneous stream-wise velocity deficit are investigated for the blades. In addition, the distribution of stream-wise velocity in upstream, between and downstream of the array, time-averaged and instantaneous values of velocity deficit and normalized Z-vorticity, as well as time-averaged turbulent kinetic energy are investigated in the near-wake. Timeaveraged stream-wise velocity, stream-wise velocity deficit and turbulence intensity are investigated in the far-wake. The analysis focuses on the following four selected cases: CO-WW-0, CO-LW-0, CN-WW-0 (diverging) and CN-LW-0 (converging), all in the same R/d = 1.5 and Φ = 30 • arrangement. The selected arrangement yields noticeable contrast in C P between the cases.
Note that the analysis is performed only for the ∆θ = 0 • phase lag since phase lag is already shown to have a marginal effect on the performance of the array. The power coefficients of the selected cases are detailed in Table 5. The differences observed in the performance of TurbI are due to the differences in the induction fields of TurbII. As a result, the incoming flow velocity experienced by TurbI varies. This was also reported by earlier studies, e.g., Sahebzadeh et al., 2020;Brownstein et al., 2019;Zanforlin, 2018. Fig . 12a shows that the least and the most experienced velocities by TurbI within −0.5 ≤ y/d ≤ 0.5 (where TurbI is positioned) belong to the CN-WW (diverging) and CN-LW (converging) cases, respectively. This trend is in line with the corresponding C TurbI P values given in Table 5. Fig. 12a also reveals a slight asymmetry in stream-wise velocity along the lateral direction in which greater induction is  Table 5). Note that the light and dark gray backgrounds indicate the y-range where TurbII is positioned in WW and LW cases, respectively.. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) observed at y/d > 0, where TurbII is positioned for WW cases and at y/d < 0, where TurbII is positioned for LW cases. An asymmetry is also present in the lateral velocity component shown in Fig. 12b. However, the lateral velocity asymmetry is dominated by the local blade effects of TurbI airfoil's motion on its immediate surrounding and relative rotational direction and positioning only slightly amplify or diminish the extent of this lateral velocity asymmetry. Note that the extent of the observed asymmetry in the wake depends on the turbine solidity (Rezaeiha et al., 2018c). Fig. 13 depicts the instantaneous momentum coefficient (C m ), the corresponding contribution of the four revolution quartiles to the generated power of the selected cases compared to that of the solo rotor, and the contribution of each quartile to the generated power for TurbI and TurbII. A marginal difference of less than 2% is observed among the C m values of TurbI throughout the upwind (45 • ≤ θ TurbI ≤ 135 • ) and leeward (135 • ≤ θ TurbI ≤ 225 • ) quartiles (see Fig. 13a). The largest variations occur at the downwind (225 • ≤ θ TurbI ≤ 315 • ) and windward (315 where the difference between the quarterly C P contributions of the CN-WW (diverging) and the CN-LW (converging) cases is about 11%. In the downwind and windward quartiles, the CN-WW (diverging) and the CN-LW (converging) cases have the lowest and highest TurbI power performances, respectively ( Fig. 13a and  c), which is consistent with the C P values presented in Table 5. As shown in Fig. 13e, the contribution of the four quartiles to generated power of TurbI is similar in all selected cases with most of the power being generated during the upwind, downwind, leeward and windward quartiles, respectively. The same order was also reported by Rezaeiha et al. (2018a). The CN-WW (diverging) case has the highest upwind quarterly contribution, but the lowest C TurbI P amongst the selected cases, while the CN-LW (converging) case with the lowest upwind quarterly contribution to the generated power, has the highest C TurbI P (Table 5). According to Fig. 13c and e, the CN-LW (converging) case gains its advantage in C TurbI P during the downwind and windward quartiles. Fig. 13  Fig. 13. (a, b) Instantaneous momentum coefficient (C m ), (c, d) corresponding contribution of the four quartiles to the C P compared to counterpart value for the solo rotor and (e, f) contribution of each quartile to power of TurbI and TurbII for selected cases (see Table 5). also shows that the values reported for the CN-LW (converging) case are close to those of the solo rotor throughout its revolution and amongst the selected cases, the CN-LW (converging) case performs the most similar to the solo rotor.
According to Table 5, the CN-WW (diverging) case has the lowest C TurbI P and C TurbII P while the CN-LW (converging) case has the highest power coefficients. To understand the associated physical mechanisms influencing the power performance of the selected cases, near wake and far wake of the turbines are further analyzed, with a focus on TurbII. Table 6 presents the range of azimuthal angles of TurbII (θ TurbII ) located in the direct wake of TurbI (downstream projection of TurbI circumference) and the corresponding quartiles. Fig. 14 shows the dimensionless instantaneous stream-wise velocity deficit experienced by TurbII for the selected cases with blade positioned at the center of the four quartiles. Figs. 15 and 16 depict the time-averaged and instantaneous values of velocity deficit and normalized Z-vorticity, respectively. Fig. 17 presents the time-averaged turbulent kinetic energy contours of the selected cases ( Table 5). The observations are categorized based  Table 5 with blade positioned at the center of (a-d) upwind, (e-h) leeward, (i-l) downwind and (m-p) windward quartiles. Note that all selected cases have a ∆θ = 0 • phase lag. on the revolution quartiles, namely upwind, leeward, downwind, and windward as follows: • Upwind quartile: TurbII in all the selected cases shows comparable performance to the solo rotor (Fig. 13d). The lowest quarterly C P value compared to the solo rotor belongs to the CN-WW (diverging) case with about 95% of the power coefficient of the solo rotor in this quartile. This is due to the greater velocity deficit in the windward side of the TurbI wake in this case (Fig. 12c, Table 6, Fig. 14a-d and Fig. 15). This is also observed by previous numerical (Alexander and Santhanakrishnan, 2019;Duraisamy and Lakshminarayan, 2014) and experimental studies (Vergaerde et al., 2020a,b;Rolin and Porté-Agel, 2018). This lateral asymmetry in the velocity deficit profile of TurbI can be explained by the different angle of attacks and wind speeds that a blade experiences when moving against the wind (in the windward side of the rotor) as opposed to when moving along the wind (in the leeward side of the rotor).
• Leeward quartile: the CO-LW and CN-WW (diverging) cases show comparable performance to that of the solo rotor, while the power performances of the CN-LW (converging) and CO-WW cases are 77% and 54% of the solo rotor in this quartile, respectively. This is because during this quartile,  Table 5.
TurbII in CN-LW (converging) and CO-WW cases, moves in the low-velocity wake of TurbI (Fig. 12c, Table 6, Fig. 14e-h and Fig. 15). The different performance of the CO-WW and CN-LW (converging) cases is because of the greater velocity deficit in the windward side of the TurbI wake compared to its leeward side (Fig. 12b, Fig. 14e and h and Fig. 15a, d).
In addition, TurbII in the CO-WW case is passing through shear-layer flow, experiencing stronger blade-wake interactions than the other cases (Fig. 15a, Fig. 16a and Fig. 17a). Less power performance deficit of the CN-LW (converging) compared to the CO-WW case is because it receives a more turbulent incoming flow. Fig. 17d shows that in this quartile, TurbII in the CN-LW (converging) case receives a more turbulent incoming flow due to higher turbulent kinetic energy generated by the upstream rotor. According to the literature, the rotor benefits from these increased turbulence levels at such azimuthal angles, resulting in higher power performance (Rezaeiha et al., 2018a).
• Downwind quartile: the CN-LW (converging) and CN-WW (diverging) cases show comparable performance to that of the solo rotor, while the power performances of the CO-LW and CO-WW cases are about 107% and 95% of the solo rotor in this quartile, respectively. The decreased power performance of the CO-WW case in this quartile is due to its TurbII being partially positioned in the low-velocity wake of the TurbI wake (Figs. 12c and 15). Furthermore, Fig. 16e shows that the CO-WW case experiences the strongest blade-wake interactions compared to other cases, resulting in its lower power performance in this quartile. The enhanced power performance in the downwind quartile of the CO-LW case is due to the deviation of the flow by TurbI towards its leeward side where TurbII of the CO-LW case is positioned.
Consequently, TurbII of the CO-LW experiences less velocity deficit compared to the solo rotor in this quartile (Fig. 12c, Fig. 14j and Fig. 15b and f).
• Windward quartile: A significant difference is observed for the performance of TurbII in the windward quartile where the CO-LW and CN-WW (diverging) cases show 59% and 92% power performance reductions compared to the solo rotor. For the CO-WW and CN-LW (converging) cases, however, a 5% increase is observed. This is especially interesting since the windward quartile contributes the least to the power performance of TurbII, but the difference between the selected cases in this quartile has a significant impact on their C TurbII P . The decreased power performance of the CO-LW and CN-WW cases in this quartile is because for both cases the entirety of TurbII windward quartile is located in the low-velocity wake region of TurbI (Table 6, Figs. 14m-p and 15) and TurbII experiences stronger bladewake interactions compared to the other two cases (Fig. 16). The higher velocity deficit in the windward side of the TurbI wake accounts for the dissimilar performances of the CO-LW and CN-WW cases. In the CO-LW case, TurbII is positioned in the LW side of TurbI and receives a flow with less velocity deficit, while in the CN-WW case, TurbII is positioned in the WW side of TurbI and receives a flow with a higher velocity deficit. This velocity deficit is especially greater in 0.25 < y/d < 0.5 where TurbII in the CN-WW case passes through its windward quartile compared to the −0.5 < y/d < −0.25 where TurbII in the CO-LW case passes through its windward quartile (Fig. 12c).
According to Fig. 13f, similar to TurbI, most of the power of TurbII is generated during the upwind and downwind quartiles, Table 6 Range of θ TurbII located in direct wake of TurbI (downstream projection of TurbI circumference, indicated with gray background in the schematics) and the corresponding quartiles for the selected cases in Table 5 respectively. In these two quartiles, all selected cases have a fairly similar contribution to C TurbII P . However, significant differences are observed at the leeward and windward quartiles, i.e., the quartiles with the least contribution to C TurbII P . Accordingly, even though the windward quartile contributes the least to the overall power production of TurbII (Fig. 13f), the superior performance of the CN-LW (converging) case in this quartile (Fig. 13d), results in CN-LW (converging) case having the highest C TurbII P among the selected cases (Table 5). In contrast, the inferior performance of the CN-WW (diverging) case in the windward quartile (Fig. 13d), results in the CN-WW (diverging) having the lowest C TurbII P among the selected cases (Table 5).
Further flow field analysis inside TurbII rotor plane better clarifies the different aerodynamic mechanisms of the selected cases and the interlinked influence of relative positioning and rotational direction on their performance (Fig. 14).
According to Fig. 14a-d, during the CN-WW (diverging) case's upwind quartile, TurbII experiences the most significant velocity deficit, compared to other cases, passing through a low-velocity flow region in the windward side of the TurbI wake. Fig. 14e-h show that the leeward quartile of TurbII in the CO-WW and CN-LW (converging) cases is located entirely in the wake of TurbI. The most significant velocity deficits are experienced by TurbII in the CO-WW case, followed by the CN-LW (converging) case.
According to Fig. 14i-l, TurbII in the CO-WW case experiences the most significant velocity deficit, compared to other cases, along its downwind quartile, while the CO-LW case experiences the least velocity deficit. Fig. 14m-p show that the windward quartile of TurbII in the CO-LW and CN-WW (diverging) cases is located entirely in the wake of TurbI. The most significant velocity deficits are experienced by TurbII in the CN-WW case, followed by the CO-LW (converging) case.
The observations above are reflected in the quarterly C TurbII P of the selected cases presented in Fig. 13d. Fig. 14 clearly shows that TurbII in the CN-WW (diverging) case experiences the greatest velocity deficits for the longest portion of its revolution, while for the CN-LW (converging) case, the downstream turbine experiences the smallest velocity deficits for the shortest portion of its revolution. This corresponds to the lowest and the highest C TurbII P values being obtained for the CN-WW and CN-LW cases, respectively (Table 5).
As mentioned, Fig. 15 depicts the time-averaged and instantaneous values of velocity deficit for the selected cases in Table 5. A closer look at the figure reveals that in the CN-LW (converging) case, TurbII blade enters the wake overlap region after passing through its upwind quartile and entering its leeward quartile, moving in the same direction as the flow. This results in a small velocity gradient. Consequently, a smooth transition is formed between the two individual wakes where, in contrast to the other selected cases, a comparatively weak shear layer is formed in the overlap wake region, i.e., where the two wakes are almost merged ( Fig. 15d and h). The reverse holds for the CN-WW (diverging) case in which TurbII blade enters the wake overlap region after passing through its downwind quartile and entering its windward quartile, moving against the flow. This results in a large velocity gradient and greater energy dissipation in the wake (Ahmadi-Baloutaki et al., 2016;Vergaerde et al., 2020a,b). Consequently, a steep transition is formed between the two individual wakes, resulting in a thick and strong shear layer ( Fig. 15c and g).
This shear layer is also evident in the stream-wise velocity profiles along lateral lines downstream of the array (Fig. 12e), where the sharpest extremum, i.e., velocity deficit, is observed for the CN-WW (diverging) case. The CN-WW (diverging) case has the thickest shear layer in its wake with the highest velocity deficit magnitude among the selected cases.
As mentioned, Fig. 16 depicts the time-averaged and instantaneous vorticity values for the selected cases in Table 5. The aforementioned differences in the shear layer between the wakes of the two turbines are also reflected in the vorticity field. Fig. 16a-d show a sharp vorticity gradient in the windward edge of the wake of the turbines where the magnitude of vorticity increases and then decreases back to zero within a small width compared to the leeward edge of the wake in which the region of increased vorticity is wider with a smoother gradient. The large coherent vortical structures shed by the blade can explain the difference between the windward and leeward edges of wake. These vortices are visible in Fig. 16e-h. Different vorticity fields within the shear layer of the selected cases are observed (Fig. 16), where the strongest and weakest vorticity gradients in the wake overlap region belong to the CN-WW (diverging) and CN-LW (converging) cases, respectively. This observation can be contributed to the different combinations of down-stroke and up-stroke blade movements in the selected cases (Rezaeiha et al., 2017b;Hand et al., 2017). Accordingly, while the overlap wake region of all selected cases is affected by two vorticity fields of different signs, in the CN-WW (diverging) case, the blade of both rotors are in their up-stroke (angle of attack increasing), experiencing increasing lift force, resulting in strong vortex shedding which in turn, leaves a considerable effect on the wake overlap region shear layer. In contrast, the blade of both rotors in the CN-LW (converging) case enter the wake overlap region in their down-stroke (angle of attack decreasing), experiencing a lift force reduction, resulting in weak vortex shedding, leaving a slight effect on the overlap wake region shear  Table 5. layer. For both counter-rotating cases, one rotor enters the wake overlap region in its down-stroke while the other is in its upstroke, thus, the strength of the vortex shedding and its impact on the wake overlap region shear layer is medial, that is neither as strong as the CN-WW (diverging) case nor as week as the CN-LW (converging) case. This is in line with the power performance of the selected cases reported in Table 5. Fig. 16 also indicates some similarities between the wakes of the studied VAWTs and yawed HAWTs, e.g., Rolin and Porté-Agel (2018), Ryan et al. (2016) and Bastankhah and Porté-Agel (2016). In both wakes, counter-rotating vortex pairs are generated and affect the wake structure moving downstream. Furthermore, in both VAWT and HAWT wakes, deflection of the wake towards the rotor windward is observed. This deflection is dependent on relative positioning and rotational direction (Fig. 12e). Note that due to the relatively lower solidity of the rotors in the current study (Table 2), the wake deflection is not as noticeable compared to other studies with higher solidity rotors such as Vergaerde et al. (2020b) and Rezaeiha et al. (2018c). Fig. 17 depicts the different turbulence levels experienced by TurbII in the selected cases. Increased turbulent kinetic energy (TKE) levels are observed at the leeward side of all rotors, corresponding to the blade down-stroke movement where a lower velocity deficit is observed. This is beneficial for the CO-LW and CN-LW cases (see Table 5) in which, TurbII passes through an increased TKE region generated by TurbI in its leeward side and consequently receives an approaching flow with higher turbulence. The increased TKE is also believed to contribute to the higher C P values of the LW cases, as the literature suggests potential benefits of higher TKE for the turbine power performance in this azimuthal range, see Rezaeiha et al. (2018a). It can be seen that the CN-LW (converging) case results in the least velocity deficit in the wake of the two rotors, rendering this rotor setting as the superior choice for adding subsequent turbines (third or more) to extend the array. This is especially important since the CN-LW (converging) already has the highest overall power performance in the selected cases and adding additional rotors is expected to further increase the overall power performance of this array. In contrast, the CN-WW (diverging) case has the greatest velocity deficit in its wake and, therefore, is not a favorable option for adding additional rotors. This is exacerbated by the fact that the CN-WW (diverging) case already results in the lowest overall power performance between the selected cases (Table 5).
Increased turbulence intensity (TI) levels are observed in the core near wake, immediately downstream of the rotors' circumferences. Farther downstream, a sharp decrease of turbulence within the core wake occurs, indicated with white gaps and weakened streaks of turbulence intensity in Fig. 18f-j. This is because the aforementioned large coherent vortical structures shed by the blade (see Fig. 16), are not a long-lasting property of the wake and dissipate quickly moving downstream. Farther downstream, turbulence increase back due to wake break-up, evolving into a far-wake behavior by contracting towards the core of the wake.  Table 5. Fig. 18. Contour plots of (a-e) dimensionless time-averaged (over one revolution) stream-wise velocity and (f-j) turbulence intensity for selected cases in Table 5 in the far field.
Comparison of the turbulence intensity contours in Fig. 18f-j point to less intensively generated turbulence in the CN-LW (converging) case due to the weaker shear layer developing within its wake than the CN-WW (diverging) case. This is because of greater velocity deficit in the CN-WW (diverging) case, since TurbII blade is moving against the flow in the overlap wake region, causing greater energy dissipation in the wake; while in contrast, in the CN-LW (converging) case, TurbII blade moves in the same direction as the flow within the wake overlap region and hence less energy dissipation occurs (Ahmadi-Baloutaki et al., 2016;Vergaerde et al., 2020a,b). This was also previously observed and discussed for near wake stream-wise velocity deficit (Fig. 13). Fig. 19a and b show the time-averaged stream-wise velocity along the stream-wise lines passing through the center of the two turbines. Fig. 19c shows the wake width for the selected cases (Table 5) in comparison to that of a solo rotor. The wake width is defined as the lateral width of the flow region along the rotor(s) centerline where u/U ∞ ≤ 0.97. It should be noted that the subject of investigation in this far wake analysis is the total wake width of the array, which is nevertheless not equivalent to the wake of a solo rotor. However, the dimensions of the wake of the solo rotor are used as a basis for comparison. For a comparative analysis, the wake widths of the double rotor arrays as well as the solo rotor are non-dimensionalized by their frontal width. That is the distance that the turbine(s) cover(s) in the lateral direction.
According to Fig. 19a and b, the velocity deficit along the stream-wise lines passing through the center of the turbines in the double rotor array is lower than that of the solo rotor up to about 11d downstream of TurbII and about 9d downstream of TurbI. This indicates that the wake of the double rotor arrays contains more energy since these arrays have a greater total momentum compared to the solo rotor.
For all selected cases, the stream-wise velocity deficit shows an increasing stream-wise trend (see 0 < x/d < 3 in Fig. 19a and b). This is due to the expanding width of the wake (Fig. 19c) and concurrent cross-stream transmission of the momentum associated with the shed vertical structures (Fig. 16). Moving farther downstream, the stream-wise velocity deficit starts decreasing (see 3 < x/d in Fig. 19a and b) and the shear layers at the edges of the wake start to destabilize (Fig. 18a-e). This is followed by the collapse of the edge shear layers into the core of the wake due to the positive cross-stream advection and inflow of air with higher momentum from the freestream. Fig. 19c and Fig. 18a-e show that all the selected cases experience a similar, about 15%, length increase in array wake compared to the solo rotor. Therefore, the relative rotational direction and positioning have little effect on the wake length.
According to Fig. 19c, for all the selected cases, a narrower wake per frontal width is formed compared to the solo rotor. This is in line with the findings of Zanforlin and Nishino (2016), Alexander andSanthanakrishnan (2019) andDe Tavernier et al. (2018). This narrower wake is because the flow between the two turbines is constricted from expanding freely. The wake of the solo rotor starts oscillating and losing its coherence about 8d downstream of TurbI. While this is delayed to about 11d downstream distance for the double rotor cases where the internal mixing (see Fig. 16) and momentum transfer with the freestream flow weaken the wake till the freestream velocity is reached all over the wake.
As derived by Betz (2014), the width of the wake is correlated with the energy that is harvested from the flow. The larger the extracted energy, the larger the width of the wake will be. This is also the case here (Fig. 19c), where the CN-LW-0 case with the highest C Overall P (Table 5), has a consistently wider wake (about 4%) compared to the rest of the cases. On the other hand, the CN-WW-0 case with the lowest C Overall P (Table 5), has a con- Fig. 19. Dimensionless time-averaged (over one revolution) stream-wise velocity and stream-wise velocity deficit along the stream-wise line passing through the center of (a) TurbI and (b) TurbII and (c) wake width divided by the frontal width of the array (W Frontal ) for selected cases in Table 5. sistently narrower wake (about 4%) than the rest of the cases (Fig. 19c). This is consistent with the wind-tunnel measurements of Vergaerde et al. (2020b) and Vergaerde et al. (2020c), which also indicated wider wakes for the CN-LW (converging) case compared to CN-WW (diverging) case.  Table 5. Light and dark gray backgrounds indicate y-range where TurbII is positioned in WW and LW cases, respectively.. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 20 depicts the normalized time-averaged (over one revolution) stream-wise velocity, stream-wise velocity deficit and lateral velocity along lateral lines in far wake for the selected cases in Table 5. Fig. 20a-c shows a stream-wise surge of the momentum deficit in the wake, up to a downstream distance of 5d, where the stream-wise velocity decreases consistently. This momentum deficit is dominated by the outward cross stream advection, caused by the wake expansion previously seen in Fig. 19c, according to which the wake width of selected cases expands up to a relative distance of about 8d downstream and does not start shrinking up to a distance of about 10d. Fig. 20d-f depict the consistent decrease of lateral velocity, indicating the reduction of blade effects on the flow by moving farther away from the near wake.

Discussion
2D CFD simulations entail some intrinsic limitations compared to full 3D simulations, such as the lower wake recovery rate due to the lack of vertical turbulent mixing and span-wise mean flow reported in the wind-tunnel measurements of Lam and Peng (2017), Vergaerde et al. (2020a) and Vergaerde et al. (2020b). However, such 3D impacts are thought to be systematic, not affecting the derived conclusions in the current study, which are based on a comparative assessment (Sahebzadeh et al., 2020;Zanforlin and Nishino, 2016). Nevertheless, 3D analysis of the subject in this study can be of interest for future research. For example, considering the highly sheared flow in the wake of the studied double rotor arrays (see Section 5), the analysis of vertical influx using 3D simulations may further contribute to aerodynamic understanding of these arrays.
For this research, Unsteady Reynolds-averaged Navier-Stokes (URANS) equations are solved to simulate the double rotor arrays. However, compared to more complex modeling methods including Large-eddy simulation (LES) and Scale-Adaptive Simulation (SAS), URANS suffers from some limitations mainly in capturing complex blade-wake interactions (Rezaeiha et al., 2019a). Therefore, the use of such more sophisticated scale-resolving turbulence modeling techniques can be of interest for future studies.
This study is based on one-bladed VAWTs. It should be noted that one-bladed VAWTs are reported to experience self-starting problems and imbalanced forces, negatively affecting their lifetime (Wagner and Mathur, 2018). A previously published study (Rezaeiha et al., 2018c) shows that the number of turbine blades has a systematic effect on the observed trends for the power performance of VAWTs. Therefore, the use of one-bladed rotors is expected to have a marginal effect on the conclusions of this study. Nevertheless, a comprehensive study about the impact of the number of the turbine blades on the performance of the double VAWT arrays can be of interest for future.
In this study, the focus is not only on power performance (C P ) of the array but also on the underlying physics contributing to the individual and overall power performance of the turbines. Therefore, simulations are performed based on the comprehensive test matrix detailed in Table 4. For such parametric studies, sampling methods such as the design of experiments (DoE) are useful to reduce computational costs. It should be noted that earlier studies have shown that the use of different DoE methods can result in significant variations in predicted results. While some DOE methods are less suitable for computer experiments (simulations), some methods such as space-filling methods are reported to provide a credible characterization of the processes of interest (Rennen et al., 2010;Johnson et al., 1990;Morris and Mitchell, 1995).

Conclusions
The following conclusions are made: • For counter-rotating arrays (CN), windward relative positioning (WW) results in decreased individual and overall power performance of the turbines, while leeward relative positioning (LW) results in an increase. This is due to the more significant velocity deficit on the windward side of the upstream turbine (TurbI).
• For co-rotating arrays (CO), the individual and overall power performance of the turbines marginally depend on their relative positioning.
• The impact of relative positioning on the individual and overall power performance of the turbines in counterrotating (CN) arrays fades away for relative angles of Φ > 45 • and relative distances of R > 1.75d (d = rotor diameter). • The individual and overall power performance of the VAWTs marginally depend on their relative phase lag (∆θ ).
• Within the relative angle and distance range of 0 • < Φ ≤ 30 • and R < 2.25d, the relative rotational direction has a significant impact on the individual and overall power performance of the turbines, resulting in up to 8.3% variations. Outside of this range, the impact of relative rotational direction is limited to less than 2%.
The blade aerodynamics and turbine wake analysis lead to the following conclusions: • The main reason for the different power performances of the studied arrays is that different parts of their TurbII revolution being affected by the wake of TurbI. For all rotors, the upwind, downwind, leeward and windward quartiles contribute the most to the power generation, in that order. Superior performance of the downstream turbine (TurbII) in converging array's, i.e., counter-rotating with a leeward relative positioning (CN-LW) in the least contributing quartile, i.e., windward, is the reason for its higher overall power performance.
• For the diverging array, i.e., counter-rotating with a windward relative positioning (CN-WW), the blade of TurbII moves against the flow in the wake overlap region. This results in great velocity deficit and energy dissipation and consequently, lower power performance. In contrast, for the converging cases (CN-LW) the blade of TurbII moves along the flow in the wake overlap region, resulting in a smaller velocity deficit and almost no shear layer, and consequently higher power performance.
• For leeward relative positioning (LW) arrays, the downstream rotor benefits from passing through an increased turbulent kinetic energy region generated by TurbI in its leeward side, corresponding to the down-stroke movement of the blade. This contributes to the increased power performance of the LW arrays.
• The diverging (CN-LW) and converging (CN-WW) arrays result in the least and most velocity deficit in the far wake, respectively. Therefore, these arrays can be considered as the most and least promising choice for adding additional rotors downstream, respectively.

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.