Fast Simulation of the Flow Field in a VAWT Wind Farm Using the Numerical Data Obtained by CFD Analysis for a Single Rotor

Abstract

The effects of an increase in output power owing to the close arrangement of vertical-axis wind turbines (VAWTs) are well known. With the ultimate goal of determining the optimal layout of a wind farm (WF) for VAWTs, this study proposes a new method for quickly calculating the flow field and power output of a virtual WF consisting of two-dimensional (2-D) miniature VAWT rotors. This new method constructs a flow field in a WF by superposing 2-D velocity numerical data around an isolated single VAWT obtained through a computational fluid dynamics (CFD) analysis. In the calculation process, the VAWTs were gradually increased one by one from the upstream side, and a calculation subroutine, in which the virtual upstream wind speed at each VAWT position was recalculated with the effects of other VAWTs, was repeated three times for each arrangement with a temporal number of VAWTs. This method includes the effects of the velocity gradient, secondary flow, and wake shift as models of turbine-to-turbine interaction. To verify the accuracy of the method, the VAWT rotor power outputs predicted by the proposed method for several types of rotor pairs, four-rotor tandem, and parallel arrangements were compared with the results of previous CFD analyses. This method was applied to four virtual WFs consisting of 16 miniature VAWTs. It was found that a layout consisting of two linear arrays of eight closely spaced VAWTs with wide spacing between the arrays yielded a significantly higher output than the other three layouts. The high-performance layout had fewer rotors in the wakes of the other rotors, and the induced flow speeds generated by the closely spaced VAWTs probably mutually enhanced their output power.

1. Introduction

To realize a decarbonized society, the early introduction of various forms of renewable energy is necessary. Small wind turbines have advantages over large wind turbines, such as ease of installation; however, the high cost of power generation remains an issue. In recent years, closely spaced arrangements of small vertical-axis wind turbines (VAWTs) have attracted attention because of their power output enhancement effect. If VAWTs are effectively installed with close spacing to suit the shape of the planned site and wind conditions, not only will the amount of electricity generated per unit installation area increase, but construction, maintenance, and removal will also be more efficient, which is expected to significantly reduce power generation costs. However, because the flow field around a VAWT is more complex than that around a horizontal-axis wind turbine (HAWT), the computational cost of numerical analysis is high, and the development of methods for determining the optimal VAWT layout and simple wake models of the VAWT has not progressed sufficiently.

In a pioneering study of VAWT-WF, Rajagopalan et al. [1] performed a two-dimensional (2-D) computational fluid dynamics (CFD) analysis of several cases, including up to 21 VAWTs, assuming a 2-D laminar flow. They showed that VAWTs located in a backward oblique direction could generate higher output power than front VAWTs exposed to undisturbed wind. Whittlesey et al. [2] proposed a method that included the velocity deficit effect of the turbine wake in the potential flow, and showed that a wind farm (WF) consisting of VAWTs closely arranged like a school of fish, in which rows of VAWTs rotating in opposite directions are aligned alternately, could potentially produce greater output per unit installation area than conventional large wind turbines consisting of HAWTs. Dabiri [3] conducted field experiments with small VAWTs with a diameter of 1.2 m and demonstrated that a WF with a grid-like unit or cluster arrangement of a pair of counter-rotating VAWTs could have a higher energy density than a conventional WF with large HAWTs. In addition, he presented measurement results showing the wind-direction dependence of the output of a pair of VAWTs. Kinzel et al. [4] constructed a WF by arranging three pairs × three pairs of counter-rotating VAWT pairs, using the same 1.2 m diameter rotor as that used by Dabiri, in a grid pattern, and experimentally showed that the wake recovery in VAWT-WFs could be faster than that of HAWT-WFs. Zanforlin et al. [5] performed a 2-D CFD analysis using a turbulence model for a 2-D rotor equivalent of the actual rotor used by Dabiri (diameter: 1.2 m). Their CFD analysis showed that a pair of closely spaced counter-rotating VAWTs produced more average power than a single turbine. De Tavernier et al. [6] conducted a 2-D numerical analysis based on the panel vortex method to investigate the effects of rotor loading, rotor spacing, and wind direction on the flow field and output of a pair of VAWTs (each rotor diameter was 20 m) with close spacing. Hezaveh et al. [7] performed a three-dimensional (3-D) CFD analysis using large eddy simulation (LES) on a WF consisting of three small VAWTs as a unit. More recent research using numerical calculations for VAWT clusters has also been conducted. Sahebzadeh et al. [8] conducted a high-fidelity CFD analysis of a pair of VAWTs with one blade for each turbine, and Hansen et al. [9] performed a 2-D CFD analysis of a pair of three-blade VAWTs with a diameter of 20 m and a three-unit arrangement. In each study, the layout conditions and power-increase ratios when a WF or cluster produced more average power than a single VAWT were investigated in detail. Silva et al. [10] conducted a 2-D CFD analysis of the arrangements of three VAWTs (diameter: 3 m, three-bladed) and showed that when the three turbines were arranged in a straight line, the power output was improved compared to other arrangements that were not in a straight line. Xia et al. [11] attempted a 3-D CFD analysis of a pair of nautilus-shaped drag-type VAWTs. They also conducted simple experiments with 3-D printing models and proposed several layouts for WFs. Recently, wind tunnel experiments on clusters of VAWTs have also been conducted by several research groups. Ahmadi-Baloutaki et al. [12] performed wind tunnel experiments on configurations of two or three small VAWTs with five blades and a diameter of 0.3 m, and investigated the optimal spacing of turbines in the trio arrangements where the third turbine was located midway downstream of the counter-rotating pair. Vergaerde et al. [13,14] conducted wind tunnel experiments on a pair of small VAWTs (two blades and a diameter of 0.5 m) to confirm that the power increase was consistent with numerical predictions. They also observed the synchronization of a pair of VAWTs. Su et al. [15] used a model of an H-type VAWT with a diameter of 0.24 m to conduct wind tunnel experiments on the output of two- and three-unit clusters. Jodai et al. [16] applied miniature butterfly-type VAWT models (diameter: 50 mm) with three blades to wind tunnel experiments for six-rotor configurations consisting of a VAWT pair or trio as the unit (cluster). The experiments demonstrated the possibility of increasing the output per unit installation area by reducing the cluster spacing.

Regarding the optimal placement of wind turbines in a WF, Mosetti et al. [17] conducted research on the placement of HAWT WFs in 1994 using a genetic algorithm (GA), which is a heuristic optimization method. Numerous studies have been conducted on HAWT WFs (Gonzales et al. [18], Wagner et al. [19], Gao et al. [20,21], Sun et al. [22]). For example, Chen et al. [23] developed a 3-D greedy algorithm for optimizing the WF of HAWTs and claimed that their method has a lower computational cost than the GA. Gualtieri [24] considered the optimal HAWT configuration and optimal HAWT specifications, inputting data from 200 land-based commercial HAWTs, and used a self-organizing map (SOM) to optimize the layout of a 2 × 2 km square WF. In their study, the Jensen model was used as the wake model for HAWT. To reduce the computational cost of CFD analysis of the wake of HAWTs, the actuator disk model (ADM) and actuator line model (ALM) have been studied, which can simulate a flow field that is closer to reality than the simple Jensen model [25]. ADM regards the wind turbine rotor as a disk and represents the disk as a number of discrete points, whereas ALM represents a blade as a linear combination of discrete points and rotates some blades to simulate rotor rotation. Both ADM and ALM typically use the aerodynamic data of the blade to calculate the aerodynamic forces corresponding to the local wind speed at each discrete point, which are then introduced into the CFD analysis as external forces acting on the fluid. Recently, Malecha and Dsouza [26] proposed an ADM that provides power and thrust coefficients as functions of wind speed without using aerodynamic data.

The wake modeling of wind turbines has also progressed for HAWTs. Göçmen et al. [27] described six widely used wake models developed at the Technical University of Denmark, including the simplest and most widely used Jensen model. They described four wake-effect superposition approaches (geometric sum, linear sum, energy balance, and quadratic sum) that considered the effects of other turbines located upstream of a turbine. In a recent study on HAWT wake models, Shao et al. [28] proposed a correction to the energy conservation sum when using the Jensen model to consider the increased turbulence when the distance to the upstream turbine is short. However, there have been few studies on the optimization of a VAWT-WF and wake models. Chen and Agarwal [29] used a GA to compare the optimization of both HAWT-WFs and VAWT-WFs. They used the simple Jensen model as the wake model for VAWTs and stated that their research was a preliminary attempt. Lam and Peng [30] developed a biased wake model for a five-bladed H-shaped VAWT model (diameter: 0.3 m, height: 0.3 m) through wind tunnel experiments, and presented the results of the optimization of the virtual VAWT-WFs (diameter: 3.0 m, rated power: 3.0 kW) using a covariance matrix adaptive evolution strategy (CMA-ES). Talamalek et al. [31] investigated the effects of turbulence on the wake of a small VAWT pair (two-bladed, diameter: 0.5 m) by conducting wind tunnel experiments, and found that turbulence is beneficial for improving the power output and wake recovery of the VAWT pair. Cazzaro et al. [32] performed a heuristic optimization based on a variable neighborhood search (VNS) implemented with a Gaussian wake model with bias for the Troposkien VAWT. The optimal arrangement of the WF, including 111 VAWTs (with different rotational directions), was reported.

The flow field around a VAWT differs from that around an HAWT. Not only is there a bias (or shift) in the VAWT wake owing to rotation, but an increase in the flow velocity is also observed on both sides of the VAWT. The size and location of the velocity-increased regions change because of the induced velocity associated with rotation. Buranarote et al. [33] attempted to develop a method to mimic the flow field obtained from the 2-D CFD analysis of a VAWT. They used the potential flow as a base flow and added an equation representing the velocity-increased regions of a VAWT to a wake model (Shapiro et al. [34]) of an HAWT. The method successfully reproduced the wind-direction dependence of the output of pairs or trios of VAWTs with a certain degree of accuracy; however, the number of parameters was very large, the program became complex, and the calculation costs were somewhat high, although not as high as those of CFD.

In this study, we fundamentally change the method proposed by Buranarote et al., and propose a new method that significantly reduces the number of parameters by directly using the numerical data of the flow field obtained from a CFD analysis of an isolated single 2-D VAWT rotor and represents the WF flow field by superposing the CFD-based flow velocity data. The possibilities of the new method to significantly reduce the calculation time and improve the prediction accuracy compared to the previous method are shown. Furthermore, this method was applied to four layouts of virtual WFs consisting of 16 miniature 2-D VAWT rotors, and the characteristics and possibilities of closely spaced VAWT arrangements are demonstrated.

2. Methods

First, we provide an overview of the calculation method for the WF flow field using the method proposed in Section 2.1. In Section 2.2, we explain the model representing the interference effects between multiple VAWTs and the method for determining the characteristics of each wind turbine. In Section 2.3, the processes used to build a flow field are described.

2.1. Calculation of WF Velocity Field Using Single-VAWT CFD Data

This study aims to predict the flow field and output power of a WF consisting of multiple VAWTs installed at a site in a certain area. However, for simplicity, a 2-D flow field is assumed. In the previous method developed by Buranarote et al. [33], to reproduce the flow field around a 2-D VAWT, which is equivalent to the results obtained by CFD analysis, a unique function representing the velocity-increased regions existing at both sides of the VAWT was added to the wake model (super-Gaussian function) of Shapiro et al. [34]. This model simulates the transformation in the distribution of the velocity deficit from a top-hat type to a Gaussian type in the near-wake region. Our new method does not use functions to represent the velocity field (velocity deficit) around an isolated single VAWT but directly uses the numerical flow velocity data obtained from a CFD analysis. However, the CFD analysis is performed as an unsteady calculation (Hara et al. [35]) in advance under the condition of a specific constant upstream wind speed U_∞0, and the results are averaged to obtain the 2-D velocity field (u_CFD, v_CFD) as a steady field. Figure 1 illustrates the proposed method. In the example in Figure 1, the white circles represent the 2-D VAWT rotors, and all the rotors rotate in the counterclockwise direction. The coordinate system has an x-axis in the mainstream direction and a y-axis perpendicular to the mainstream direction. Two color maps in Figure 1 show the distributions of the u and v components of flow velocity near a single miniature VAWT rotor (diameter: 50 mm) obtained by CFD analysis under the condition U_∞0 = 10 m/s. The upstream wind speed U_∞0 when performing CFD analysis does not have to be the same as the upstream wind speed U_∞ of any WF being calculated. If there are CFD analysis data for a specific wind speed, the data can be used to calculate the WF under any upstream wind speed condition. In addition, the effects of the other rotors in the WF on the surrounding flow field of the rotor are calculated based on the prepared CFD data of a single isolated VAWT.

The rotor diameter of the VAWT to be calculated is assumed to be D. The velocity data of the prepared CFD results around a single isolated VAWT are defined as ($u_{\mathrm{C}\mathrm{F}\mathrm{D}},v_{\mathrm{C}\mathrm{F}\mathrm{D}}$). When the coordinate origin in the CFD analysis is set at the center of the VAWT rotor and the normalized coordinate system divided by the rotor diameter is expressed as ($x_{\mathrm{n}0},y_{\mathrm{n}0}$), $x_{\mathrm{n}0}$ and $y_{\mathrm{n}0}$ are defined by Equations (1) and (2), respectively:

$$x_{\mathrm{n}0}={\displaystyle \frac{x}{D}}$$

(1)

$$y_{\mathrm{n}0}={\displaystyle \frac{y}{D}}$$

(2)

The normalized velocity deficit field (${du}_{\mathrm{d}\mathrm{i}\mathrm{f}},{dv}_{\mathrm{d}\mathrm{i}\mathrm{f}}$) generated around a single isolated VAWT is defined by Equations (3) and (4) when the upstream wind speed is U_∞0.

$${du}_{\mathrm{d}\mathrm{i}\mathrm{f}}\left(x_{\mathrm{n}0},y_{\mathrm{n}0}\right)={\displaystyle \frac{1}{U_{\infty 0}}}\left\{u_{\mathrm{C}\mathrm{F}\mathrm{D}}\left(x_{\mathrm{n}0},y_{\mathrm{n}0}\right)-U_{\infty 0}\right\}$$

(3)

$${dv}_{\mathrm{d}\mathrm{i}\mathrm{f}}\left(x_{\mathrm{n}0},y_{\mathrm{n}0}\right)={\displaystyle \frac{v_{\mathrm{C}\mathrm{F}\mathrm{D}}\left(x_{\mathrm{n}0},y_{\mathrm{n}0}\right)}{U_{\infty 0}}}$$

(4)

In this study, the entire computational domain, including the WF area, was assumed to be a square with side L. The origin of the coordinate system for the WF calculation was placed at the center of the entire computational domain, and the size of the computational domain was given by the parameter ROC = L/(2D), which is the length from the center to the boundary divided by the turbine diameter. The computational grid is divided into m parts, vertically and horizontally. The grid cell size Cell_{n_D}, which is nondimensionalized by the wind turbine diameter, is defined as follows:

$${Cell}_{\mathrm{n}\_\mathrm{D}}={\displaystyle \frac{L}{mD}}={\displaystyle \frac{2ROC}{m}}$$

(5)

The CFD data (as input data) used in the calculation were limited to those in the ranges shown in Equations (6) and (7). Figure 2 shows the area of the CFD data to be input using a color map (only the u component is shown as an example).

$$-2\le x_{\mathrm{n}0}\le 10$$

(6)

$$-8\le y_{\mathrm{n}0}\le 8$$

(7)

Figure 2. Extrapolation regions [A] to [H].

Figure 2. Extrapolation regions [A] to [H].

In the proposed method, when the center of a single rotor or any rotor in the WF is set as the origin of the local coordinate system, the normalized velocity deficit field at any point P(x_np, y_np) existing in the area from [A] to [H] outside the input range of the CFD data in Figure 2 is calculated by extrapolating from the value at point B(x_nb, y_nb) on the boundary of the input range of the CFD data using Equations (8) and (9): in Equations (8) and (9), disPB is the distance between a point P and a point B, and τ_out is a parameter that represents the rate of decay of the normalized velocity deficit.

$${du}_{\mathrm{d}\mathrm{i}\mathrm{f}}\left(x_{\mathrm{n}\mathrm{p}},y_{\mathrm{n}\mathrm{p}}\right)={du}_{\mathrm{d}\mathrm{i}\mathrm{f}}\left(x_{\mathrm{n}\mathrm{b}},y_{\mathrm{n}\mathrm{b}}\right)\times e^{-disPB/{\tau }_{\mathrm{o}\mathrm{u}\mathrm{t}}}$$

(8)

$${dv}_{\mathrm{d}\mathrm{i}\mathrm{f}}\left(x_{\mathrm{n}\mathrm{p}},y_{\mathrm{n}\mathrm{p}}\right)={dv}_{\mathrm{d}\mathrm{i}\mathrm{f}}\left(x_{\mathrm{n}\mathrm{b}},y_{\mathrm{n}\mathrm{b}}\right)\times e^{-disPB/{\tau }_{\mathrm{o}\mathrm{u}\mathrm{t}}}$$

(9)

Figure 2 shows three examples (points P, P′, and P″) where point P exists in different regions. When a point P is in the regions [A] to [D] shown in Equations (10)–(13) (e.g., point P or point P′ in Figure 2), point B is a point on the boundary of the input range of the CFD data (e.g., point B or point B′ in Figure 2).

$$\mathrm{region}\mathrm{A}:x_{\mathrm{n}0}<-2\mathrm{and}-8\le y_{\mathrm{n}0}\le 8$$

(10)

$$\mathrm{region}\mathrm{B}:x_{\mathrm{n}0}>10\mathrm{and}-8\le y_{\mathrm{n}0}\le 8$$

(11)

$$\mathrm{region}\mathrm{C}:-2\le x_{\mathrm{n}0}\le 10\mathrm{and}{y}_{\mathrm{n}0}<-8$$

(12)

$$\mathrm{region}\mathrm{D}:-2\le x_{\mathrm{n}0}\le 10\mathrm{and}{y}_{\mathrm{n}0}>8$$

(13)

When a point P is in the regions [E] to [H] shown in Equations (14)–(17) (e.g., point P′′ in Figure 2), point B is a point on the boundary of region [C] or [D] (e.g., point B″ in Figure 2).

$$\mathrm{region}\mathrm{E}:x_{\mathrm{n}0}<-2\mathrm{and}{y}_{\mathrm{n}0}<-8$$

(14)

$$\mathrm{region}\mathrm{F}:x_{\mathrm{n}0}>10\mathrm{and}{y}_{\mathrm{n}0}<-8$$

(15)

$$\mathrm{region}\mathrm{G}:x_{\mathrm{n}0}<-2\mathrm{and}{y}_{\mathrm{n}0}>8$$

(16)

$$\mathrm{region}\mathrm{H}:x_{\mathrm{n}0}>10\mathrm{and}{y}_{\mathrm{n}0}>8$$

(17)

Assuming that N VAWTs of the same size exist in a WF, the proposed method assumes that the flow velocity components u and v at any position (x, y) in the WF can be calculated using Equations (18) and (19), respectively:

$$u\left(x,y\right)=U_{\infty }+\sum _{j=1}^N{U}_{\mathrm{F}}\left(j\right){du}_{\mathrm{d}\mathrm{i}\mathrm{f}}\left(x_{\mathrm{n}}\left(j\right),y_{\mathrm{n}}\left(j\right)\right)$$

(18)

$$v\left(x,y\right)=\sum _{j=1}^N{U}_{\mathrm{F}}\left(j\right){dv}_{\mathrm{d}\mathrm{i}\mathrm{f}}\left(x_{\mathrm{n}}\left(j\right),y_{\mathrm{n}}\left(j\right)\right)$$

(19)

where U_∞ is the upstream wind speed and U_F(j) is the virtual upstream wind speed of the j-th VAWT located at position (x_j, y_j). The decision on the value of U_F(j) is explained in the following section. The arguments for the normalized velocity deficits x_n(j) and y_n(j) are given by Equations (20) and (21), respectively.

$$x_{\mathrm{n}}(j)={\displaystyle \frac{x-x_j}{D}}$$

(20)

$$y_{\mathrm{n}}(j)={\displaystyle \frac{y-y_j}{D}}$$

(21)

2.2. Interference Model and Wind Turbine Characteristics

In the proposed method, the virtual upstream wind speed U_F(j) of each VAWT is predicted using Equation (22).

$$U_{\mathrm{F}}\left(j\right)=u_{\mathrm{a}\mathrm{v}\mathrm{e}}\left(j\right)+∆u1\left(j\right)-∆u2\left(j\right)$$

(22)

where u_ave(j) is the u component of the average wind speed at the position where the j-th VAWT is installed before installation. As for the VAWTs other than the subject VAWT, some were installed and others were not, and their situations differed in the calculation process. In this study, u_ave(j) is defined as the average of the u component at 11 points on a line segment of the same length as the diameter in the y direction, which includes the center of the j-th VAWT. ∆u1(j) and ∆u2(j) in Equation (22) represent the influence of the surrounding VAWTs on the j-th VAWT and are defined in Equations (23) and (24), respectively.

$$∆u1\left(j\right)=k_{\mathrm{c}}\left\{1-E_{\mathrm{s}\mathrm{s}}{S}_{\mathrm{r}\mathrm{d}}(j)k_{\mathrm{u}}{\displaystyle \frac{du}{dy}}\right\}\left\{1-S_{\mathrm{r}\mathrm{d}}\left(j\right)k_{\mathrm{v}}{v}_{\mathrm{a}\mathrm{v}\mathrm{e}}\left(j\right)\right\}$$

(23)

$$∆u2\left(j\right)=u_{\mathrm{a}\mathrm{v}\mathrm{e}}\left(j\right)\times e^{-\left\{{\displaystyle \frac{\left|v_{\mathrm{a}\mathrm{v}\mathrm{e}}\left(j\right)\right|}{k_1}}+{\displaystyle \frac{d_{\mathrm{m}\mathrm{i}\mathrm{n}}}{k_{2}D}}\right\}}$$

(24)

v_ave(j) in the above equations is the v component of the average wind speed at the position where the j-th VAWT is to be installed before installation. In this study, v_ave(j) is defined as the average of the v component at 11 points on a line segment of the same length as the diameter in the y direction, which includes the center of the j-th VAWT.

∆u1(j) in Equation (23) is introduced as the increment in the virtual upstream wind speed U_F(j) to represent the resulting increase in the power output of the j-th VAWT due to the effects of other VAWTs in the vicinity, which cause a velocity gradient du⁄dy or a secondary flow component v_ave(j) at the position where the j-th VAWT is installed. In this study, wind turbines with a large solidity σ are considered to be of the drag type, and those with a small solidity σ are considered to be of the lift type. Figure 3 illustrates the predicted effects of the velocity gradient du/dy and the secondary flow component v on the output power of each type of wind turbine.

As shown in Figure 3, the drag- and lift-type outputs are expected to respond differently to the velocity gradient, whereas it is assumed that both types of wind turbines provide the same response to secondary flow. In Equation (23), the effects of the velocity gradient and secondary flow are not treated independently, but are given as the product of both, because the results of a preliminary comparison with the CFD analysis of a pair of VAWTs, described below, showed that Equation (23) was better than the case in which the two effects were treated independently. Although Figure 3 shows only the cases of counterclockwise (CCW) rotation, the parameter S_rd(j), which represents the rotation direction of the j-th VAWT and is defined by Equations (25) and (26), is introduced to properly handle clockwise (CW) rotation.

$$S_{\mathrm{r}\mathrm{d}}\left(j\right)=1\mathrm{for}\mathrm{CCW}\mathrm{rotation}$$

(25)

$$S_{\mathrm{r}\mathrm{d}}\left(j\right)=-1\mathrm{for}\mathrm{CCW}\mathrm{rotation}$$

(26)

The results of the CFD analyses previously conducted in our laboratory on four differently sized rotor pairs of VAWTs with different solidities are shown in Figure 4.

Table 1. Specification and calculation conditions of four rotor pairs.

Symbol	L	ML	M	S
Solidity: σ[–]	0.025	0.057	0.095	0.382
Diameter: D [m]	10	5	2	0.05
Number of blade: B [–]	2	3	3	3
Chord length: c [m]	0.4	0.3	0.2	0.02
gap/D [–]	0.2	0.2	0.2	0.2
Wind velocity: U∞ [m/s]	6	10	10	10

summarizes the specifications and calculation conditions of the four cases (L, ML, M, and S).

Although the details are not shown, Figure 4b shows the solidity dependence of the three types of paired VAWT arrangements (counter-down (CD), counter-up (CU), and co-rotation (CO)) shown in Figure 4a, with the vertical axis representing the output power averaged between the two rotors (P_{n_ave}), which is normalized by dividing by an isolated single-rotor power. The horizontal axis represents the solidity (σ = Bc/(πD)) of each VAWT rotor with a logarithmic coordinate axis. CFD analyses were performed only under the condition of gap/D = 0.2. As shown in Figure 4b, when the solidity is large (S-size), the output of the CD arrangement is larger than that of the CU arrangement, and the CO arrangement is biased but lies between the CD and CU arrangements. However, when the solidity was small (L-size), the output of the CU arrangement was larger than that of the CD arrangement, indicating the opposite tendency. This tendency has also been shown in the results of numerical analyses conducted by other researchers. For example, in the study by Zanforlin et al. [5], CD > CU results were obtained for a VAWT pair with high solidity (B = 3, c = 0.128 m, D = 1.2 m, σ = 0.102), whereas, in the research of De Tavernier et al. [6], CU > CD results were obtained for a VAWT pair with low solidity (B = 2, c = 1 m, D = 20 m, σ = 0.0318) under the condition of gap/D = 0.2.

Referring to Figure 4b, the reversal of the CD and CU outputs occurs around the solidity (σ = 0.057) of the ML-size VAWT pair. Therefore, in the present model, we define a parameter E_ss in Equation (27) that expresses the effects of solidity, including the reversal of the CD-CU output power. This parameter is multiplied by a factor indicating the effects of the velocity gradient in ∆u1(j) of Equation (23), which is expected to produce different responses between the drag type and lift type.

$$E_{\mathrm{s}\mathrm{s}}={\displaystyle \frac{\sigma -0.057}{0.325}}$$

(27)

k_c, k_u, and k_v in Equation (23) are the fitting parameters that were determined in this study by comparison with a previous CFD analysis of a pair of VAWTs (Hara et al. [35]).

Next, ∆u2(j) in Equation (24) is explained. ∆u2(j) is a correction introduced from a comparison of a CFD analysis of paired VAWTs (Hara et al. [35]) and the predictions for 16 wind directions by the present method. In the preliminary comparison, it was found that when a VAWT was located diagonally behind another VAWT located upstream, the method tended to predict a higher output power than the CFD analysis. Because the state of overprediction corresponded to a region where the absolute value of the secondary velocity was small to some extent, the virtual upstream wind speed was reduced by the amount calculated using Equation (24). Note that d_min in Equation (24) is the distance to the nearest VAWT upstream of the j-th VAWT (center-to-center), and if this value is large, the amount of correction by ∆u2(j) will be small. k₁ and k₂ in Equation (24) are the fitting parameters determined by comparison with the CFD analysis for a pair of VAWTs.

Induced velocities occur around the VAWTs owing to the rotation of the rotors. Consequently, the wake of the VAWT shifts in a direction perpendicular to the main flow. The CFD analysis results (input data) for a single isolated VAWT include this shift; however, the amount of wake shift changes owing to the influence of the other VAWTs. To consider this effect, in the proposed method, argument y_n(j) of the normalized velocity deficit for the j-th VAWT introduced in Equation (21) is replaced by Equation (28).

$$y_{\mathrm{n}}(j)={\displaystyle \frac{y-y_j-k_{\mathrm{s}}{v}_{\mathrm{a}\mathrm{v}\mathrm{e}}\left(j\right)\left\{1-e^{-x_{\mathrm{n}}\left(j\right)}\right\}}{D}} (x_{\mathrm{n}}\left(j\right)>0)$$

(28)

The above equation is applied only to the wake region where x_n(j) is positive, and Equation (21) is used upstream. Equation (28) assumes that a shift proportional to the secondary flow component v_ave(j) at the center of the j-th VAWT occurs in the wake region; however, an exponential function is introduced to prevent a discontinuous shift near the rotor. k_s is a fitting parameter that adjusts the amount of the shift.

In this method, once the virtual upstream wind speed U_F(j) is obtained using Equation (22), it is converted into the characteristics of the j-th VAWT in the WF based on the pre-entered power curve and the relationship between the operating rotor speed and wind speed to obtain the prediction results. At this time, the input characteristics, such as the power curve, may be experimental data instead of theoretically predicted values. In this study, we focused on the 2-D rotor of a miniature VAWT with a diameter of 50 mm, for which we previously conducted a CFD analysis. Therefore, the input power curve and the relationship between the rotation speed and wind speed are given as a curve and a straight line passing through the power (P_s__CFD = 177 mW) and rotation speed (N_{s_CFD} = 3483 rpm) obtained in a previous CFD analysis of an isolated single 2-D rotor under the condition of a wind speed U_∞= 10 m/s. Figure 5 shows the power curve and the relationship between the rotation speed and wind speed, which were used as input data in this study.

The target wind turbine in this study is the same as that assumed in reference [35], and the torque characteristics have been obtained as the results of CFD analysis for four upstream wind speeds U_∞ = 6, 8, 10, and 12 m/s, as shown in Figure 2 of reference [35] (see Appendix A). In this study, the load torque characteristics are defined as a curve that is proportional to the square of the angular velocity ω passing through the state of 95% of the maximum output at U_∞ = 10 m/s, and are given by Equation (1) in reference [35] or Equation (A1) in Appendix A. The rotation speed N_s [rpm] and power P_s [mW] of the isolated single wind turbine shown in Figure 5 are approximated by Equations (29) and (30), respectively, as functions of the upstream wind speed U_∞.

$$N_{\mathrm{s}}=376.34U_{\mathrm{\infty }}-280.58 \left[\mathrm{rpm}\right]$$

(29)

$$P_{\mathrm{s}}=0.2459U_{\mathrm{\infty }}^3-0.8344U_{\mathrm{\infty }}^2+1.4757U_{\mathrm{\infty }}-0.2466 \left[\mathrm{mW}\right]$$

(30)

From Equation (29), when the wind speed is below 0.7455 m/s, the rotational torque of the single turbine obtained by the CFD analysis is smaller than the load torque. For simplicity, in this study, it is assumed that the target wind turbine cannot rotate at wind speeds below 2 m/s, and the input power curve data are limited to U_∞ = 2 m/s or more.

The proposed method was developed as an engineering model with low computational costs that can be applied to explore the optimal layout of VAWTs. Therefore, ideally, it is desirable to fix the upstream wind speed for the prepared CFD analysis to a single value of the rated wind speed of the target wind turbine or the average wind speed of the site (U_∞ = 10 m/s in this study). The extent to which this method remains reliable when the upstream wind speed changes must be confirmed by conducting additional CFD analyses for various configurations; however, this is a topic for future work.

2.3. Flow Field Construction Process

In the previous section, we explained the outline of the analytical formula used in the calculations of the proposed method; however, the calculation procedure affected the results. In this method, the flow field of a WF with multiple VAWTs was constructed by placing the turbines individually in the flow field, starting with the furthest upstream and gradually constructing the entire flow field. Figure 6 illustrates the case in which the total number of VAWTs in the WF is N = 4. As shown in Figure 6a, initially, only one VAWT is placed at the most upstream position (temporary number of VAWTs: N_temp = 1), and a flow field proportional to the CFD analysis obtained by substituting U_∞ into U_F(1) in Equations (18) and (19) is assumed. Using the assumed flow field, the average wind speed components (u_ave, v_ave) at the location where the second VAWT from the upstream was to be installed were calculated, and the virtual upstream wind speed U_F(2) was obtained from Equation (22). Once U_F(2) is determined, temporarily assuming N = N_temp = 2 in Equations (18) and (19), the flow field, including the two VAWTs in the WF, is calculated by superimposing the effects of each VAWT. Using the resultant flow field, the average wind speed (u_ave, v_ave) at the third VAWT position upstream was calculated, and the virtual upstream wind speed U_F(3) was obtained using Equation (22). Thus, the number of VAWTs was gradually increased to determine the virtual upstream wind speed, U_F(j), for all the wind turbines. Finally, using the power curve data shown in Figure 5, all U_F(j) were converted into the output power and rotation speed of each VAWT to obtain a prediction result.

Although the process of building the flow field is described in Figure 6, in the case where some VAWTs have the same x coordinates, there is a degree of freedom in the ordering from the upstream side. In addition, the effects of the downstream turbines on the upstream turbines must be considered. Therefore, in the proposed method, although the calculation cost increases, once all virtual upstream wind velocities U_F(j) for the temporary wind turbine number (N_temp) have been obtained for the first time, several iterations of the calculation routines, as shown in the example in Figure 7, are inserted, and the calculations are repeated until the flow field for the temporary number state of the VAWTs converges. Figure 7 shows the case in which N_temp = 3. After the virtual upstream wind speed U_F(3) is calculated for the first three VAWTs from upstream, each VAWT turbine is removed individually from upstream, and the virtual upstream wind speed U_F(j) (j = 1, 2, 3) is recalculated. This recalculation process was performed for all intermediate states, and the flow field was gradually constructed. Although this iterative calculation is a drawback of the proposed method, it is important for constructing a reliable flow field.

3. Results and Discussion

In this section, we apply the proposed method (the new model) to specific clusters of VAWTs and demonstrate their accuracy and features. Moreover, we apply the method to a WF consisting of many miniature VAWT rotors in specific arrangements to illustrate the expected output power of the WF. In Section 3.1, we compare the results of a previous CFD analysis with the predictions calculated using the new model regarding the 16-wind-direction dependence of the output of paired VAWTs. In Section 3.2, the CFD analysis results of a four-turbine cluster are compared with those of the proposed method. In Section 3.3, the method is applied to four layouts of a virtual WF consisting of 16 VAWTs, and the annual power generation for each layout is forecasted. Section 3.4 describes the computational cost of the proposed method.

3.1. A Comparison of Both the Results of CFD and the Proposed Method on Paired VAWTs

Figure 8a,b show the flow field calculated by applying the proposed method to a pair of miniature 2-D VAWTs with each diameter of D = 50 mm rotating in the counterclockwise direction in a CO configuration. The miniature VAWT considered in this study had three blades with an NACA 0018 airfoil cross section and the S-size specifications listed in Table 1. The performance of an isolated single turbine was assumed to be a power curve, as shown in Figure 5. In Figure 8, the upstream wind speed is U_∞ = 10 m/s, and the gap length between two rotors is gap = 25 mm (gap/D = 0.5). The angle θ, which represents the wind direction, is defined in Figure 8a. In the calculation program, the mainstream always flows from left to right, and a condition with different wind direction is created by rotating the arrangement of VAWTs counterclockwise by an angle θ. Figure 8a shows the state where θ = 22.5°, and Figure 8b shows the state where θ = 337.5° (−22.5°). In Figure 8a,b, the distribution of the u component of the flow velocity is shown only in the vicinity of the paired VAWTs. However, the entire calculation domain is a square region of L × L = 1 × 1 m. The number of grids in the vertical and horizontal directions was m = 400 and the normalized grid size was Cell_{n_D} = 0.05. Figure 8c compares the variation in the output power of Rotor 1 (R1) when the wind direction was changed to 16 directions with the results of a previous CFD analysis (Hara et al. [35]). The radar graph shows the normalized output P_R1 of R1 divided by the output of a single wind turbine P_SI = 177 mW (at U_∞ = 10 m/s). Similarly, Figure 8d shows the 16-direction dependence of the normalized power (P_ave/P_SI) of the average power P_ave for R1 and R2 and is compared with previous CFD analysis results. As shown in Figure 8c, the increase in the output power of Rotor 1 at θ = 180° by the induced velocity caused by Rotor 2 and the decrease at θ = 270° resulting from the wake effects of Rotor 2 are also predicted by the new method.

Figure 9 shows the results of performing the same calculations as in Figure 8 for an inverse-rotation (IR) arrangement in which the rotation directions are opposite to each other. However, Figure 9a is for the CD arrangement with θ = 0°, and Figure 9b corresponds to the CU arrangement with θ = 180°.

In this study, the fitting parameters in the proposed method were determined mainly by comparison with the CFD analysis for the 16-wind-direction output distribution shown in Figure 8 and Figure 9. In all calculations presented in this paper, we set k_c = 0.000103, k_u = 5.4, k_v = 10.4, k_s = 0.012, k₁ = 1.2, and k₂ = 0.64. The τ_out introduced in Equations (8) and (9) does not have any effect when the VAWTs are close to each other in a small area, as in Figure 8 and Figure 9. In this study, the value of τ_out was determined to be 10 from the results of the tandem arrangement of four VAWTs described later. The iterative calculation routine performed in each state, consisting of a temporary number of VAWTs described in Figure 7, was fixed at three times, at which convergence was confirmed.

As shown in Figure 8c,d and Figure 9c,d, the degree of agreement (or, in other words, the error) between the CFD analysis results and predictions using the new model varies depending on the wind direction. Even if the values of the aforementioned fitting parameters are fixed, the degree of match varies depending on the arrangement of the target VAWT pairs, gap length, the size of the entire calculation domain, and normalized grid size (resolution). Therefore, in this study, the difference between results by CFD analysis and the new model for the output of R1 in the wind direction θ_i is defined as Equation (31), and the average value of the difference ∆P_{n_R1}(θ_i) over 16 wind directions is defined as Equation (32), and the degree of agreement (error) is quantified and compared using those values. Similarly, the difference between the CFD and our method of the average power P_ave of R1 and R2 for each wind direction θ_i and the average value over the 16 wind directions is defined by Equations (33) and (34), respectively.

$${∆P}_{\mathrm{n}\_\mathrm{R}1}\left({\theta }_i\right)={\displaystyle \frac{\left|P_{\mathrm{R}1}\left({\theta }_i\right)-P_{\mathrm{R}1\_\mathrm{C}\mathrm{F}\mathrm{D}}\left({\theta }_i\right)\right|}{P_{\mathrm{S}\mathrm{I}}}}$$

(31)

$${Err}_{\mathrm{n}\_\mathrm{R}1}={\displaystyle \frac{1}{16}}\sum _{i=1}^{16}∆P_{\mathrm{n}\_\mathrm{R}1}\left({\theta }_i\right)\times 100\left[\%\right]$$

(32)

$${∆P}_{\mathrm{n}\_\mathrm{a}\mathrm{v}\mathrm{e}}\left({\theta }_i\right)={\displaystyle \frac{\left|P_{\mathrm{a}\mathrm{v}\mathrm{e}}\left({\theta }_i\right)-P_{\mathrm{a}\mathrm{v}\mathrm{e}\_\mathrm{C}\mathrm{F}\mathrm{D}}\left({\theta }_i\right)\right|}{P_{\mathrm{S}\mathrm{I}}}}$$

(33)

$${Err}_{\mathrm{n}\_\mathrm{a}\mathrm{v}\mathrm{e}}={\displaystyle \frac{1}{16}}\sum _{i=1}^{16}∆P_{\mathrm{n}\_\mathrm{a}\mathrm{v}\mathrm{e}}\left({\theta }_i\right)\times 100\left[\%\right]$$

(34)

The differences Err_{n_R1} and Err_{n_ave} between the CFD and model predictions defined in Equations (32) and (34) when the grid number is m = 400 are shown in Figure 10a,b. Six cases are shown with gap lengths of 25, 50, and 100 mm in the CO and IR configurations. For each case, the normalized cell size, Cell_{n_D}, was set to 0.05, 0.1, 0.2, 0.4, and 0.75 (corresponding to an ROC of 10, 20, 40, 80, and 150), as shown in the color-coded bar graphs. Similarly to Figure 10, the values of Err_{n_R1} and Err_{n_ave} when the number of grids is m = 600 are shown in Figure 11. The normalized cell sizes, Cell_{n_D}, were set to 0.0333, 0.0667, 0.1333, 0.2667, and 0.5 (corresponding to ROC = 10, 20, 40, 80, and 150).

Figure 10 and Figure 11 show that the normalized errors, Err_{n_R1} and Err_{n_ave}, decrease as the gap increases. Err_{n_ave} is smaller than Err_{n_R1} because the normalized error of the output of R1 or R2 alone is partially canceled out owing to the symmetry in the average error of the two rotors. Although this was not observed in all cases, there was a tendency for errors to increase as the normalized cell size increased. In addition, there was little difference in the trends owing to the differences in the number of grids m.

3.2. Comparison of the Predictions by CFD and the Proposed Model on a Four-VAWT Cluster

Previous predictions by CFD analysis (Buranarote et al. [33]) of a four-VAWT cluster, in which all the rotors rotated in the same direction and were arranged in line with a gap of 150 mm (gap/D = 3), were compared with predictions using the proposed method. In subsequent model calculations, we set m = 400 and Cell_{n_D} = 0.1 (ROC = 20). Figure 12a is the same as Figure 20a in the paper by Buranarote et al. and shows the unsteady flow field (u-component distribution) obtained by CFD analysis. The upstream wind speed is U_∞ = 10 m/s. The values shown in the figure for P_R1, P_R2, P_R3, and P_R4 are the predicted output power values for the VAWT rotors (mW). The flow field of the four VAWTs arranged in parallel, as shown in Figure 12b, was predicted using the proposed model. The VAWT output predicted by the model was in good agreement with the CFD analysis, and the errors, which were calculated as the difference between the predicted outputs of the CFD and the model divided by the output of an isolated single VAWT (P_SI), were 0.73%, 2.1%, 1.1%, and 8.6% for R1 to R4, respectively.

Figure 13 shows a similar comparison for a four-VAWT tandem arrangement, where Figure 13a is the same as Figure 20b in Buranarote et al. [33]. The difference between the CFD and model predictions for the tandem arrangement was somewhat larger; however, the differences based on the output P_SI of a single VAWT were 0.11%, 15.5%, 3.8%, and 5.1% for R1 to R4, respectively. The error for R1, the most upstream VAWT, is small because, as mentioned above, the value of τ_out was determined so that the output of the lead turbine in this tandem arrangement of four VAWTs matches between the CFD and the model predictions.

3.3. Application to Virtual 16-VAWT Wind Farms

3.3.1. Flow Field Predictions and WF Output Dependence on 16 Wind Directions

In this section, we apply the method to a virtual WF consisting of 16 miniature 2-D VAWTs. The area where the VAWTs were placed was assumed to be a square with dimensions of 500 × 500 mm. The following four layouts of VAWTs were assumed:

• CO–4 × 4 layout
• CO–8–pairs layout
• IR–8–pairs layout
• CO–8 × 2 layout

The flow field and output of each VAWT obtained when the wind direction was θ_i = 0° are shown in Figure 14a, Figure 15a, Figure 16a and Figure 17a, in which the details of each layout are also illustrated. The upstream wind speed was assumed to be U_∞ = 10 m/s. Figure 14b, Figure 15b, Figure 16b and Figure 17b show the distribution of the average power of the WF, which was normalized by dividing the WF average power P_16ave over 16 wind directions by an isolated single VAWT power P_SI. For all layouts, the entire computational domain was a square of L × L = 2 × 2 m, and the values of the other parameters were the same as those in the four-VAWT case described in Section 3.2.

Comparing the results of four layouts, only in the case of the CO–8 × 2 layout of Figure 17 did the WF average power P_16ave exceed the single VAWT power P_SI in the wind-direction range of –22.5° ≤ θ_i ≤ 22.5° or 157.5° ≤ θ_i ≤ 202.5°. For the other three layouts, P_16ave < P_SI for all wind directions. In the case of Figure 17a, the eight VAWTs are arranged in a line perpendicular to the mainstream (parallel arrangement) with very narrow spacing (gap1/D = 0.286), so many of the turbines (actually, half the number) in the WF are not in the wake, and the effects of increasing the output by the induced velocity are also significant. In particular, R1 turbine, which is located at the bottom end of the upwind parallel arrangement, had an exceptionally high output. A similar power output improvement was reported by Mereu et al. [36] in a CFD analysis of a parallel arrangement of eight Savonius turbines; although the turbine spacing dependence differs owing to different solidities, their study and the results of our study may be due to the same effects.

3.3.2. Annual Energy Production Prediction

The proposed method assumes the existence of 2-D flow field data for an isolated single VAWT analyzed by CFD at a specific wind speed U_∞0. Assuming the availability of the measurement data of wind-direction occurrence probability (WDP), the probability of obtaining wind direction in a certain wind direction θ_i among the 16 direction segments is defined as WDP(θ_i) in this study. For simplicity, we assume that the wind speed occurrence rate is independent of wind direction and is given by a Rayleigh distribution; the probability density function of any wind speed V according to the Rayleigh distribution is expressed as f_R(V). Once the average output P_ave(U_∞0, θ_i) [W] of a WF consisting of N VAWTs at a specific wind speed U_∞0 and a wind direction θ_i is obtained by the proposed method, the annual energy production (AEP), AEP(U_∞0) [kWh], in the bin at the specific wind speed U_∞0 when the wind speed bin width is ∆V = 1 m/s, is calculated using Equation (35). The constant 8.76 is 60 × 60 × 24 × 365/(3.6 × 10⁶).

$$AEP\left(U_{\infty 0}\right)=\sum _{i=1}^{i=16}{P}_{\mathrm{a}\mathrm{v}\mathrm{e}}\left({U_{\infty 0},\theta }_i\right)\times N\times WDP\left({\theta }_i\right)\times f_{\mathrm{R}}\left(U_{\infty 0}\right)\times 8.76$$

(35)

If the cut-in wind speed is V₂ = 2 m/s and the cut-out wind speed is V₂₀ = 20 m/s, the AEP of the WF can be easily predicted using the previously calculated average power P_ave (U_∞0, θ_i) or AEP(U_∞0) using Equation (36) or Equation (37) without any additional flow field calculations.

$$AEP={\sum }_{j=2}^{j=20}{\sum }_{i=1}^{i=16}{P}_{\mathrm{a}\mathrm{v}\mathrm{e}}\left({U_{\infty 0},\theta }_i\right)\times {\displaystyle \frac{P_{\mathrm{S}\mathrm{I}}\left(V_j\right)}{P_{\mathrm{S}\mathrm{I}}\left(U_{\infty 0}\right)}}\times N\times WDP\left({\theta }_i\right)\times f_{\mathrm{R}}\left(V_j\right)\times 8.76$$

(36)

$$AEP={\sum }_{j=2}^{j=20}{\displaystyle \frac{P_{\mathrm{S}\mathrm{I}}\left(V_j\right)}{P_{\mathrm{S}\mathrm{I}}\left(U_{\infty 0}\right)}}\times {\displaystyle \frac{f_{\mathrm{R}}\left(V_j\right)}{f_{\mathrm{R}}\left(U_{\infty 0}\right)}}\times AEP\left(U_{\infty 0}\right)$$

(37)

Table 2. AEPs for four layouts consisting of 16 VAWTs for uniform wind-direction distribution.

Layout	AEP [kWh]	AEP/Area [kWh/m2]	AEP/AEP16SI
CO–4 × 4	25.857	103.428	0.673
CO–8–pairs	21.725	86.900	0.565
IR–8–pairs	22.121	88.484	0.576
CO–8 × 2	31.083	124.332	0.809

shows the AEPs calculated for the four layouts of 16 VAWTs shown in Section 3.3.1, assuming a hypothetical uniform wind-direction distribution (WDP(θ_i) = 0.0625) where the wind-direction occurrence probability is the same in all 16 directions. The wind speed distribution was assumed to be Rayleigh with an average wind speed of 10 m/s. Table 2 lists the AEP per unit area and AEP based on the total AEP of the 16 isolated VAWTs. Even in the CO–8 $\times$ 2 layout, which had the highest AEP among the four layouts, the AEP was only 80.9% of the total AEP of the 16 isolated VAWTs.

Figure 18 shows an example of the nonuniform wind-direction occurrence probability, where the prevailing wind direction appears in the direction of 157.5°. The probability of the wind-direction occurrence was measured at the Arid Land Research Center (ALRC) of Tottori University. For the four 16-VAWT layouts shown in Section 3.3.1, in a uniform wind-direction distribution, the CO–4 × 4 layout showed the maximum average WF power in four wind directions, θ_i = 67.5°, 157.5°, 247.5°, 337.5° (−22.5°), the CO–8–pairs layout showed the maximum average power in two wind directions, θ_i = 157.5° and 337.5° (−22.5°), the IR–8–pairs layout showed the maximum average power in one wind direction, θ_i = 337.5° (−22.5°), and the CO–8 × 2 layout showed the maximum average power in two wind directions, θ_i = 0° and 180°. If the condition that maximizes the average WF output in a uniform wind-direction distribution is aligned with the prevailing wind direction in a nonuniform wind-direction distribution, the maximum AEP of the WF must be obtained. Assuming an extremely nonuniform wind-direction distribution, as shown in Figure 18, and a Rayleigh distribution with an average wind speed of 10 m/s and matching the maximum output direction of each layout of the 16 VAWTs to the prevailing wind direction, the AEPs were calculated and are tabulated in

Table 3. AEPs for four layouts consisting of 16 VAWTs for nonuniform wind-direction distribution shown in Figure 18.

Layout	AEP [kWh]	AEP/Area [kWh/m2]	AEP/AEP16SI
CO–4 × 4	27.603	110.412	0.718
CO–8–pairs	25.868	103.472	0.673
IR–8–pairs	26.170	104.68	0.681
CO–8 × 2	46.508	186.032	1.210

.

In all layouts, when the maximum average output direction was aligned with the prevailing wind direction, the output increased compared to the case of a uniform wind-direction distribution. By calculating the percentage increase in AEP for each layout based on the total AEP of the 16 isolated VAWTs, the results were 6.7% for CO–4 × 4, 19.1% for CO–8–pairs, 18.2% for IR–8–pairs, and 49.6% for CO–8 × 2. The CO–8 × 2 layout shows an outstanding increase rate, and as a result, it is interesting that we can expect an increase in the WF output of approximately 20% compared with 16 isolated VAWTs.

3.4. Calculation Cost of the Developed Simulation Program

Table 4. The calculation cost of a program implemented using the proposed method.

Layout	N	Time (1-dir.) [s]	Time (16-dir.) [s]
SI	1	0.085	0.385
CO–50 mm	2	0.220	1.107
4–rotor parallel	4	1.079	2.93
CO–4 × 4	16	54.425	114.739

lists the computational costs of the simulation program for calculating the WF flow field and the output powers of the VAWTs using the proposed method. Four layouts were considered: one, two, four, and sixteen VAWTs. In this computational cost analysis, calculations were performed under the same conditions: ROC = 20, m = 400, and Cell_{n_D} = 0.1. The program was written in C#, and the PC used for the calculations had a 10-core Intel Core (TM) i9-10900 2.81 GHz processor. Table 4 lists the time required for the calculation of one wind direction and the time required for the continuous calculation of sixteen wind directions. However, to calculate the 16 wind directions, parallel calculations are performed using the “Parallel.For” method in the C# language. As shown in Section 2.3, in this calculation method, the flow field is constructed by increasing the number of VAWT rotors one by one, and the sub-calculation routine is repeated three times for each temporary number state of the rotors. Therefore, as shown in the calculation results for one wind direction in Table 4, the calculation time increased at a rate of approximately N(N + 1) with increasing N, the number of VAWTs. However, owing to the effect of parallel computing, the calculation time for 16 wind directions was approximately twice that for one wind direction, even in the case of 16 VAWTs.

4. Conclusions

The ultimate goal of this study was to determine the optimal layout of VAWTs for a WF. As preliminary research, we proposed a method to predict the flow field and output of the VAWT-WF in a short time. The proposed method assumes the existence of numerical data obtained by CFD analysis of the velocity field around an isolated single two-dimensional VAWT rotor. By superimposing the data of all VAWT rotors in a WF, a significant reduction in calculation time and an improvement in accuracy were achieved compared with previous algebraic wake modeling methods. The turbine-to-turbine interaction model in this method includes the effects of velocity gradient, secondary flows, and wake shift. We also proposed a parameter (E_ss) to represent the effect of solidity, but further research is needed to extend this to solidities other than the miniature wind turbine (σ = 0.382) targeted in this study.

This method is also unique in the WF flow field building process. When constructing the flow field of a WF with a total VAWT number of N, the temporary number of VAWTs is gradually increased from one to N, while at each temporary rotor number state, one turbine is removed to re-calculate the virtual upstream wind speed U_F of the removed turbine repeatedly. Therefore, the computational cost for a particular wind direction is approximately proportional to N (N + 1). However, the calculation time for 16 wind directions was reduced by parallelizing the calculations.

Comparing the previous CFD analysis of a pair of VAWTs with the results predicted by the proposed method, the normalized error tended to decrease as the gap length between the two turbines increased. In addition, the error tended to increase as the normalized cell size increased; however, there was little difference in the trends owing to the differences in the number of grids.

A four-VAWT linear arrangement (parallel and tandem) was also compared between the CFD and model analyses, and showed relatively good agreement.

Current model parameter values vary depending on the size, structure, and solidity of the wind turbine. If an operation curve is determined for a single wind turbine, the model parameters can be obtained by conducting a CFD analysis for a minimum of 17 cases, including eight wind directions for each (from which the power dependence of each rotor pair on the 16 wind directions is derived from symmetry) of representative closely spaced rotor pairs of CO and IR arrangements (e.g., gap = 0.5D) and one wind direction for a representative tandem layout (e.g., four rotors with gap = 3D). In the future, if a CFD analysis of various wind turbines is performed to clarify the dependence of each parameter on the turbine size and solidity, it may be possible to express the parameter in a universal mathematical formula.

When this method was applied to four virtual WFs with 16 2-D miniature VAWTs arranged, it was found that a layout (CO–8 × 2) in which the spacing between neighboring turbines (gap1/D = 0.286) was shortened and eight turbines were arranged in a line in two rows that were sufficiently separated (gap2/D = 8) showed significantly better performance than the other layouts. This could be attributed to the fact that there were many wind turbines that were not in the wake, and the close spacing resulted in a significant increase in the output owing to the induced velocity on each turbine.

In the case of the CO–8 × 2 layout of the 2-D VAWTs with large solidities discussed in this paper, an approximately 20% increase in AEP could be expected compared to the total AEP of isolated turbines by matching the condition that maximizes the average WF output to the prevailing wind direction in a nonuniform wind-direction distribution.

It may be difficult to obtain results similar to those described above for practical-sized lift-type VAWTs because of their low solidity. The optimal layout for lift-type VAWTs should be investigated in the future by combining the proposed method with a GA and other methods.