A Fast Two-Dimensional Numerical Method for the Wake Simulation of a Vertical Axis Wind Turbine

: In the array design of the vertical axis wind turbines (VAWT), the wake effect of the upstream VAWT on the downstream VAWT needs to be considered. In order to simulate the velocity distribution of a VAWT wake rapidly, a new two-dimensional numerical method is proposed, which can make the array design easier and faster. In this new approach, the ﬁnite vortex method and vortex particle method are combined to simulate the generation and evolution of the vortex, respectively, the fast multipole method (FMM) is used to accelerate the calculation. Based on a characteristic of the VAWT wake, that is, the velocity distribution can be ﬁtted into a power-law function, a new correction model is introduced to correct the three-dimensional effect of the VAWT wake. Finally, the simulation results can be approximated to the published experimental results in the ﬁrst-order. As a new numerical method to simulate the complex VAWT wake, this paper proves the feasibility of the method and makes a preliminary validation. This method is not used to simulate the complex three-dimensional turbulent evolution but to simulate the velocity distribution quickly and relatively accurately, which meets the requirement for rapid simulation in the preliminary array design.


Introduction
Since the energy crisis and global warming are getting more serious, the development of renewable energy has received increasing attention. As common renewable energy, wind energy has been well developed due to its wide distribution and nonpollution. A review published by the Intergovernmental Panel on Climate Change (IPCC) shows that the contribution of wind power to the global electricity supply in 2050 will reach 13-14% [1]. There are two common types of wind turbines: horizontal axis wind turbine (HAWT) and vertical axis wind turbine (VAWT). Because of those advantages, such as no need for a yaw control mechanism and low manufacturing costs [2], VAWT has become a research hotspot in recent years. Although the power coefficient of a single VAWT is lower than that of a HAWT, however, some previous research showed that the performance of VAWTs in array configuration would improve significantly. The power coefficient (cp) will increase in a counter-rotating VAWT pairs configuration [3], and the power density of the VAWT farms is potentially more than that of the HAWT farms [4]. To make the VAWT array design more efficient, it is necessary to calculate the velocity field quickly. Therefore, a fast simulation method of the VAWT wake is proposed in this paper so that the wake effect can be measured quickly according to the velocity deficit.
Compared with the HAWT, the study on the wake of a VAWT is still limited. Moreover, in the limited experiments on the wake-field, most of them only focus on the near wake within 3 diameter distances [5][6][7]. In general, the VAWT array is designed based on the The process of the wake simulation can be divided into three steps. In the following examples, the whole simulation process can be represented by the combination of the label mentioned above: Li [27] used the experimental data of static airfoil in the wake research, the evolution of the wake was simulated in the form of the discrete vortex with an empirical decay function (The label of this method is 1T-2V-3VE). The actuator line model of Mendoza [28] used the experimental data modified by the dynamic stall model, and the force coefficient was projected onto the grid to simulate the evolution of the flow field (The label of this method is 1T-2F-3G). There is a hybrid method [32] which was used to calculate the forces on the Euler grid, then the vorticity source term was transformed into the vortex particle in a near-wall transition region, and the far-field was simulated in the Lagrangian framework (The label of this method is 1G-2V-3VA). Based on the compromise between accuracy and computational efficiency, the vortex panel method was applied to calculate the blade force and vortex shedding because of its affordable computational costs. As a type of the vortex panel method, the finite vortex method [19] is more applicable in a wider range of working conditions due to its improved Kutta condition. However, the simulation of a VAWT wake is a challenging task; the wake evolution of the VAWT is much more complex than that of HAWT, such as the breakup and merging behavior of the vortex swarm. Moreover, there is also a blade-vortex interaction (BVI) problem in the VAWT wake; that is, the blade will encounter the shedding vortex at the downstream half revolution. Nevertheless, compared with the grid-based method, the vortex particle method will be more competent for this job because of its low numerical dissipation and easy-to-parallel. In addition, according to the characteristics of the VAWT wake-field, high vorticity is mainly distributed in the upper and lower rows of the wake instead of the whole wake area. In contrast, the whole domain needs to be calculated in the grid-based method, while only a small region where the vorticity is relatively high needs to be calculated in the vortex particle method. In other words, the calculation domain of the vortex particle method is smaller than that of the grid-based method. Thus, it can be concluded that the evolution of the VAWT wake can be calculated quickly and accurately by the combination of the finite vortex method and the vortex particle method. The corresponding label of this method proposed in this paper is 1 V-2 V-3VA.
To save computing costs, the two-dimensional model was often used in the preliminary array design [33][34][35]. Even in the atmospheric boundary layer, the wake can also be simplified as the two-dimensional Gaussian wake model [36]. In addition, the aspect ratio (height/diameter) of the VAWT in real engineering practice is usually big, so that more wind energy can be captured by increasing the height without increasing the land area. Therefore, the influence of the three-dimensional effect, which is caused by the blade end, will become smaller in the case of a large aspect ratio. The above reasons make the two-dimensional numerical method more suitable in the VAWT array design.
When the two-dimensional numerical model is used in the preliminary design of the VAWT array, both the blade force and the velocity field should be modified by the threedimensional effect correction model. There are some three-dimensional effect correction models for the force and performance [37,38], but there are few three-dimensional effect correction models for the VAWT wake. By comparing the three-dimensional and twodimensional wake simulation [39], it can be found that the length of the three-dimensional VAWT wake is shorter than that of the two-dimensional VAWT wake because velocity deficit will recover in the third dimension. In order to consider the three-dimensional effect in the velocity recovery, a new correction model of the VAWT wake is proposed, which is the projection coefficient from the two-dimensional to three-dimensional along the streamwise direction.
It is true that the turbulent evolution of the VAWT wake is a three-dimensional problem. However, the purpose of this paper is to propose a numerical method that can calculate the velocity distribution quickly and make the VAWT array design more efficient; the complex details of the three-dimensional turbulence evolution can be ignored. Therefore, it is not necessary to pay too much attention to the turbulent problems such as the dynamic stall and the turbulent inflow conditions because these turbulence details are not important for the preliminary array design and will be covered by the three-dimensional effect correction model.
In general, according to the tradeoff between accuracy and efficiency, a two-dimensional numerical method with necessary three-dimensional effect correction is proposed; it can calculate the wake-field quickly and relatively accurately, the first-order approximation of the velocity distribution is sufficient, which meets the requirement in the preliminary array design. More details of this method will be introduced in Section 3.

Methodology
The whole process of this new method is like this: First, the finite vortex method is used to calculate the blade load and vortex shedding. Then, the shedding vortex is discretized into the vortex particles to simulate the wake evolution. Finally, it is necessary to introduce a three-dimensional effect correction model into the two-dimensional velocity field so that the modified simulation results can be used in the preliminary array design. Details are given in the following subsections.

Finite Vortex Method
The forces and circulation of the blade can be calculated accurately by the vortex panel method, which was confirmed in the published papers [18,20]. The finite vortex method is one type of the vortex panel method [19,40]. A schematic of the vertical axis turbine in the reference coordinate is shown in Figure 1. A close-up view of the blade is shown in Figure 2. The motion of the blade can be described by the translational velocity U and angular velocity Ω relative to the local coordinate system o. The inflow velocity is V A. There are N sources and sinks distributed on the surface of the blade and a linear vortex lineγon the chord line. The schematic diagram is shown in Figure 3.
Energies 2021, 14, x FOR PEER REVIEW field so that the modified simulation results can be used in the preliminary a Details are given in the following subsections.

Finite Vortex Method
The forces and circulation of the blade can be calculated accurately b panel method, which was confirmed in the published papers [18,20]. The method is one type of the vortex panel method [19,40]. A schematic of the turbine in the reference coordinate is shown in Figure 1. A close-up view of shown in Figure 2. The motion of the blade can be described by the translatio U and angular velocity Ω relative to the local coordinate system o. The inflo VA. There are N sources and sinks distributed on the surface of the blade and a lineγon the chord line. The schematic diagram is shown in Figure 3.   field so that the modified simulation results can be used in the preliminary array Details are given in the following subsections.

Finite Vortex Method
The forces and circulation of the blade can be calculated accurately by th panel method, which was confirmed in the published papers [18,20]. The finit method is one type of the vortex panel method [19,40]. A schematic of the vert turbine in the reference coordinate is shown in Figure 1. A close-up view of the shown in Figure 2. The motion of the blade can be described by the translational U and angular velocity Ω relative to the local coordinate system o. The inflow ve VA. There are N sources and sinks distributed on the surface of the blade and a linea lineγon the chord line. The schematic diagram is shown in Figure 3.       The total velocity potential φ satisfies the Laplace equation in the global coordinate system, as shown in Equation (1) below: The surface of the blade satisfies no through-flow boundary condition, as shown in Equation (2), where r is the vector from the local coordinate system o to the blade surface, n b is the unit normal vector that points to the interior of the blade, the symbol S b denotes the surface of the blade.
The symbol p is the field point; in the far-field, the induced velocity potential ϕ satisfies lim p→∞ ∇ϕ = 0. At the initial time t = 0, the induced velocity potential satisfies the relationship shown in Equation (3): Some experimental studies [41,42] on the flow field showed that there is a finite pressure difference at the trailing edge under the unsteady conditions, which makes the streamline roll-up. In order to solve the stall condition at a high angle of attack, the improved Kutta condition set up a finite pressure difference at the trailing edge. Therefore it is called the finite vortex method [19,43]. The subscripts u and d represent the up and downsides of the trailing edge. The improved Kutta condition of the finite vortex model is shown in Equation (4): In order to determine the value of this pressure difference, set a parameter l to satisfy the relation: w is the circulation of the vortex shedding at the step k. By adjusting the value of l, the lift coefficient of numerical results would be kept consistent with the lift coefficient obtained by the experiment at different angles of attack [19]. The value of l is close to 0 at the small angle of attack, which indicates that the computational error of conventional Kutta condition (pressure difference at trailing edge is equal to 0) is not large at the small angle of attack. With the increase of angle of attack, the absolute value of l is also increasing, which indicates that the pressure difference of the trailing edge is also increasing gradually. Therefore, due to the existence of the correction factor l and the finite value of pressure difference at the trailing edge, the aerodynamic loads calculated by the finite vortex method can still be consistent with the experimental data at a high angle of attack. It also means that the finite vortex method can be more accurate than the conventional vortex panel method in the case of a high angle of attack. According to the unsteady Bernoulli equation (as shown in Equation (5)), the pressure in the flow field is: As the up and downsides of the trailing edge are very close to each other, the vector of the trailing edge can be written as Equation (6): Therefore, after the linearization, the Kutta condition can be rewritten as Equation (7): where mean disturbance velocity vector V TE of the trailing edge is shown in Equation (8): According to Kelvin's theorem, the total circulation does not change with time under the potential flow theory; in other words, the sum of circulation of the bound vortex (Γ f ) and free vortex (Γ w ) does not change with time. Hence, there is Equation (9). At different azimuthal angles, the lift and circulation of the airfoil will change. In each time step, the circulation variation of the bound vortex and that of the free vortex are equal in value, but the sign of them are opposite. It can be expressed as Equation (10). Here, the subscript f represents the circulation of the foil, and the subscript w represents the vortex in the wake, the superscript k represents the time step. Combining Equation (10) with Kelvin's theorem (Equation (9)), Equation (11) is obtained. Then, the circulation of the vortex shedding γ can be calculated by inserting Equation (11) into Equation (7) above. Moreover, the γ obtained here is exactly the circulation that will be discretized into the vortex particles. The detailed process is shown in the references [18][19][20].

Vortex Particle Method
The circulation of the shedding vortex is calculated by the finite vortex method in the potential flow. Then, it will be discretized into a vortex particle swarm to calculate the wake evolution in viscous flow. The vortex particle program is divided into six subsections as follows, and the definition of symbols in each equation is only applicable to the relevant subsection.

Discretization of the Circulation
In this paper, the Vatistas model [44] is used as the vortex core model, which is a common smoothing function to remove the singularity, and it is also similar to the Gaussian smoothing function. The Vatistas model also has been applied in the research of a HAWT wake [13]. The regularization function K is shown in Equation (12): (where ρ = r/r c and r c is the radius of the vortex core): Which satisfies the normalization condition, as shown in Equation (13): The function of vorticity distribution is shown in Equation (14):

Vorticity Transport Equation
The vortex particle method is a viscous method, the governing equation as follows is the vorticity transport equation, which is obtained by taking the curl of the incompressible NS equation, as shown in Equation (15): It is worth mentioning that in the two-dimensional case, there is no stretching term, so the control equation is simplified as Equation (16): It can be seen that the governing equation of the vortex particle method in the twodimensional case is much simpler than the NS equation of the grid-based method.

The Calculation of the Velocity of Vortex Particles
The convection of the vortex particles is calculated by the Biot-Savart Equation.
The induced velocity of the Vatistas model in Equation (17): It can be seen that Equation (17) is very close to the induced velocity equation of a two-dimensional singular vortex, as shown in Equation (18): In this paper, the fast multipole method (FMM) [45] is introduced to accelerate the calculation of the induced velocity of the vortex particles. Here, the induced velocity can be written in the form of induced velocity potential ϕ, z 0 is the target position, z is the position of the vortex particle, and z c is the centroid position of the block. The derivation of the Green function is shown in Equation (19): ξ k k to the second logarithmic term on the right side of Equation (19), Equation (19) can be transformed into the following form (Equation (20)): There are two auxiliary functions (I k (z − z c ) and O k (z 0 − z c )) defined here: The final expression of velocity potential ϕ (Equation (23)) and accuracy criterion of order p (Equation (24)) can be written as follows: (Γ is the circulation of each vortex particle, p is the order of Taylor expansion, N is the number of the vortex particles, the error ε is set to 10 −5 in this paper, R is the diagonal length of the block. Moreover, the black dot here represents the scalar multiplication.) There is an error criterion to judge whether the FMM can replace the direct calculation. The order of the FMM program is NlogN; by contrast, the order of the direct calculation method is N 2 [46,47]. In the FMM program, it is necessary to allocate the whole vortex particles into a tree structure and store them in the leaf nodes of the tree structure. The schematic diagram of the fast multipole method is shown in Figure 4.
The final expression of velocity potential φ (Equation (23)) and accuracy criterion of order p (Equation (24)) can be written as follows: (Γ is the circulation of each vortex particle, p is the order of Taylor expansion, N is the number of the vortex particles, the error ε is set to 10 −5 in this paper, R is the diagonal length of the block. Moreover, the black dot here represents the scalar multiplication.) There is an error criterion to judge whether the FMM can replace the direct calculation. The order of the FMM program is NlogN; by contrast, the order of the direct calculation method is 2 N [46,47]. In the FMM program, it is necessary to allocate the whole vortex particles into a tree structure and store them in the leaf nodes of the tree structure. The schematic diagram of the fast multipole method is shown in Figure 4.

The Calculation of the Vorticity of Vortex Particles
When the position of the vortex particles is updated, it is necessary to update the vorticity of the vortex particles. There are two common methods to simulate the viscous dissipation. One is called the random walk algorithm; the other is the particle strength exchange method (PSE). In this paper, the latter method (PSE) is used. In order to calculate the vorticity dissipation term, a transformation Equation (25) is introduced here to transform the second derivative form of the NS equation into the integral or summation form [48] (The symbol f in the following equation represents the function of vorticity ω).

The Calculation of the Vorticity of Vortex Particles
When the position of the vortex particles is updated, it is necessary to update the vorticity of the vortex particles. There are two common methods to simulate the viscous dissipation. One is called the random walk algorithm; the other is the particle strength exchange method (PSE). In this paper, the latter method (PSE) is used. In order to calculate the vorticity dissipation term, a transformation Equation (25) is introduced here to transform the second derivative form of the NS equation into the integral or summation form [48] (The symbol f in the following equation represents the function of vorticity ω). where: The governing equation is transformed into the summation form (Equation (27)): Subscript p represents the target point, q is the vortex particles, and the symbol σ is the radius of the vortex particle.

The Calculation of the Turbulent Viscosity Coefficient
It is worth mentioning that the viscosity coefficient ν in Equation (27) includes physical viscosity and turbulent viscosity, and the turbulent viscosity is often larger than the physical viscosity. As mentioned in Section 2, in the two-dimensional simulation, no matter what numerical method is used, the three-dimensional effect correction model is needed. However, the turbulence model is still necessary, so here is a brief introduction. The large eddy simulation (LES) in the Lagrangian frame [13] is introduced to calculate the turbulent viscosity in the flow field. The method proposed by Mansfield et al. [49] in Lagrangian large eddy simulation can establish the relationship between filter scale and vortex particle scale. The smoothing function K is related to the particle filter G; the relationship is shown in Equation (28): where c is a constant, the filtering scale ∆ and the radios of the vortex particle σ satisfy the relationship: ∆ = cσ. The particle filter and circular filter meet the energy conservation equation in the wavenumber domain. By numerical calculation, the solution c = ∆/σ = 2.83 can be obtained. More numerical processes can be found in the work of Liu [13]. When the filtering scale ∆ is obtained, the turbulent viscosity can be calculated from the following Equation (29) according to the LES Smagorinsky model [50].
3.2.6. The Redistribution of the Vortex Particles In order to ensure the convergence of the vortex particle method, a redistribution algorithm is in need [51]. Especially when the distance between particles is too large, the redistribution of particles can ensure a certain overlap rate between particles. There are some interpolation stencils, which distribute the vortex particle on the surrounding nodes with different weights. In this paper, the third-order interpolation function (Equation (30)) is used where λ = x/h, x is the distance from the center of mesh to the old particles, h is the mesh size. Each old vortex particle will be replaced by nine new particles on the surrounding mesh nodes, as shown in Figure 5.

Discrete Parameters and Relevant Convergence Criteria
In this paper, the PIV experiment [10] is used as the benchmark test case for validation. The experimental parameters are shown in the Table 2 below. The relevant experience of the discrete parameters and the corresponding convergence criteria will be introduced one by one: The finite vortex method is not a time-consuming algorithm, and in order to ensure convergence accuracy, in this paper, there are 200 discrete elements on the blade surface. The time step of the finite vortex method is the time for 1° rotation, and it is a common discrete scale, which is proved by another paper [52]. The numerical parameters in the finite vortex method are conservative, and the accuracy of the program also has been verified by the previously published paper [18].
The diameter of the vortex core in the vortex particle method should be slightly larger than the chord length of VAT. In this case, the diameter of the vortex core is 1.2 times the chord length. This experience is also similar to the selection of the smooth length in the actuator line model [53].
It is worth mentioning that the vortex particle method is a mesh-free method. Therefore, there is no need to verify the independence of mesh size, and it is diffusion free in theory [54]. As long as the number of the vortex particles is enough, the wake evolution could be accurately simulated. The redistribution program in the vortex particle method redistributes each vortex particle to nine new particles through the interpolation function, which satisfies the conservation of the circulation from the zero-order moment to the second-order moment. In this way, the number of the vortex particles will be enough at the place where the vorticity is high, and the vortex particle methods will be automatically adaptive and will concentrate on the regions of physical interest [55]. Since the number of the vortex particles is adjusted adaptively by the redistribution program, the initial number of the vortex particles has little impact on the simulation of the wake-field. Therefore, it is not necessary to verify the independence of the number of vortex particles. In general, the numerical viscosity of the grid-based method will produce a big diffusive error in high

Discrete Parameters and Relevant Convergence Criteria
In this paper, the PIV experiment [10] is used as the benchmark test case for validation. The experimental parameters are shown in the Table 2 below. The relevant experience of the discrete parameters and the corresponding convergence criteria will be introduced one by one: The finite vortex method is not a time-consuming algorithm, and in order to ensure convergence accuracy, in this paper, there are 200 discrete elements on the blade surface. The time step of the finite vortex method is the time for 1 • rotation, and it is a common discrete scale, which is proved by another paper [52]. The numerical parameters in the finite vortex method are conservative, and the accuracy of the program also has been verified by the previously published paper [18].
The diameter of the vortex core in the vortex particle method should be slightly larger than the chord length of VAT. In this case, the diameter of the vortex core is 1.2 times the chord length. This experience is also similar to the selection of the smooth length in the actuator line model [53].
It is worth mentioning that the vortex particle method is a mesh-free method. Therefore, there is no need to verify the independence of mesh size, and it is diffusion free in theory [54]. As long as the number of the vortex particles is enough, the wake evolution could be accurately simulated. The redistribution program in the vortex particle method redistributes each vortex particle to nine new particles through the interpolation function, which satisfies the conservation of the circulation from the zero-order moment to the second-order moment. In this way, the number of the vortex particles will be enough at the place where the vorticity is high, and the vortex particle methods will be automatically adaptive and will concentrate on the regions of physical interest [55]. Since the number of the vortex particles is adjusted adaptively by the redistribution program, the initial number of the vortex particles has little impact on the simulation of the wake-field. Therefore, it is not necessary to verify the independence of the number of vortex particles. In general, the numerical viscosity of the grid-based method will produce a big diffusive error in high Reynolds number. In contrast, the redistribution program does not bring discernible numerical dissipation.
Redistribute the vortex particles in 2D simulations, which shows substantially lower errors than infrequent redistribution [46]. The redistribution program could provide a potential solution for the accurate simulation of the vorticity [56]. It can be seen that the redistribution is of great significance to the convergence of the vortex particle method.
The mesh size of the redistribution program is h = 0.0025 m. In order to ensure that there is a certain overlap rate between adjacent particles, the convergence criterion required that the vortex particle radius should be larger than the mesh size [46]. Hence, the radius of the vortex particle is r 0 = 0.0033 m.
According to the equation πD/c ≈17.4, the circumference is 17.4 times of the chord length. Therefore, it is set to discharge the vortex 18 times in one rotation period. That is 360 • /18 = 20 • ; the time interval of the vortex shedding is the time for 20 • rotation. The corresponding time step here is 0.0028 s. According to the stability criterion [46], the time step should be less than the reciprocal of the maximum vorticity. In this view, the maximum allowable vorticity under the convergence criterion should be 1/0.0028 s ≈350 s −1 . Although the value of vorticity changes with time due to the redistribution, the magnitude of vorticity is within ten in this case, which obviously satisfies the stability criterion.
It can be seen from the benchmark experiment that the blocking effect can be ignored. Therefore, the wake evolution in this paper is simulated under unbounded conditions, and the blocking effect has not been involved. In the vortex particle method, there is no need to set special boundary conditions for unbounded flow. This is a feature different from the grid method.

The Three-Dimensional Effect Correction Model of Wake
This new correction model is a projection function along the streamwise direction to correct the three-dimensional effect of the velocity field. As we all know, the threedimensional effect of the wake is that the third dimension will accelerate the recovery of the velocity deficit in the two-dimensional wake. By multiplying the two-dimensional velocity deficit field with a projection function, the influence of the third dimension on the velocity recovery can be equivalent.
The correction process is based on some empirical fitting functions of the VAWT wake. According to the analysis of PIV experiment results [57], the velocity deficit of the wake can be fitted into the power-law function of the downstream position X: where X t is the downstream transition location which is related to the dynamic solidity σ D , and it can be obtained by this equation: X transition = D(4.78 − 4.93σ D ), where D is the diameter of the VAWT. These three coefficients, c 1 c 2 c 3 , are the empirical fitting coefficients which are influenced by solidity and tip speed ratio. Moreover, this form of the power-law function is also applicable to the HAWT wake and the two-dimensional or three-dimensional bluff body wake [57]. Therefore, the process of the three-dimensional effect correction proposed in this paper is based on the three-dimensional velocity distribution functions summarized by the previous study [57]. Then, we divided these three-dimensional velocity distribution functions by a two-dimensional velocity distribution function, and the proportion obtained by the division is the projection coefficient from the two-dimensional to the three-dimensional, which is exactly what we want to correct the three-dimensional effect of the wake-field. Finally, the two-dimensional velocity field is multiplied by the projection coefficient along the streamwise direction, and the two-dimensional simulation results considering the three-dimensional effect correction can be used in the preliminary array design.
In order to get the appropriate three-dimensional velocity distribution function which is close to the benchmark case in this paper, two models are selected from the previous study [57] because the parameters of these two models are relatively close to that of the benchmark case, and they were marked as VAWT_1 and VAWT_2, respectively. The velocity distribution of the three-dimensional wake models (VAWT_1 and VAWT_2) and the general form of the velocity distribution of the two-dimensional wake model are described in Table 3 below. Table 3. The parameters of these three vertical axis wind turbine (VAWT) wake models.

Solidity The Downstream
Transition Location

The Power Function of Wake Velocity Deficit
Benchmark case (two-dimensional wake model) It can be seen from these functions that the power-law exponent of the two-dimensional wake is −0.5, which is significantly larger than that of the other two three-dimensional wake models (−0.684 and −0.622). That is to say, due to the three-dimensional effect, the decay rate of the two-dimensional wake model is slower than that of the three-dimensional model.
In addition, it should be noted that the two models (VAWT_1 and VAWT_2) selected here are relatively close to the model in this paper; this can be seen in Figure 25 in the reference [57] that there is little difference in the power function between different VAT conditions. Therefore, there is no need to be too concerned about whether the selection of the experimental model is completely consistent with the validation model used in this paper.
In Figure 6, the curves formed by the square symbol ( ) represent the power functions of the three-dimensional wake models (VAWT_1 and VAWT_2), the curve formed by the closed circle symbol (•) represents the general form of the power function of the twodimensional wake model. The functions of these three curves are described above. Then, these two power functions of the three-dimensional wake models (symbol ) are divided by the power function of the two-dimensional wake model (symbol •), and the ratio functions obtained (symbol ×) are the correction coefficient along the streamwise, which could be used to correct the velocity deficit from the two-dimensional to the three-dimensional. There are two other considerations: The first is to eliminate the unreasonable situation that the ratio is much greater than one within one diameter. Second, considering that the three-dimensional effect will make the power coefficient of VAT smaller than that of the two-dimensional model [33], and there is a positive correlation between the power coefficient and the velocity deficit, that is to say, the velocity deficit and power coefficient of the three-dimensional model are smaller than those of the two-dimensional model. Hence, the correction coefficient should be slightly less than 1 in the near field. Based on these two considerations, combine the two ratio functions (symbol ×), what we get is the empirical coefficient (symbol♦) related to the streamwise position, and the combined function (symbol♦) can be approximately regarded as the final three-dimensional effect Figure 6. The correction coefficient along the streamwise direction. The X-axis is the downstream position normalized by the rotor diameter (X/D). The Y-axis represents two normalized parameters; one is the normalized velocity deficit ∆U/U 0 (symbol • and ) of the wake model, the other is the correction coefficient (symbol × and ) obtained by the division.
There are two other considerations: The first is to eliminate the unreasonable situation that the ratio is much greater than one within one diameter. Second, considering that the three-dimensional effect will make the power coefficient of VAT smaller than that of the two-dimensional model [33], and there is a positive correlation between the power coefficient and the velocity deficit, that is to say, the velocity deficit and power coefficient of the three-dimensional model are smaller than those of the two-dimensional model.
Hence, the correction coefficient should be slightly less than 1 in the near field. Based on these two considerations, combine the two ratio functions (symbol ×), what we get is the empirical coefficient (symbol ) related to the streamwise position, and the combined function (symbol ) can be approximately regarded as the final three-dimensional effect projection function which is used in this paper.
It can be found that the curve is a function of the distance along the streamwise direction; by multiplying the two-dimensional velocity field with the projection coefficient (symbol ) along the streamwise direction, the wake-field corrected by the threedimensional effect can be obtained. Moreover, this is the whole process of the threedimensional effect correction for two-dimensional wake simulation in this paper.
In this paper, the power-law function is proposed to correct the three-dimensional effect, which will be a universal idea. Not only the wake of HAWT and VAWT but also the wake of the two-dimensional or three-dimensional bluff body conform to the form of the power-law function [57,58]. Moreover, the two-dimensional Gaussian power function was also used as the basic wake model for array layout design [35,59]. Even if the influence of the boundary layer on the wake is considered, the velocity distribution function could also be fitted into the form of the power-law function [60].
In fact, the correction model is widely used in real engineering practice to simplify the simulation, and it can be equivalent to the effect of some complicated turbulence details, such as dynamic stall and complex turbulent inflow conditions. Especially for the array design, the first-order approximation of the velocity distribution is enough to meet the demand, and the three-dimensional effect correction of the two-dimensional wake simulation is a good choice. The correction model proposed in this paper is derived from the power-law function of velocity distribution, which could make the numerical results in good agreement with the three-dimensional results. The validation of this numerical method is shown in Section 4.

Validation and Discussion
In order to verify the reliability of the new numerical method proposed in this paper, a PIV experiment is selected to compare with the simulation results. The parameters of this PIV experiment have been described in detail in Table 2 in Section 3.3.
The comparison results are shown in Figure 7. Relatively speaking, the average velocity profiles calculated by the proposed method match well with the PIV experimental data [10]. At least for the preliminary array design, it is an encouraging result to achieve at least first-order accuracy in the two-dimensional wake simulation. It also proves that the two-dimensional numerical method with a reasonable three-dimensional correction is an acceptable and practicable approach that can make the wake simulation faster and the array design easier. In Figure 8, the rotation direction of the VAWT is counterclockwise. According to the analysis of velocity distribution of a VAWT wake, the peak value of velocity deficit is found not on the centerline. This asymmetry can be seen in the velocity field ( Figure 8).
It needs to be clarified here, in the preliminary array design, only the time-averaged velocity field needs to be considered, and the unsteady flow details can be ignored. In addition, since the validation data of the reference is time-averaged, and the benchmark test case did not provide the transient experimental data at different times. In order to keep consistent with the benchmark test case, the time-averaged velocity field was used in this paper.
In the vorticity field (Figure 9), due to the velocity deficit caused by the turbine and the conservation of momentum, the wake vortex will expand in the near field, and the expansion phenomenon of the stream tube in the numerical simulation is similar to the flow field observed by another PIV experiment [5]. It can also be seen that there is a skewness in the distribution of the wake vortex. According to the Magnus effect, a lateral force perpendicular to the inflow direction will be produced on a spinning object. In addition, according to Newton's third law, the forces acting on the turbine and the forces acting on the flow field belong to action and reaction. In this case, the VAWT rotates counterclockwise and is subjected to a downward lateral force. Therefore, the wake is subjected to the upward force and will move upward. The similar phenomena of the wake skewness have also been found in some other VAT wake experiments [8,11]. It can be seen that there is a skewness in the VAT wake caused by the Magnus effect. Moreover, the faster the rotation, the greater the wake skewness. The explanation that wake skewness is related to the Magnus effect can also be found in other studies [5,9,61]. However, from the micro point of view, the generation of the wake skewness is caused by the wind-blade interaction, and it can also be explained from the macro point of view that the wake skewness is caused by the Magnus effect. Both explanations are correct. By comparison, there is usually no skewness in the HAWT wake, this consideration will make the array layout of VAT different from that of HAT. Moreover, more suggestions can be provided for the array arrangement of VATs, such as, the arrangement distance for the downstream VAT should be different considering the rotation direction of the upstream VAT.

Validation and Discussion
In order to verify the reliability of the new numerical method proposed in this paper, a PIV experiment is selected to compare with the simulation results. The parameters of this PIV experiment have been described in detail in Table 2 in Section 3.3.
The comparison results are shown in Figure 7. Relatively speaking, the average velocity profiles calculated by the proposed method match well with the PIV experimental data [10]. At least for the preliminary array design, it is an encouraging result to achieve at least first-order accuracy in the two-dimensional wake simulation. It also proves that the two-dimensional numerical method with a reasonable three-dimensional correction is an acceptable and practicable approach that can make the wake simulation faster and the array design easier. In Figure 8, the rotation direction of the VAWT is counterclockwise. According to the analysis of velocity distribution of a VAWT wake, the peak value of velocity deficit is found not on the centerline. This asymmetry can be seen in the velocity field ( Figure 8).
It needs to be clarified here, in the preliminary array design, only the time-averaged velocity field needs to be considered, and the unsteady flow details can be ignored. In addition, since the validation data of the reference is time-averaged, and the benchmark test case did not provide the transient experimental data at different times. In order to keep consistent with the benchmark test case, the time-averaged velocity field was used in this paper.
. Figure 7. The comparison of average velocity profiles at six different cross-sections. The X-axis represents the sum of the two normalized parameters; one is the streamwise position, which is normalized by the turbine diameter (X/D), the other is the velocity deficit, which is normalized by the inflow velocity (ΔU/U0). The Y-axis is the lateral position, which is normalized by the turbine diameter (Y/D). The data of black curves are from a published paper [10]; the color curves are calculated by the numerical method proposed in this paper. Figure 8. The time-average velocity field of the VAWT wake. The X-axis is the streamwise position, which is normalized by the turbine diameter (X/D). The Y-axis is the lateral position, which is normalized by the turbine diameter (Y/D). The color bar is the velocity deficit, which is normalized by the inflow velocity (ΔU/U0). The blue dashed lines indicate that the velocity profiles at that location are used for comparison with the PIV experimental results. Figure 7. The comparison of average velocity profiles at six different cross-sections. The X-axis represents the sum of the two normalized parameters; one is the streamwise position, which is normalized by the turbine diameter (X/D), the other is the velocity deficit, which is normalized by the inflow velocity (∆U/U 0 ). The Y-axis is the lateral position, which is normalized by the turbine diameter (Y/D). The data of black curves are from a published paper [10]; the color curves are calculated by the numerical method proposed in this paper.

Validation and Discussion
In order to verify the reliability of the new numerical method proposed in this paper, a PIV experiment is selected to compare with the simulation results. The parameters of this PIV experiment have been described in detail in Table 2 in Section 3.3.
The comparison results are shown in Figure 7. Relatively speaking, the average velocity profiles calculated by the proposed method match well with the PIV experimental data [10]. At least for the preliminary array design, it is an encouraging result to achieve at least first-order accuracy in the two-dimensional wake simulation. It also proves that the two-dimensional numerical method with a reasonable three-dimensional correction is an acceptable and practicable approach that can make the wake simulation faster and the array design easier. In Figure 8, the rotation direction of the VAWT is counterclockwise. According to the analysis of velocity distribution of a VAWT wake, the peak value of velocity deficit is found not on the centerline. This asymmetry can be seen in the velocity field ( Figure 8).
It needs to be clarified here, in the preliminary array design, only the time-averaged velocity field needs to be considered, and the unsteady flow details can be ignored. In addition, since the validation data of the reference is time-averaged, and the benchmark test case did not provide the transient experimental data at different times. In order to keep consistent with the benchmark test case, the time-averaged velocity field was used in this paper.
. Figure 7. The comparison of average velocity profiles at six different cross-sections. The X-axis represents the sum of the two normalized parameters; one is the streamwise position, which is normalized by the turbine diameter (X/D), the other is the velocity deficit, which is normalized by the inflow velocity (ΔU/U0). The Y-axis is the lateral position, which is normalized by the turbine diameter (Y/D). The data of black curves are from a published paper [10]; the color curves are calculated by the numerical method proposed in this paper. Figure 8. The time-average velocity field of the VAWT wake. The X-axis is the streamwise position, which is normalized by the turbine diameter (X/D). The Y-axis is the lateral position, which is normalized by the turbine diameter (Y/D). The color bar is the velocity deficit, which is normalized by the inflow velocity (ΔU/U0). The blue dashed lines indicate that the velocity profiles at that location are used for comparison with the PIV experimental results. Figure 8. The time-average velocity field of the VAWT wake. The X-axis is the streamwise position, which is normalized by the turbine diameter (X/D). The Y-axis is the lateral position, which is normalized by the turbine diameter (Y/D). The color bar is the velocity deficit, which is normalized by the inflow velocity (∆U/U 0 ). The blue dashed lines indicate that the velocity profiles at that location are used for comparison with the PIV experimental results. can also be explained from the macro point of view that the wake skewness is caused by the Magnus effect. Both explanations are correct. By comparison, there is usually no skewness in the HAWT wake, this consideration will make the array layout of VAT different from that of HAT. Moreover, more suggestions can be provided for the array arrangement of VATs, such as, the arrangement distance for the downstream VAT should be different considering the rotation direction of the upstream VAT. Figure 9. The vorticity field at the near wake. The X-axis is the streamwise position which is normalized by the turbine diameter (X/D). The Y-axis is the lateral position which is normalized by the turbine diameter (Y/D). The color bar is the nondimensional vorticity which is normalized by the maximum vorticity in the field (ω/ωmax).
There is another interesting phenomenon in the far wake in Figure 10. It can be seen from the following eight moments in one period (the period here is not the rotation period, but the period of the vortex street oscillation), due to the complex evolution, the vortex pairs in the far-field are similar to the Karman vortex street in the flow around a cylinder. This characteristic can also be found in another reference [62]. The similar wake evolution will occur in the far-field of a HAT wake as well [39]. It is still unknown whether this complex evolution is similar to the meandering of the HAWT wake. However, at least, this interesting phenomenon in these eight images shown here can lead to a new research direction in the future. There is another interesting phenomenon in the far wake in Figure 10. It can be seen from the following eight moments in one period (the period here is not the rotation period, but the period of the vortex street oscillation), due to the complex evolution, the vortex pairs in the far-field are similar to the Karman vortex street in the flow around a cylinder. This characteristic can also be found in another reference [62]. The similar wake evolution will occur in the far-field of a HAT wake as well [39]. It is still unknown whether this complex evolution is similar to the meandering of the HAWT wake. However, at least, this interesting phenomenon in these eight images shown here can lead to a new research direction in the future.
It is also worth mentioning that the vorticity of the vortex particle swarm is different at each time step, and the redistribution program makes the number and vorticity of the vortex particles change all the time. Therefore, the peak value of the vorticity field (ω max ) and the range of the color bar (ω) also change at each time step. In order to unify the color bar for observation, the vorticity field (Figures 9 and 10) is normalized by the maximum vorticity at that moment. The color bars of the wake-field are represented by the dimensionless vorticity (ω/ ω max ). In general, although the redistribution program makes the range of vorticity field change all the time, for the vortex particle method, the whole field still satisfies the conservation condition, so it may not be necessary to pay too much attention to the absolute value of the vorticity at each moment.
However, the purpose of this paper is to propose a fast simulation method for the VAWT wake, so the complex wake evolution phenomenon will not be discussed in depth. The simulation results are calculated by the code written by the authors, and there is no commercial software to use this numerical method so far. In terms of the calculation efficiency of the algorithm, the vortex particle method takes less than 2 h on 4-core CPUs, and the number of the vortex particles is about a hundred thousand, when the focus is only on the flow field within 10D, as we all know, the calculation cost here is much less than that of the grid-based method.  Figure 10. The eight moments of vorticity field in one period at a wide view. The interval of each moment is approximately equal to one rotation period. The X-axis is the streamwise position, which is normalized by the turbine diameter (X/D). The Y-axis is the lateral position, which is normalized by the turbine diameter (Y/D). The color bar is the nondimensional vorticity, which is normalized by the maximum vorticity in the field (ω/ωmax).
It is also worth mentioning that the vorticity of the vortex particle swarm is different at each time step, and the redistribution program makes the number and vorticity of the vortex particles change all the time. Therefore, the peak value of the vorticity field (ωmax) and the range of the color bar (ω) also change at each time step. In order to unify the color bar for observation, the vorticity field (Figures 9 and 10) is normalized by the maximum vorticity at that moment. The color bars of the wake-field are represented by the dimensionless vorticity (ω/ωmax). In general, although the redistribution program makes the range of vorticity field change all the time, for the vortex particle method, the whole field still satisfies the conservation condition, so it may not be necessary to pay too much attention to the absolute value of the vorticity at each moment.
However, the purpose of this paper is to propose a fast simulation method for the VAWT wake, so the complex wake evolution phenomenon will not be discussed in depth. The simulation results are calculated by the code written by the authors, and there is no commercial software to use this numerical method so far. In terms of the calculation efficiency of the algorithm, the vortex particle method takes less than 2 h on 4-core CPUs, and the number of the vortex particles is about a hundred thousand, when the focus is only on the flow field within 10D, as we all know, the calculation cost here is much less than that of the grid-based method. Figure 10. The eight moments of vorticity field in one period at a wide view. The interval of each moment is approximately equal to one rotation period. The X-axis is the streamwise position, which is normalized by the turbine diameter (X/D). The Y-axis is the lateral position, which is normalized by the turbine diameter (Y/D). The color bar is the nondimensional vorticity, which is normalized by the maximum vorticity in the field (ω/ ω max ).

Conclusions
In this paper, by combining the finite vortex method and the vortex particle method, a fast and relatively accurate two-dimensional numerical method is proposed to simulate the wake-field of the vertical axis wind turbine. Based on the power-law function of velocity distribution, a new correction model is proposed to correct the three-dimensional effect of the VAWT wake. Then, the feasibility of the numerical method is preliminarily verified by the PIV experiment. Finally, according to the flow field analysis, the skewness of the VAWT wake is discussed. In the far-field, the VAWT wake evolves into the Karman vortex street, which is similar to the flow around a cylinder.
Although the published experiments on the VAT wake, which can be used for the verification, are relatively limited, however, the numerical method proposed in this paper is feasible in theory and has been preliminarily verified. It has been proven that this new method can simulate the velocity distribution quickly and relatively accurately. The first-order approximation can be achieved by this two-dimensional numerical method, which meets the requirement for efficiency in the preliminary array design. In the future, with more experimental data published, the numerical method proposed in this paper will be systematically studied and improved. Nevertheless, we believe that it can be a promising approach to make the VAWT array design easier and more reasonable.

Conflicts of Interest:
The authors declare no conflict of interest.