Method to Predict Outputs of Two-Dimensional VAWT Rotors by Using Wake Model Mimicking the CFD-Created Flow Field

Abstract

Recently, wind farms consisting of clusters of closely spaced vertical-axis wind turbines (VAWTs) have attracted the interest of many people. In this study, a method using a wake model to predict the flow field and the output power of each rotor in a VAWT cluster is proposed. The method uses the information obtained by the preliminary computational fluid dynamics (CFD) targeting an isolated single two-dimensional (2D) VAWT rotor and a few layouts of the paired 2D rotors. In the method, the resultant rotor and flow conditions are determined so as to satisfy the momentum balance in the main wind direction. The pressure loss of the control volume (CV) is given by an interaction model which modifies the prepared information on a single rotor case and assumes the dependence on the inter-rotor distance and the induced velocity. The interaction model consists of four equations depending on the typical four-type layouts of selected two rotors. To obtain the appropriate circulation of each rotor, the searching range of the circulation is limited according to the distribution of other rotors around the rotor at issue. The method can predict the rotor powers in a 2D-VAWT cluster including a few rotors in an incomparably shorter time than the CFD analysis using a dynamic model.

1. Introduction

Wind power generation is one of the alternatives to fossil energy. The obvious advantage of wind power is almost no pollution to the environment. To realize carbon neutrality, a large amount of renewable energy is expected to be introduced [1,2]. Therefore, wind power generators are becoming increasingly larger in size, and their application to offshore generation is increasing [3,4]. However, the efficiency of a wind farm can be reduced owing to the wake effects [5,6,7]. One of the challenges in the sector of wind power is to exactly evaluate the influence of the wakes, so as to improve the efficiency of a wind farm [8,9,10,11,12,13,14]. Additionally, to maximize the power output from a wind farm, it is necessary to deploy wind turbines with the optimal layout for the wind condition of the planned site [15,16,17,18,19]. The wake control methods using the blade’s pitch [20] or the yaw [21,22] of the propeller-typed wind turbines were studied to optimize the wind farm.

According to the studies by Whittlesey et al. [23] and Dabiri [24], the output per unit land area of a wind farm consisting of small-size vertical-axis wind turbines (VAWTs) with a high aspect ratio (rotor height/diameter), which are closely arranged using a unit of counter-rotating paired rotors, can be much greater than that of a conventional wind farm consisting of large-size horizontal-axis wind turbines (HAWTs), which are deployed with the inter-rotor intervals of several multiples of the rotor diameter (in general, the interval in the dominant wind direction is about 10 times as long as the diameter). Since then, many researchers have been interested in the closely arranged VAWT wind farms and, especially, the interaction effects between two VAWT rotors. For example, Zanforlin and Nishino [25] performed a two-dimensional (2D) computational fluid dynamics (CFD) analysis of a pair of inversely rotating VAWTs to show the greater averaged output than the output of an isolated single VAWT. De Tavernier et al. [26] carried out the 2D-CFD based on the panel vortex method of a closely arranged VAWT pair, each of which had a 10 m rotor radius, to show the effects of the load and rotor spacing on the paired rotor performance. Bangga et al. [27] proposed two layouts of a VAWT array based on their CFD study of rotor pairs arranged side by side. Sahebzadeh et al. [28] numerically analyzed the output performance of a co-rotating rotor pair by widely changing the rotor spacing and the relative angle to the main stream. Peng et al. [29] investigated the effects of configuration parameters such as airfoil section, solidity, pitch angle, rotational direction, and turbine spacing on twin VAWTs by CFD analysis. The effects of a three-rotor cluster of VAWTs were also numerically studied by Hezaveh et al. [30] and Silva and Danao [31].

As examples of experimental studies of paired VAWTs, Vergaerde et al. [32] conducted the wind tunnel test using two H-type Darrieus rotors (rotor diameter: 0.5 m, rotor height: 0.8 m [33]). Their turbines were placed side-by-side against the main flow and were adequately controlled by DC motors. They observed the power increase up to 16% for the counter-rotating VAWTs and reported the stable synchronized operation of twin rotors. Jodai and Hara [34] studied the interaction between two closely spaced VAWTs by using miniature 3D printed rotors (diameter: 50 mm, low aspect ratio of 0.87, and high solidity of 0.382) arranged side-by-side. Their experimental results showed a maximum 15% increase in power in the case of the counter-down layout when the inter-rotor space became the shortest (gap space: 10% of the rotor diameter).

Hara et al. [35] applied the dynamic fluid body interaction (DFBI) model to the CFD analysis to simulate the closely arranged paired 2D rotors corresponding to the equator-cross section of the experimental model used in Ref. [34]. Their CFD analysis considering the time-varying rotor speed, for the first time, simulated the synchronization operation of twin rotors and showed the alternation in the angular velocities of two rotors. Furukawa et al. [36] developed an analytical model considering the pressure fluctuation (or increase in flow velocity) observed in the gap region between twin rotors in the above experiments and CFD analyses. The model successfully demonstrated the alternation in the angular velocities of two rotors and showed that the period of the variation in rotor speed depended on the strength of the interaction between the two rotors.

Although the increase in the averaged power of a closely spaced side-by-side VAWT pair is clear, the effects of the distribution of wind direction on the VAWT cluster consisting of many rotors must be investigated more extensively to search for the optimal layout of VAWTs. The CFD analysis, especially the simulation using the DFBI model, can give reliable results, but it needs a long calculation time. Although the experiments can also give useful information, the cost is high and the time for preparation is long. If the number of rotors in a target wind farm increases, the simulation by CFD or the experiment using a lot of rotors is non-realistic. Therefore, a method that can simulate precisely and in a short time the flow field of a wind farm including a large number of VAWTs is necessary.

Buranarote et al. [37] proposed a wake model of a 2D-VAWT rotor, which was named the ultra-super-Gaussian function because the model improved the super-Gaussian function proposed by Shapiro et al. [38] which could express the transformation of the wake profile from top-hat shape to Gaussian. The ultra-super-Gaussian function includes a correction function to express the acceleration regions and the deflection of a VAWT wake. The method proposed by Buranarote et al. [37] was based on the potential flow and included the velocity deficit artificially, like the method by Whittlesey et al. [23]. Buranarote et al. [39] improved their method by including the modification of the y-component (cross-flow) of the flow velocity to mimic the CFD results. Moreover, they introduced the Biot-Savart law to consider the effects of the interaction between the rotors [40] on the wake shift and width. However, in the previous method of Buranarote et al. [39,40], the circulation around each rotor, which was used to estimate the power output, was newly calculated with the modified flow field after adding the velocity deficits; the calculated circulation was sensitive to the flow condition resulted in failure in the reliable prediction of the power outputs of rotors.

This study proposes a new method to predict outputs of two-dimensional VAWT rotors by using a wake model mimicking the CFD-created flow field. The new method is based on the previous method but the circulation that is used to estimate the output of each rotor is the same as the input value to calculate the potential flow in a wind farm. The decision on the appropriate flow and rotor conditions is conducted by evaluating the momentum balance, which is calculated using the momentum transports and the pressure forces at the boundaries of the control volume (CV) and thrust forces of rotors. The necessary CFD data, or the available and reliable experimental data, are the power performance of an isolated single VAWT and the averaged flow velocity distribution in the CV, and the pressure distributions at the boundaries under several wind speed conditions. In addition, the power output data of closely spaced paired VAWTs in typical four layouts in the case of a specific inter-rotor distance are necessary. The interaction effects are considered by modifying the given pressure loss of the isolated single rotor, according to the relative layout of the selected two rotors and considering the distance and the induced velocity.

The method will be validated with the CFD results for two or three rotors studied by Hara et al. [35] and Okinaga et al. [41], in which the rotor height was considered so as to correspond to the experimental rotor used in the experimental study by Jodai et al. [34]. The 2D-CFD analysis of an isolated single rotor is outlined in Appendix B. The same 2D-CFD rotor model is used for the CFD analysis of four-rotor arrays conducted in this study. The present method does not include the three-dimensional effects caused by the finite rotor height because our target in the future is a wind farm consisting of small-scale VAWTs of 14 m diameter with a low aspect ratio.

The final goal of our project is to provide a cost-effective and relatively short-time method to optimize the layout of VAWTs in an arbitrary wind farm. In this paper, at the early stage of the project, we show the possibility to predict a reasonable condition of a VAWT cluster. Therefore, the round robin, which needs a long calculation time when the number of rotors is large, is utilized in the search for adequate conditions. The maximum number of rotors in a VAWT cluster considered in this study is four due to the problem of calculation time. However, if some advanced optimization method like the genetic algorithm (GA) is adopted, the problem of the computation time will be mitigated.

The detail of the new method is described in the next section. The application of the method to an isolated single rotor, paired rotors, three-rotor clusters, and four-rotor layouts is discussed in Section 3. Finally, Section 4 concludes the discussion.

2. Model

2.1. Method

Figure 1 shows a schematic image of a VAWT wind farm, where a VAWT rotor is shown by a circle and our method does not need detailed information on the configuration of a turbine such as the number of blades and the cross-section. In our method, the VAWTs are dealt with as 2D rotors with each having a diameter, D. Let us assume the wind farm consists of N VAWTs. In Figure 1, the coordinate axis x is defined as parallel to the upstream wind speed U_∞; the coordinate axis y is perpendicular to the dominant wind direction. The center position (rotational axis) of the k-th rotor is expressed as (x_k, y_k). Our method assumes that the circulation around an isolated single rotor (Γ) is given as a linear function of wind speed (U_∞) in advance. Therefore, the blockage effect (dipole μ), thrust force (Th), and power output (P) are given as functions of the circulation, respectively, for each rotor.

If an arbitrary set of the conditions of circulation (Γ_k) and blockage effect (μ_k) is given for a wind farm, the complex velocity potential W(z), where z = x + iy is expressed by Equation (1) [23].

$$W\left(z\right)=U_{\infty }z+{\displaystyle \sum }_{k=1}^N\left\{-i\frac{{\Gamma }_k}{2\pi }\ln \left(z-z_k\right)+{\mu }_k{\left(z-z_k\right)}^{-1}\right\}$$

(1)

Using the above potential W(z), the potential flow (u_p, v_p) at an arbitrary position P:(x, y) in the wind farm can be calculated as follows:

$$u_{\mathrm{p}}\left(x,y\right)=U_{\infty }-{\displaystyle \sum }_{k=1}^N\left[\frac{{\Gamma }_k}{2\pi }\frac{\left(y-y_k\right)}{{\left(x-x_k\right)}^2+{\left(y-y_k\right)}^2}+{\mu }_k\frac{{\left(x-x_k\right)}^2-{\left(y-y_k\right)}^2}{{\left\{{\left(x-x_k\right)}^2+{\left(y-y_k\right)}^2\right\}}^2}\right]$$

(2)

$$v_{\mathrm{p}}\left(x,y\right)={\displaystyle \sum }_{k=1}^N\left[\frac{{\Gamma }_k}{2\pi }\frac{\left(x-x_k\right)}{{\left(x-x_k\right)}^2+{\left(y-y_k\right)}^2}-{\mu }_k\frac{2\left(x-x_{\mathrm{k}}\right)\left(y-y_k\right)}{{\left\{{\left(x-x_k\right)}^2+{\left(y-y_k\right)}^2\right\}}^2}\right]$$

(3)

However, the potential flow is an ideal flow and cannot express the actual velocity deficit or wake generated by each rotor. Therefore, as Whittlesey et al. introduced [23], the component in the x-direction of the potential flow is modified by a wake function du_k showing the velocity deficit of each rotor; the resultant flow u_w(x, y) is expressed by Equation (4). The subscript “w” means the flow field obtained by applying the wake function (or wake model) in this study.

$$u_{\mathrm{w}}\left(x,y\right)=u_{\mathrm{p}}\left(x,y\right)\left\{1-{\displaystyle \sum }_{k=1}^{N}d{u}_k\left(x,y\right)\right\}$$

(4)

In our method, the wake function du_k is given by Equation (5) including the ultra-super-Gaussian function f_{USG_k} defined in our previous study [37]. The function f_{USG_k} modifies the super-Gaussian function 𝑓_{SG_k} proposed by Shapiro et al. [38] expressing the wake profile transformation from top-hat shape to Gaussian by adding the correction function 𝑓_{COR_k} expressing the acceleration regions and the deflection δ_k of the wake. μ_ref is a reference value of the blockage effect and C_{W_k} is a fitting parameter to the prepared flow field of an isolated single rotor.

$$d{u}_k\left(x,y\right)=C_{\mathrm{w}\_k}\frac{{\mu }_k}{{\mu }_{\mathrm{ref}}}{f}_{\mathrm{USG}\_k}=C_{w\_k}\frac{{\mu }_{\mathrm{k}}}{{\mu }_{\mathrm{ref}}}\left\{f_{\mathrm{SG}\_k}-f_{\mathrm{COR}\_k}\right\}$$

(5)

The present method modifies the wake deflection δ_k and the wake width d_W, which are used in the ultra-super-Gaussian function 𝑓_{USG_k}, by using the induced velocity given by the Biot–Savart law. The details are described in Appendix A.

The component in the y-direction of the potential flow is also modified by Equation (6) in our method [40].

$$v_{\mathrm{w}}\left(x,y\right)=v_{\mathrm{p}}\left(x,y\right)+U_{\infty }{\displaystyle \sum }_{k=1}^N\frac{\left|{\Gamma }_k\right|}{\left|{\Gamma }_{\mathrm{SI}}\right|}d{v}_k\left(x,y\right)$$

(6)

where dv_k shows the difference in y-component velocity between the potential flow and the prepared flow (by CFD or experiments) around an isolated single rotor, the circulation of which is defined as Γ_SI and depends on U_∞ (see Equation (A17) in Appendix B). The dv_k is approximated using the sum of four Gaussian-type functions f_{G_i} and four resonant-type functions f_{R_i} as shown in Equation (7). The details of functions f_{G_i} and functions f_{R_i} are described in Appendix A.

$$d{v}_k\left(x,y\right)={\displaystyle \sum }_{i=1}^4{f}_{\mathrm{G}\_i}+{\displaystyle \sum }_{i=1}^4{f}_{\mathrm{R}\_i}$$

(7)

If an arbitrary set of circulations (Γ_k) of N rotors is given, a flow field and a set of rotor outputs are calculated by using Equations (2)–(7) and the prepared relation between the circulation and rotor power. However, the result is not always the actual correct condition. Therefore, we have to find out the set of circulations that give the condition satisfying the conservation of momentum expressed by Newton’s law of motion.

Figure 2 shows a schematic image of a control volume (CV) used for the calculation of the flow field around an isolated single rotor. The prepared data of an isolated single rotor, regardless of CFD data or experimental data, must satisfy Equation (8). The left-hand side of Equation (8) is the total force acting on the fluid in the x-direction (main flow direction) and the right-hand side shows the variation in the momentum in the x-direction per unit time. The variation in the momentum $\mathsf{\Delta }\left(\dot{m}{u}_x\right)$ is calculated from the flow velocity at the boundaries of CV by Equation (9), in which the mass flow rate ${\dot{m}}_x$ and ${\dot{m}}_y$ are expressed by $\rho udy$ and $\rho vdx$, respectively. Here, $\rho$ is the air density; dx and dy are small boundary elements. The total force $F_{\mathrm{total}\_x}$ is given by Equation (10), in which the pressure forces Fp_in and Fp_out acting on the inlet and outlet boundaries are considered and are calculated by the integration along each boundary. The Fp in Equation (10) shows the total pressure force in the x-direction and means the pressure loss. The forces ($F_{\_\tau \_\mathrm{top}} , F_{\_\tau \_\mathrm{bottom}}$) caused by shear stress on the top and bottom boundaries are neglected in this study. From Equations (8) and (10), the relation of Equation (11) is obtained for an isolated single rotor. In our method, the pressure loss Fp has to be prepared as a function of the upstream speed or the corresponding circulation of a single rotor for an applied CV. Note that the pressure loss Fp cannot be calculated from the model flow field obtained by using Equations (4) and (6).

$$F_{\mathrm{total}\_x}=\mathsf{\Delta }\left(\dot{m}{u}_x\right)$$

(8)

$$\mathsf{\Delta }\left(\dot{m}{u}_x\right)=-\int {\dot{m}}_x{u}_{\mathrm{in}}+\int {\dot{m}}_x{u}_{\mathrm{out}}-\int {\dot{m}}_y{u}_{\mathrm{bottom}}+\int {\dot{m}}_y{u}_{\mathrm{top}}$$

(9)

$$F_{\mathrm{total}\_x}=-Th+\int p_{\mathrm{in}}dy-\int p_{\mathrm{out}}dy=-Th+F{p}_{\mathrm{in}}-F{p}_{\mathrm{out}}=-Th+Fp$$

(10)

$$\mathsf{\Delta }\left(\dot{m}{u}_x\right)=-Th+Fp$$

(11)

To apply the momentum balance given by Equation (11) to a wind farm, an evaluation function Err is defined by Equation (12).

$$Err=\mathsf{\Delta }\left(\dot{m}{u}_x\right)+{\displaystyle \sum }_{k}T{h}_k-{\displaystyle \sum }_{k}f{p}_{k}F{p}_k$$

(12)

where Th_k and Fp_k are the thrust force and pressure loss caused by the k-th rotor under the condition of an isolated single rotor, respectively. Th_k and Fp_k are calculated using the input value of circulation (Γ_k) based on the prepared information for an isolated single rotor. fp_k is the correction function introduced to express the interaction between the k-th rotor and other rotors and is defined by Equation (13), in which $I_j$ is a function expressing the interaction effects from the j-th rotor to the k-th rotor.

$$f{p}_k=1-{\displaystyle \sum }_{j\ne k}{I}_j$$

(13)

The layouts between selected two 2D-VAWT rotors can be roughly categorized into four kinds according to the relative rotation condition against the main flow U_∞ as shown in Figure 3; i.e., co-rotating (CO), counter-down (CD), counter-up (CU), and tandem (TD).

In our method, the interaction function $I_j$ is defined separately for each category of the selected two-rotor layout as shown in Equations (14)–(17).

$$I_j={\alpha }_1\frac{\left|{\Gamma }_j\right|}{{\Gamma }_{\mathrm{SI}}}\left|\sin {\phi }_j\right|{\left(\frac{D}{r_{jk}}\right)}^2 ;\mathrm{for} \mathrm{CO}-\mathrm{like}$$

(14)

$$I_j={\alpha }_2\frac{{\Gamma }_j}{{\Gamma }_{\mathrm{SI}}}\sin {\phi }_j{\left(\frac{D}{r_{jk}}\right)}^2 ;\mathrm{for} \mathrm{CD}-\mathrm{like}$$

(15)

$$I_j={\alpha }_3\frac{{\Gamma }_j}{{\Gamma }_{\mathrm{SI}}}\sin {\phi }_j{\left(\frac{D}{r_{jk}}\right)}^2 ;\mathrm{for} \mathrm{CU}-\mathrm{like}$$

(16)

$$I_j={\alpha }_4\frac{\left|{\Gamma }_j\right|}{{\Gamma }_{\mathrm{SI}}}\cos {\phi }_j{\left(\frac{D}{r_{jk}}\right)}^2 ;\mathrm{for} \mathrm{TD}-\mathrm{like}$$

(17)

${\phi }_j$ is defined as the angle between the direction seen from the k-th rotor to the j-th rotor and the upwind direction as shown in Figure 4. $r_{jk}$ is the distance between the centers of the two rotors. A constant angle ${\gamma }_1$ divides the CV into two zones, i.e., the wake zone (in gray) and the out-of-wake zone (in white). If the j-th rotor (counterpart to the k-th rotor) exists in the wake zone, Equation (17) is used for the calculation of the interaction function $I_j$ as the TD-like layout. When the counterpart rotor exists in the out-of-wake zone, the equation used for the calculation of the interaction function $I_j$ is determined as one of three equations (Equations (14)–(16)) according to the corresponding layout category of the two rotors. The fitting parameters of ${\alpha }_1$,${\alpha }_2, {\alpha }_3$, and ${\alpha }_4$ are introduced so as to obtain a sufficient correspondence between the results of rotor power prediction using the model and the CFD analysis, in the specific four conditions of paired rotors, i.e., CO, CD, CU, and TD layouts.

The flow chart of the method is shown in Figure 5. Before starting the calculation, the relations, such as between the circulation and rotor output, and the fitting parameters must be prepared from the information of an isolated single rotor and four specific conditions of two rotors. The first step of the actual model calculation is to decide the calculation conditions such as the upstream wind speed (U_∞), the positions of the rotors, and the rotational directions. The input parameters are the fitting parameters obtained by comparison of the velocity distribution (x- and y-components) of the model with that of the CFD in the case of an isolated single rotor. For example, the fitting parameters include C_N0, C_N1, C_N2, C_N3, C_P0, C_P1, C_P2, and C_P3, which are used to define the correction function f_COR of the ultra-super-Gaussian function and are defined in Equations (A4)–(A9) in Appendix A. The fitting parameters must be determined at several positions in the flow field according to the complexity of the local velocity profiles. The present study uses 35 parameters. All the fitting parameters are shown in

Table A1. The fitting parameters in the ultra-super-Gaussian function, f_USG, (x_n < 0).

xn	Cw	kw	fp	CN0	CN1	CN2	CN3	CP0	CP1	CP2	CP3
−10	0.0091	0.6582	0.7010	0.2542	0.9057	0.2054	0.1845	0.4361	2.1068	0.2156	0.1880
−8	0.0114	0.6905	0.0993	−0.4745	4.7472	2.6759	1.1221	0.2326	0.7857	0.4652	0.1646
−6	0.0168	0.7932	1.3848	−0.5536	5.5377	3.1215	1.3090	0.1482	0.9012	0.5873	0.0326
−4	0.0274	0.9092	0.8591	−0.3956	7.7506	1.4080	0.5688	−0.0644	1.6639	3.3820	0.4386
−2	0.0768	0.7449	0.3266	−0.6307	1.8931	2.0142	0.7956	0.1582	15.5537	5.7160	2.4223
−1	0.1826	0.7898	0.0532	−0.2961	1.0724	1.4556	0.8866	0.1070	13.8469	3.3443	1.4025
−0.75	0.2772	0.5971	0.0423	−0.3472	0.7244	1.4133	0.7344	0.0878	9.8289	2.8971	1.5263
−0.7	0.3000	0.6001	0.0161	−0.4380	0.5724	1.6683	0.7654	0.0849	12.3316	3.0527	0.4652
−0.6	0.2882	0.6330	0.0552	−0.2869	0.6256	1.3519	0.7741	0.1008	9.7876	2.6759	0.5537
−0.5	0.3282	0.5557	0.0002	−0.3087	0.5715	1.2986	0.6040	0.1042	7.4638	2.5409	0.7335
−0.4	0.2887	0.5019	0.0501	−0.3079	0.4466	1.1460	0.6543	0.1325	6.5829	1.9413	0.8771
−0.3	0.4235	0.2310	0.0501	−0.2739	0.3718	0.9942	0.5110	0.1035	5.9861	1.6361	0.7088
−0.2	0.6754	0.2062	0.4530	−0.2698	0.3480	0.1711	0.1711	0.0797	5.2713	1.2941	0.3692
−0.1	0.6232	0.0334	0.1954	−0.0935	0.2087	0.5196	0.5196	0.1074	4.5032	0.2434	0.2433

,

Table A2. The fitting parameters in the ultra-super-Gaussian function, f_USG, (x_n $\ge$ 0).

xn	Cw	kw	fp	CN0	CN1	CN2	CN3	CP0	CP1	CP2	CP3
0	0.5428	0.0002	0.0503	0.0376	3.2977	1.8589	0.4196	0.1575	3.6227	0.2587	0.2586
0.1	0.6577	0.0002	0.3936	0.0636	2.3979	1.7094	0.0005	0.1713	3.0411	0.2162	0.2161
0.2	0.7026	0.0002	0.8119	0.0569	2.7857	0.9608	0.0870	0.1869	2.7598	0.2382	0.2343
0.3	0.5986	0.2500	2.4680	0.0966	2.3364	0.1319	0.0618	0.2673	2.5117	0.4382	0.1426
0.4	0.7250	0.2348	2.5906	0.0879	2.0929	0.7942	0.4227	0.2282	2.4431	0.7059	0.2423
0.5	0.9100	0.2414	2.0028	0.0762	1.8844	0.7625	0.7146	0.2028	2.2853	0.5807	0.2810
0.6	0.9017	0.2699	1.4722	0.1082	1.9000	0.2629	0.2554	0.2149	2.0951	0.6432	0.4851
0.7	0.8632	0.2434	2.0263	0.1524	1.7128	0.4202	0.1218	0.2324	2.1085	0.6836	0.4465
0.75	0.8382	0.2401	2.2507	0.1566	1.7334	0.5260	0.1615	0.2393	2.1085	0.6402	0.4687
0.8	0.8504	0.2306	2.3917	0.1651	1.7208	0.5245	0.1327	0.2456	2.0862	0.6898	0.4399
0.9	0.8672	0.2166	2.4366	0.1739	1.7442	0.5439	0.0802	0.2452	2.1060	0.6962	0.4270
1	0.8762	0.2012	2.2059	0.2108	1.5943	0.8172	0.0553	0.2418	2.1336	0.7414	0.4284
1.1	0.8914	0.1931	2.1107	0.2157	1.6410	0.8202	0.0080	0.2379	2.1930	0.7759	0.3961
1.2	0.8879	0.1777	1.9914	0.2077	1.6346	0.8721	0.1508	0.2282	2.2008	0.8100	0.5021
1.3	0.9051	0.1636	1.8485	0.2067	1.7130	0.9247	0.0650	0.2064	2.2741	0.7501	0.6179
1.4	0.9128	0.1570	1.7436	0.2088	1.7008	0.9034	0.0864	0.2033	2.2947	0.7559	0.6180
1.5	0.9219	0.1498	1.6507	0.2090	1.7150	0.8972	0.0837	0.2025	2.3173	0.8285	0.5823
1.6	0.9303	0.1494	1.3965	0.2043	1.7353	0.8768	0.0966	0.1889	2.4271	0.7589	0.6133
1.7	0.9296	0.1336	1.4457	0.1977	1.7899	0.9275	0.1035	0.1875	2.4322	0.8566	0.6182
1.8	0.9378	0.1281	1.2945	0.1940	1.8314	0.9296	0.0805	0.1801	2.4907	0.8648	0.6216
1.9	0.9436	0.1284	1.1287	0.1920	1.8235	0.8183	0.0833	0.1748	2.5667	0.9122	0.5759
2	0.9740	0.1333	1.3191	0.1849	1.8677	0.3759	0.0424	0.1971	2.6559	0.9214	0.0583
2.1	0.9750	0.1277	1.1529	0.1757	1.9256	0.3570	0.0216	0.1867	2.7499	0.9503	0.0442
2.2	0.9441	0.1243	1.1448	0.1661	2.0006	0.4031	0.0826	0.1743	2.8841	0.9444	0.1459
2.3	0.9400	0.1194	1.0655	0.1572	2.0721	0.3583	0.0821	0.1666	2.9870	0.9797	0.1179
2.4	0.9266	0.1147	0.9121	0.1495	2.1548	0.3858	0.0739	0.1531	3.0981	0.9936	0.2387
2.5	0.9105	0.1129	0.8607	0.1441	2.2145	0.3849	0.0788	0.1480	3.2481	1.0166	0.1923
3	0.7977	0.1058	0.6963	0.1077	2.8396	0.3861	0.0758	0.1184	4.2106	1.1886	0.1689
4	0.5657	0.1063	0.4007	0.3131	2.1568	6.6251	0.1685	0.3866	2.6141	6.8336	0.3954
6	0.3571	0.1143	0.3388	0.1253	3.5691	7.9674	1.5521	0.2249	3.9554	7.7057	1.1499
8	0.3064	0.1051	0.0357	0.0941	4.8910	8.3139	0.8055	0.1531	5.3229	7.8979	1.2282
10	0.2920	0.0905	0.3596	0.1307	3.5955	8.6259	1.2900	0.1461	5.2149	7.7429	1.0662

,

Table A3. The fitting parameters in Gaussian-type functions, f_{G_i}, (x_n < 0).

xn	C1fg1	C2fg1	C2fg1	C1fg2	C2fg2	C2fg2	C1fg3	C2fg3	C2fg3	C1fg4	C2fg4	C3fg4
−10	0.0235	0.1686	−0.3653	−0.0231	0.1752	−0.3619	−0.0597	0.0938	1.5479	0.0597	0.0952	1.5486
−8	0.0108	0.1888	−0.2683	−0.0081	0.1372	−0.3694	−0.0103	0.1976	0.1644	0.0078	0.1417	0.2737
−6	0.0135	0.1554	−0.3257	−0.0159	0.1723	−0.2669	−0.0427	0.1721	0.1944	0.0448	0.1779	0.1726
−4	0.0229	0.4217	−1.3449	−0.0224	0.4184	−1.3278	−0.0531	0.0965	0.4337	0.0540	0.0973	0.4337
−2	0.0082	0.3241	−1.1351	−0.0044	0.3085	−1.0243	0.0142	0.2679	0.3599	−0.0101	0.7884	0.4118
−1	0.0250	0.4050	−1.2649	−0.0250	0.3950	−1.1774	0.0156	0.4366	0.2983	−0.0092	1.1056	0.5103
−0.75	−0.0212	0.2538	−0.5622	0.0362	0.1196	−0.2382	0.0133	0.1835	0.9382	−0.0112	0.4931	0.8562
−0.7	−0.0198	0.1009	−0.9601	−0.0194	0.0485	−0.5937	−0.0272	0.2134	0.5007	0.0221	0.7050	0.1507
−0.6	−0.0421	0.5509	−0.2054	0.0235	0.0507	−0.5007	0.0202	0.2292	0.7502	−0.0123	1.0304	0.4118
−0.5	−0.1213	1.1116	−0.4966	0.0780	1.1743	−0.7618	0.1627	0.1053	0.8753	−0.1438	0.1041	0.8753
−0.4	−0.1678	0.9429	−0.7003	0.1329	1.0688	−0.8001	0.0942	0.1765	0.3735	−0.1902	0.1649	0.2601
−0.3	−0.0451	0.6440	−0.7994	0.0202	1.1743	−1.0633	−0.0198	0.1237	0.5260	−0.0968	0.0549	0.4159
−0.2	−0.0870	0.3531	−0.4945	0.0488	0.1115	−0.6442	−0.0277	0.0174	0.4979	−0.1414	0.1458	0.2983
−0.1	−0.0111	0.1642	−1.1036	0.0371	0.0441	−0.5663	0.0461	0.2134	0.8890	−0.1499	0.5382	0.2635

,

Table A4. The fitting parameters in Gaussian-type functions, f_{G_i}, (x_n $\ge$ 0).

xn	C1fg1	C2fg1	C2fg1	C1fg2	C2fg2	C2fg2	C1fg3	C2fg3	C2fg3	C1fg4	C2fg4	C3fg4
0	−0.0282	0.2644	−0.7399	0.0321	0.2067	−0.5096	−0.0188	0.1606	0.6415	−0.1988	0.1055	0.3195
0.1	−0.0101	0.1185	−1.0749	0.0648	0.0471	−0.5048	−0.0347	0.0912	0.5759	−0.1326	0.1690	0.1507
0.2	−0.0086	0.1905	−0.8931	0.0829	0.0559	−0.4788	0.1798	0.0244	0.5212	−0.1983	0.0467	0.4911
0.3	−0.0207	0.4419	−0.3236	0.0560	0.1071	−0.4590	0.0998	0.1044	0.3243	−0.0432	0.0118	0.7140
0.4	−0.0301	0.2881	−0.8644	0.0458	0.5590	−0.4720	0.0930	0.0640	0.4118	−0.0321	0.0200	0.7194
0.5	−0.0782	0.4384	−0.6955	0.1074	0.5869	−0.4583	0.0617	0.1149	0.4528	−0.0458	0.0221	0.7659
0.6	−0.0075	0.1009	−0.8507	0.0620	0.1972	−0.3284	0.0980	0.0508	0.4549	−0.0153	0.1649	0.3134
0.7	−0.0116	0.2116	−0.5554	0.0341	0.0302	−0.4679	0.0823	0.6071	0.2744	−0.0290	0.8212	0.6627
0.75	−0.0559	0.4797	−0.5663	0.0822	0.6792	−0.2601	0.0606	0.8796	0.3729	−0.0480	1.1384	0.5431
0.8	−0.0641	0.4331	−0.5472	0.0928	0.5810	−0.2874	0.0702	0.8866	0.4494	−0.0569	1.0591	0.6319
0.9	−0.0324	0.2019	−0.6852	0.0464	0.5345	−0.4195	0.2196	0.8986	0.2774	−0.1847	1.1113	0.2855
1	−0.0195	0.0280	−0.6645	0.0231	0.0927	−0.3909	0.1242	0.8241	0.7064	−0.1016	0.6778	0.9532
1.1	−0.0107	0.2134	−0.5882	0.0321	0.0485	−0.3264	0.1567	0.9455	0.5062	−0.1289	1.0222	0.6326
1.2	−0.0219	0.1079	−0.4754	0.0590	0.0661	−0.3271	0.1216	0.9341	0.6162	−0.1017	0.9743	0.7481
1.3	−0.0558	0.0745	−0.5868	0.0794	0.1665	−0.3756	0.0442	0.1009	0.5868	−0.0063	1.0618	1.4181
1.4	−0.0680	0.0872	−0.5884	0.0808	0.1703	−0.4052	0.0379	0.0807	0.5984	−0.0082	0.8341	1.3434
1.5	−0.0694	0.0877	−0.5999	0.0810	0.1966	−0.4046	0.0406	0.0807	0.6152	−0.0074	1.1495	1.6060
1.6	−0.0722	0.0915	−0.5991	0.1011	0.2151	−0.3756	0.0540	0.1458	0.6705	−0.0273	0.7772	0.5111
1.7	−0.0551	0.0708	−0.5765	0.0867	0.1519	−0.3620	0.0340	0.0817	0.6651	−0.0060	1.3904	1.0672
1.8	−0.0155	0.0552	−0.5759	0.0634	0.0961	−0.2594	0.1763	0.0736	0.6538	−0.1443	0.0795	0.6524
1.9	−0.0372	0.0745	−0.5528	0.0779	0.1463	−0.2945	0.0737	0.2138	0.7138	−0.0492	0.4508	0.6401
2	−0.0160	0.0745	−0.5143	0.0693	0.0894	−0.2747	0.1489	0.0990	0.7382	−0.1203	0.1036	0.7424
2.1	−0.0602	0.0877	−0.5690	0.0880	0.1811	−0.3981	0.0455	0.3259	0.6873	−0.0301	0.8226	0.5003
2.2	−0.0496	0.0991	−0.5868	0.0715	0.2031	−0.3920	0.0619	0.1132	0.6155	−0.0466	0.1882	0.5007
2.3	−0.0718	0.1378	−0.5964	0.0937	0.2646	−0.4255	0.0416	0.1352	0.6401	−0.0276	0.5382	0.2819
2.4	−0.0469	0.1202	−0.5868	0.0655	0.2646	−0.3790	0.0355	0.1202	0.6436	−0.0234	0.4507	0.3257
2.5	−0.0363	0.1149	−0.6060	0.0529	0.2998	−0.3735	0.0284	0.1158	0.6360	−0.0172	0.5382	0.2382
3	−0.0111	1.4491	−1.3743	0.0219	0.8227	−0.7618	0.0533	0.3804	0.3011	−0.0405	0.7993	0.2054
4	0.0096	0.1062	−1.1569	−0.0037	0.2368	−1.0284	0.0125	0.1835	0.3284	−0.0099	0.2470	1.3101
6	0.0224	0.6458	−0.8438	−0.0217	0.7993	−0.8329	−0.0067	0.0921	0.2259	0.0139	0.1854	0.1944
8	0.0231	0.4278	−0.6121	−0.0225	0.3979	−0.6374	−0.0258	0.4929	0.5649	0.0301	0.4274	0.5417
10	0.0247	0.4050	−1.0838	−0.0246	0.3745	−1.0626	−0.0270	0.4929	0.4884	0.0312	0.4425	0.4774

,

Table A5. The fitting parameters in resonance-type functions, f_{R_i}, (x_n < 0).

xn	C1fr1	C2fr1	C3fr1	C1fr2	C2fr2	C3fr2	C1fr3	C2fr3	C3fr3	C1fr3	C2fr3	C3fr3
−10	0.2350	0.0101	−0.0006	−0.2396	0.0101	0.3111	1.6824	1.4410	−0.0001	−1.6824	1.4410	0.0001
−8	0.2405	0.0087	−0.2206	−0.2456	0.0087	0.2860	1.0335	1.3209	−0.0001	−1.0335	1.3209	0.0001
−6	0.2470	0.0161	−0.0988	−0.2538	0.0163	0.4355	1.2991	0.6408	−0.0004	−1.2991	0.6408	0.0004
−4	0.2460	0.0237	−0.4757	−0.2534	0.0237	0.2733	1.7798	0.8008	0.0007	−1.7798	0.8008	−0.0007
−2	0.3042	0.0648	−0.4639	−0.3168	0.0652	0.3341	2.6455	1.2809	0.0024	−2.6455	1.2809	−0.0024
−1	0.2864	0.0647	−0.6180	−0.3014	0.0660	0.1831	4.6927	0.6408	0.0159	−4.6927	0.6408	−0.0159
−0.75	0.2980	0.0607	−0.6671	−0.3132	0.0622	0.0819	5.0415	0.8016	0.0214	−5.0415	0.8016	−0.0214
−0.7	0.2742	0.0568	−0.2993	−0.2964	0.0632	0.4784	5.5835	0.8906	0.0209	−5.5835	0.8906	−0.0209
−0.6	0.2550	0.0575	−0.4743	−0.2768	0.0636	0.3882	6.0181	1.0445	0.0213	−6.0181	1.0445	−0.0213
−0.5	0.2636	0.0556	−0.5181	−0.2821	0.0600	0.3116	6.0327	1.1109	0.0240	−6.0327	1.1109	−0.0240
−0.4	0.2591	0.0565	−0.5181	−0.2791	0.0619	0.3062	6.0938	1.1992	0.0259	−6.0938	1.1992	−0.0259
−0.3	0.2593	0.0562	−0.4743	−0.2776	0.0613	0.3062	6.0913	1.0992	0.0268	−6.0913	1.0992	−0.0268
−0.2	0.2771	0.0544	−0.2200	−0.2958	0.0603	0.4757	6.1938	1.0992	0.0275	−6.1938	1.0992	−0.0275
−0.1	0.2784	0.0538	−0.1899	−0.2959	0.0594	0.4757	6.1182	0.9992	0.0273	−6.1182	0.9992	−0.0273

and

Table A6. The fitting parameters in resonance-type functions, f_{R_i}, (x_n $\ge$ 0).

xn	C1fr1	C2fr1	C3fr1	C1fr2	C2fr2	C3fr2	C1fr3	C2fr3	C3fr3	C1fr3	C2fr3	C3fr3
0	0.2802	0.0552	−0.1681	−0.2976	0.0613	0.4873	6.1572	0.9992	0.0267	−6.1572	0.9992	−0.0267
0.1	0.2839	0.0538	−0.3431	−0.2976	0.0577	0.3075	6.0498	0.9992	0.0267	−6.0498	0.9992	−0.0267
0.2	0.2344	0.0547	−0.4962	−0.2447	0.0578	0.3116	5.8423	1.0453	0.0269	−5.8423	1.0453	−0.0269
0.3	0.2285	0.0538	−0.5071	−0.2355	0.0557	0.3021	5.7690	0.9570	0.0257	−5.7690	0.9570	−0.0257
0.4	0.2608	0.0557	−0.2966	−0.2665	0.0575	0.3896	5.6689	0.9008	0.0244	−5.6689	0.9008	−0.0244
0.5	0.2589	0.0510	−0.2316	−0.2619	0.0519	0.4347	5.3271	0.8008	0.0244	−5.3271	0.8008	−0.0244
0.6	0.2679	0.0500	−0.1260	−0.2689	0.0505	0.5311	5.3979	0.8008	0.0214	−5.3979	0.8008	−0.0214
0.7	0.2399	0.0501	−0.0382	−0.2391	0.0505	0.7026	4.7998	0.8008	0.0217	−4.7998	0.8008	−0.0217
0.75	0.2698	0.0542	0.0493	−0.2668	0.0539	0.7136	4.1235	0.8039	0.0238	−4.1235	0.8039	−0.0238
0.8	0.2462	0.0580	0.0493	−0.2421	0.0575	0.7806	4.5215	0.8227	0.0196	−4.5215	0.8227	−0.0196
0.9	0.2860	0.0580	0.0518	−0.2794	0.0566	0.6763	4.8775	0.7918	0.0158	−4.8775	0.7918	−0.0158
1	0.2571	0.0519	0.0542	−0.2494	0.0501	0.7363	5.2356	0.7008	0.0131	−5.2356	0.7008	−0.0131
1.1	0.2719	0.0569	0.0542	−0.2620	0.0544	0.7061	4.3723	0.7208	0.0117	−4.3723	0.7208	−0.0117
1.2	0.1681	0.0445	0.0710	−0.1605	0.0426	1.0612	2.2775	0.5074	0.0210	−2.2775	0.5074	−0.0210
1.3	0.3056	0.0543	−0.2436	−0.2901	0.0497	0.3166	2.4834	0.6488	0.0115	−2.4834	0.6488	−0.0115
1.4	0.3011	0.0707	0.0688	−0.2866	0.0673	0.6694	3.0307	1.0250	0.0035	−3.0307	1.0250	−0.0035
1.5	0.3249	0.0551	−0.2494	−0.3071	0.0503	0.2907	2.3580	1.6010	0.0019	−2.3580	1.6010	−0.0019
1.6	0.2896	0.0539	−0.2406	−0.2722	0.0488	0.3316	2.5045	1.6788	−0.0028	−2.5045	1.6788	0.0028
1.7	0.2894	0.0482	−0.1825	−0.2739	0.0442	0.3709	2.3087	1.5857	−0.0061	−2.3087	1.5857	0.0061
1.8	0.2839	0.0425	−0.3007	−0.2681	0.0386	0.2556	1.9482	1.1992	−0.0108	−1.9482	1.1992	0.0108
1.9	0.3027	0.0386	−0.0691	−0.2893	0.0361	0.4328	2.6024	1.3038	−0.0101	−2.6024	1.3038	0.0101
2	0.2838	0.0284	−0.1701	−0.2721	0.0264	0.3567	2.1167	1.6008	−0.0142	−2.1167	1.6008	0.0142
2.1	0.2867	0.0319	−0.0874	−0.2739	0.0297	0.4251	4.2651	1.0992	−0.0097	−4.2651	1.0992	0.0097
2.2	0.2430	0.0313	−0.1366	−0.2300	0.0285	0.4306	4.7534	1.2008	−0.0100	−4.7534	1.2008	0.0100
2.3	0.2664	0.0291	−0.0942	−0.2542	0.0270	0.4251	5.3882	1.0398	−0.0098	−5.3882	1.0398	0.0098
2.4	0.2539	0.0294	−0.0901	−0.2412	0.0271	0.4306	5.8081	1.0281	−0.0100	−5.8081	1.0281	0.0100
2.5	0.2307	0.0250	−0.1487	−0.2185	0.0228	0.4444	5.9830	0.9945	−0.0104	−5.9830	0.9945	0.0104
3	0.2566	0.0226	−0.1334	−0.2425	0.0203	0.3124	4.3969	0.8562	−0.0176	−4.3969	0.8562	0.0176
4	0.2470	0.0164	0.0075	−0.2347	0.0147	0.3545	5.3616	0.5986	−0.0135	−5.3616	0.5986	0.0135
6	0.2499	0.0098	0.0062	−0.2410	0.0089	0.1392	2.5100	0.4103	−0.0108	−2.5100	0.4103	0.0108
8	0.2223	0.0098	0.0395	−0.2160	0.0094	0.1740	2.9275	0.4103	−0.0014	−2.9275	0.4103	0.0014
10	0.2540	0.0253	0.0462	−0.2489	0.0254	0.1833	3.1876	0.3155	0.0001	−3.1876	0.3155	−0.0001

of Appendix A. After the initial condition of the rotors and a set of parameters are input, the process of searching for a reasonable combination of circulations is initiated starting from the lowest value of the circulation in a given searching range for each rotor (determining the searching ranges will be described in Section 2.2). For a given combination of circulations, the potential flow is calculated first and it is modified by our wake model to include the effects of the velocity deficits. With the flow field (u_w, v_w) obtained, the momentum changes in the x-direction per unit of time $\mathsf{\Delta }\left(\dot{m}{u}_x\right)$ is calculated from the four boundary conditions. The thrust force Th_k and the (isolated single rotor) pressure loss Fp_k are calculated by Equations (A19) and (A20) in Appendix B by using the given circulation value for each rotor. Using the interaction function $I_j$, the correction function fp_k which expresses the effects of the inter-rotor interaction is obtained for each rotor by Equation (13). Then, the momentum balance Err is evaluated by Equation (12). The “error” of momentum balance is defined in Equation (18) in the present study. The denominator in Equation (18) is the momentum per unit time entering into the CV from the inlet boundary.

$$error=\frac{\sqrt{Er{r}^2}}{\int {\dot{m}}_x{u}_{\mathrm{in}}}$$

(18)

The above calculation is repeated by changing the combination of circulations with the interval of ΔΓ throughout all the searching ranges. After the round robin, a combination of circulation Γ_k which gives the smallest “error” is obtained. In the present study, the similar searching process (sub-searching) is repeated twice by narrowing the searching ranges as (Γ_k − ΔΓ) $\le$ Γ_{k_i} $\le$ (Γ_k + ΔΓ). Finally, for the smallest “error” condition, the power output of each rotor can be obtained by using Equation (A21) in Appendix B.

2.2. Specific Calculation and Searching Range of Circulation

In this section, we consider a specific case assuming small rotors in order to explain the selection of searching range to obtain adequate circulation. The target rotors are the same as those investigated in the previous study [35] and the rotor configuration is depicted in Figure A4 in Appendix B. The necessary relations on the performance of the single rotor are given in Equation (A17)–(A20). The values of ${\alpha }_1$,${\alpha }_2, {\alpha }_3$, and ${\alpha }_4$ were determined as 0.2304, 0.19, −0.2853, and 1.07, respectively, from the preliminary calculations of the specific four layouts of paired rotors. The CV was defined as 20D × 20D × 0.868D. The thickness of the CV given by 0.868D corresponds to the rotor height of the target rotor; by considering the rotor height, our 2D method can be applied to the prediction of the power output of an actual 3D-VAWT rotor.

Figure 6a,d show the averaged distributions of the x- and y-component velocities calculated by the CFD analysis around two VAWT rotors rotating CCW direction in the CO-like layout. The commercial code STAR-CCM+ was utilized for the CFD simulation. The circulation of the upper rotor (rotor-k) is 0.334 m²/s, which was obtained by the integration of the flow field along a circular path at r/D = 0.75. The circulation of the lower rotor (rotor-j) is 0.345 m²/s. The details of the rotor shape are not shown but are replaced by circles in white in Figure 6a. The black solid lines around the rotors are path lines.

Figure 6b,e show the resultant flow field obtained by the proposed method using an in-house code, in which the circulation of each rotor was changed step by step in a round robin over a wide searching range from 0.1${\Gamma }_{\mathrm{SI}}$ to 1.1${\Gamma }_{\mathrm{SI}}$. Here, ${\Gamma }_{\mathrm{SI}}$ is the circulation of an isolated single rotor obtained from the CFD analysis and the value is 0.326 m²/s in the case of U_∞ = 10 m/s. The obtained circulations of two rotors in Figure 6b are Γ_k = 0.072 m²/s and Γ_j = 0.037 m²/s, respectively, and they are different from the CFD results in Figure 6a. To demonstrate the variation in the error of momentum balance graphically, we define a new indicator (1 − error) that approaches 1 when the error becomes small. The distribution of the values of (1 − error) obtained in the searching process of the smallest error in the case of Figure 6b is shown in Figure 7, which includes 625 results corresponding to the combinations of circulation values (Γ_j, Γ_k) with the interval of (1.1 − 0.1)${\Gamma }_{\mathrm{SI}}$/25. Figure 7 shows that there are a lot of combinations of circulation which might give a small error.

In the parallel layouts, it is expected that two rotors have large circulations close to the value of ${\Gamma }_{\mathrm{SI}}$. Therefore, the search range can be narrowed. Figure 6c,f are is the flow field obtained using a searching range limited from 0.95${\Gamma }_{\mathrm{SI}}$ to 1.1${\Gamma }_{\mathrm{SI}}$ for each rotor in the CO layout. The circulation combination giving the smallest error, shown in Figure 6c, is Γ_k = 0.337 m²/s and Γ_j = 0.349 m²/s, respectively. Figure 8 is the distribution of the 625 values of (1 − error) obtained in the searching process with a limited range from 0.95${\Gamma }_{\mathrm{SI}}$ to 1.1${\Gamma }_{\mathrm{SI}}$ for each rotor. The results of our method with appropriately limited searching ranges can approach the CFD results.

Incorrect results are also obtained when the rotor at issue exists in the wake of other rotors if the wide searching range of circulation is applied. In our method, we propose a procedure to search for the appropriate result close to the CFD result, in which the searching range of circulation of each rotor in a VAWT wind farm is limited in the manner described below.

Let us define a circular region of the radius R_int around a rotor (k) at issue, and assume there are N_int rotors in the region. We consider that only the rotors (j) included in the interaction region defined by R_int are used to decide the searching range of circulation of the rotor (k).

First of all, the angle ${\phi }_j$ defined in Figure 4 is checked for all rotors in the interaction region. If the absolute value of one of the angles is less than a constant angle ${\gamma }_2$, the searching range expressed by Equation (19) is applied. Note that the angle ${\gamma }_2$ is different from the angle ${\gamma }_1$, and ${\gamma }_2$ is less than ${\gamma }_1$. The interaction region is determined to be a circle of R_int = 10D in this specific calculation; therefore, the maximum of the inter-rotor distance (between rotor centers) r_jk is 10D in Equation (19) in this study.

$$0.3{\Gamma }_{\mathrm{SI}}\le {\Gamma }_k<\left(0.075\frac{r_{jk}}{D}+0.35\right){\Gamma }_{\mathrm{SI}} ; \mathrm{for} \left|{\phi }_j\right|<{\gamma }_2$$

(19)

In the second step of the procedure, the average of the distance d_j between the center of the j-th rotor and the stream-wise center line through the k-th rotor at issue is calculated with all the counterpart rotors of the k-th rotor in the interaction region using Equation (20).

$$\overline{d_j}=\frac{{{\displaystyle \sum }}_{j\ne k}\left|d_j\right|}{N_{\mathrm{int}}-1}$$

(20)

In the third step of the procedure, the searching range of the circulation Γ_k of the k-th rotor is determined as one of Equations (21)–(23) according to the averaged lateral distance $\overline{d_j}$. Figure 9 schematically shows the definition of the angle ${\gamma }_2$ and three zones determining the searching range of circulation of the k-th rotor. The delimiting lateral distances are determined as d₁ = 0.9D and d₂ = 1.5D in this study.

$$0.85{\Gamma }_{\mathrm{SI}}\le {\Gamma }_k<1.0{\Gamma }_{\mathrm{SI}} ; \mathrm{for} 0\le \overline{d_j}<d_1$$

(21)

$$0.9{\Gamma }_{\mathrm{SI}}\le {\Gamma }_k<1.0{\Gamma }_{\mathrm{SI}} ; \mathrm{for} d_1\le \overline{d_j}<d_2$$

(22)

$$0.95{\Gamma }_{\mathrm{SI}}\le {\Gamma }_k<1.1{\Gamma }_{\mathrm{SI}} ; \mathrm{for} \overline{d_j}\ge d_2$$

(23)

3. Results and Discussion

The method described in Section 2 is validated in this section by applying it to several layouts consisting of the same small 2D rotors as used in the previous section (Section 2.2). Firstly, in Section 3.1, the power dependency of the isolated single rotor on the wind speed is confirmed. Then, in Section 3.2, the power dependency of the paired rotors on the 16 wind directions is compared between the present method and the CFD analysis obtained in the previous study. As for the more complicated layouts, the power dependency of the three-rotor cluster on the 12 wind directions is investigated in Section 3.3. Finally, to show the applicability of the proposed method to VAWT wind farms, the power prediction in the four-rotor layouts arranged in a line is tried in Section 3.4.

3.1. An Isolated Single Rotor

Table 1. Predicted power output of an isolated single rotor by the present method and the comparison with the CFD results (reproduced with permission from Hara et al. [35]).

U∞ (m/s)	P_Model (mW)	P_CFD (mW)	% Error
4	8.14	7.99	1.88
6	32.87	31.88	3.11
8	85.50	83.76	2.08
10	177.62	177.62	0.0
12	314.23	322.17	−2.46

shows the comparison between the CFD analysis [35] and the model simulation using the present method based on the momentum balance in the cases of an isolated single rotor in the five different upstream wind speeds (U_∞). Although the original CFD analysis was conducted using a wide calculation domain of 40 D × 50 D, the square region of 20 D × 20 D enclosing the single rotor is extracted as the CV, which is the same size as that of the model simulation, in order to acquire the dependence of the pressure loss Fp on the circulation of the isolated single rotor. The output power shown in the unit of “mW” is equivalent to that of the small rotor, the size of which is 50 mm in diameter × 43.4 mm in height. The parameters necessary to rebuild the flow field calculated by the CFD were obtained at the reference wind speed of U_∞ = 10 m/s. Therefore, the smallest percentage difference (% Error) between the model and CFD results is obtained in the case of U_∞ = 10 m/s. The output power predicted for a single rotor in other wind speeds agrees well with the CFD result; this means the method using momentum conservation works well.

3.2. Paired Rotors

In this section, the prediction of the powers of paired rotors is investigated. Figure 10 shows the definition of relative wind direction ($\theta$) to the rotor pair which is categorized into two configurations in terms of the relative rotational direction of two rotors. The configuration in Figure 10a is defined as the co-rotation (CO), which includes the CO and TD layouts where two rotors rotate in the same direction (see Figure 3). On the other hand, the configuration in Figure 10b is defined as the inverse rotation (IR), which includes the CD, CU, and TD layouts where two rotors rotate mutually in opposite directions. The dotted lines in red in Figure 10 correspond to the boundaries of the wake zone defined by the angle ${\gamma }_1$ in Figure 4. The angle ${\gamma }_1$ is set at 30° in this study.

The power outputs of the upper (Rotor 1: R1) and lower (Rotor 2: R2) rotors and the averaged power of both rotors were predicted in each of the 16 wind directions by the present method under the condition of the upstream wind speed U_∞ of 10 m/s. The results in the case the inter-rotor space (gap) is equal to the rotor diameter are shown in Figure 11a–c for CO and Figure 12a–c for the IR configurations. Each of the predicted rotor power is normalized by the power of the isolated single rotor (SI) and compared with the normalized CFD results in Figure 11 and Figure 12. The conditions of $\theta$ = 0° and 180° in Figure 11c are the same CO layout and their output powers were adjusted to agree with the CFD results to get the adequate value of the parameter ${\alpha }_1$. The conditions of $\theta$ = 0° and 180° in Figure 12c correspond to the CD and CU layouts, respectively, and the values of the parameters ${\alpha }_2$ and ${\alpha }_3$ were adjusted so as to agree with the CFD results. Similarly, the conditions of $\theta$ = 90° and 270° in Figure 11c and Figure 12c correspond to either of the two TD layouts shown in Figure 3d, and the parameter ${\alpha }_4$ can be determined by fitting the resultant circulations obtained by the model simulation to those of the CFD analysis in one of the four TD conditions. The actual values of ${\alpha }_1$,${\alpha }_2, {\alpha }_3$, and ${\alpha }_4$ were already shown in Section 2.2. Although the difference between the model simulation and the CFD analysis is somewhat large in a few wind directions as shown in Figure 11 and Figure 12, the method can predict the power of a rotor pair with the gap = 1.0D very well.

Figure 13 and Figure 14 are the results of the rotor power distribution of the paired rotors with the gap = 0.5D. The model simulation seems to underestimate the rotor powers compared with the CFD analysis in many wind directions.

Figure 15 and Figure 16 show the results in the case of gap = 2.0D. Except for a few directions, the model and CFD predictions agree well in the long inter-rotor distance case. This fact may suggest the necessity to modify the interaction function $I_j$ defined in Equations (14)–(17) in the case of short inter-rotor distance.

Although further improvements may be needed in the present method, it is worth noting that the prediction of the power of paired rotors in a specific wind direction can be performed by the method in about 40 min, which is 500 times shorter than the calculation time with the CFD analysis using the DFBI model (about 2 weeks).

3.3. Three-Rotor Cluster

The three-rotor cluster is also categorized into the CO and IR configurations. Figure 17 shows the definition of the 12 wind directions for the two configurations in the cases of three-rotor clusters which are arranged like a triangular shape with an inter-rotor space of gap = 1.0D. Under the condition of the upstream wind speed U_∞ of 10 m/s, the prediction of the averaged output power of the rotor clusters was carried out. All the parameters and necessary relations such as the pressure loss function are the same as those used in the case of paired rotors except for the division of the searching range of circulation. The number of divisions of the searching range was set as 10 in the three-rotor case instead of the 25 used in paired rotors case to save the calculation time. Therefore, the error in the prediction is anticipated to be a little worse.

Figure 18 is the prediction results of the averaged power of the three rotors. The model calculation underestimates the averaged power of the rotor clusters. However, the dependence on the wind direction is well simulated and the model simulation gives the same trend as the CFD analysis [41].

3.4. Four-Rotor Layout

Finally, to show the applicability of the proposed method to VAWT wind farms, four-rotor layouts arranged in line with an inter-rotor space of 3.0D were selected as the targets of the output power prediction. To cover a large number of rotors, the CV was set to be a large size of 40D × 40D × 0.868D. As the size of the CV was changed, the function expressing the pressure loss Fp of the CV was newly prepared from the CFD flow field of the isolated single rotor; i.e., the reference value Fp_ref in Equation (A20) in Appendix B was changed from 229.61 to 347.95 mN. The upstream wind speed was assumed to be 10 m/s. The number of divisions of the searching range of the circulation of each rotor was set as 8. In this case, the calculation time is about 50 times as long as the calculation time in a two-rotor case with the 25 division because the subdivision process is executed twice in the actual calculation using the in-house code.

Figure 19a–d show the results of the power prediction and the distributions of x- and y-direction velocity components in two kinds of four-rotor layouts, i.e., parallel and tandem, respectively. The fitting parameters ${\alpha }_1$,${\alpha }_2, {\alpha }_3$, and ${\alpha }_4$ are set to be the same values as those used in the paired rotors (see Section 2.2). The model predicts the decrease in the rotor power in the order of R1, R2, R3, and R4 in the parallel array of four rotors. Similarly, in the tandem array of four rotors, the output power decreases in the same order.

In this study, the CFD analysis of the two kinds of four-rotor layouts was carried out by using the DFBI model. The mesh size and whole domain size of each simulation of the four-rotor layouts are the same as those used in the previous CFD analysis of one-, two-, and three-rotor arrangements. The details of the meshes created for the four-rotor layouts in the present study are shown in Appendix C. The upstream wind speed U_∞ is set at 10 m/s in the CFD analysis of four-rotor layouts.

Figure 20a,b show the distributions of x-component of unsteady velocity obtained by the CFD for the parallel and tandem layouts, respectively. Figure 20c,d illustrate the unsteady flow field shown by the y-component velocity for the two layouts of four rotors. The condition shown in Figure 20 corresponds to the state at 4 s from the beginning of the simulation. The time history of the angular velocity of each rotor is shown in Appendix C. The CFD analysis was conducted by a high-performance PC with 28 cores; it took about 1 week for each calculation of two layouts.

Table 2. Comparison of each rotor power in the parallel and tandem layouts between the present model and the CFD analysis.

Layout	Parallel				Tandem
Rotor	R1	R2	R3	R4	R1	R2	R3	R4
Power (model) (mW)	181.0	191.3	205.6	209.3	170.0	135.6	47.7	42.5
Power (CFD) (mW)	185.3	211.8	216.0	236.2	147.8	47.4	22.8	3.0
Err.P (%)	−2.4	−11.5	−5.8	−15.1	12.5	49.5	14.0	22.2

shows the comparison of each rotor power in the parallel and tandem layouts between the present model and CFD analysis. The percentage error of the rotor power (Err.P) is defined by Equation (24) in this study. P_Model is the rotor power predicted by our model using Equation (A21) and the P_CFD is the CFD result of the rotor power which is obtained by multiplying the averaged angular velocity and the averaged rotor torque during the last 0.5 s. The P_ref is the output power of an isolated single rotor.

$$Err.P=\frac{P_{\mathrm{Model}}-P_{\mathrm{CFD}}}{P_{\mathrm{ref}}}\times 100$$

(24)

The comparison shown in Table 2 does not give satisfactory results. In particular, the difference between the model and CFD is large in the tandem layout of four rotors. One of the reasons for the disagreement is the small number of initial divisions of the searching range used for the model prediction. An advanced method like the genetic algorithm (GA) should be used to give the possible combinations of circulations at random. Further improvement in the model of the interaction among rotors also is necessary. Nevertheless, the present model can predict the same trend in the rotor power in the parallel and tandem layouts of four rotors as that of the CFD analysis.

4. Conclusions

We developed a method to predict the output powers of the vertical-axis wind turbine (VAWT) rotors of an arbitrary layout; the calculation time can be incomparably shorter than computational fluid dynamics (CFD) using the dynamic fluid/body dynamics (DFBI) model, which enables the time variation in rotor speed. However, our method needs to know information, such as the flow velocity distributions and pressure loss and circulation, and so on, around an isolated single rotor in several wind speed conditions. In addition, four fitting parameters must be determined in advance to express the interaction effects between two rotors in the CO, CD, CU, and TD layouts having a specific inter-rotor space (gap). By applying the law of the momentum conservation with including the interaction effects through the modification of the pressure loss of the control volume (CV), the method could predict the power dependence of the paired rotor with the rotor gap of 1.0D on the 16 wind directions, which agreed with the CFD analysis well. However, the dependence on the inter-rotor distance (between rotor centers) should be improved. The application of the method to the three-rotor cluster showed almost the same trend as the CFD in the averaged power distribution over 12 wind directions. The method was applied to the four-rotor layouts arranged in line of parallel or tandem to the main flow to show the applicability to VAWT wind farms. The important advantages of our method include the very small number of the calculation grid (actually 400 × 400 used) and the random positioning of rotors as well as the significantly short calculation time compared with the CFD using the DFBI model. Unlike the conventional wake models, our wake model can express the acceleration regions existing on both sides of the velocity deficit. Once the necessary parameters are prepared, the averaged flow field around an arbitrary rotor layout of the VAWT cluster can be reproduced with fidelity as high as the CFD analysis. To apply our method to the problem of finding the optimal layout of rotors in a VAWT wind farm, it is necessary to increase the accuracy of the prediction and to introduce other optimization methods such as the genetic algorithm (GA) instead of the round robin.