The study on performance and aerodynamics of micro counter-rotating HAWT

Energy 161 (2018) 939e954

Contents lists available at ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

The study on performance and aerodynamics of micro counter-rotating HAWT Li Zhiqiang, Wu Yunke*, Hong Jie, Zhang Zhihong, Chen Wenqi National Key Laboratory of Science and Technology on Aero-Engine Aero-thermodynamics, School of Energy and Power Engineering, Beihang University, Beijing 100191, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 6 February 2018 Received in revised form 15 June 2018 Accepted 9 July 2018 Available online 11 July 2018

Compared with large HAWT, the power coefficient and economical performance of the micro HAWT are inferior. However, due to flexible deployment, rapid installation, easy-to-use and low cost, it is still necessary to carry out research on micro HAWT to improve their performance. As a structure proven to be able to increase the power coefficient of HAWT effectively, the aerodynamic characteristics of the micro counter-rotating HAWT (CRWT) is studied in detail by experiments and LBM-LES simulation in this paper. Firstly, the results show the micro CRWT can enhance the power coefficient by two schemes. One is that the equal diameter counter rotating rotors with certain pitch angle can achieve a maximum power coefficient in a wider range of speed. The other is to use the vice rotor to enhance the power output of the main rotor, whether the vice rotor is in front rotor or rear, the rotation of the vice rotor can significantly improve power coefficient of the main rotor, so as to enhance the power coefficient of the total system. Then the aerodynamic mechanism of the Vice-Main and Main-Vice rotor schemes is further discussed with LBM-LES result in detail. © 2018 Elsevier Ltd. All rights reserved.

Keywords: Wind turbine Counter-rotating HAWT Aerodynamics LBM-LES Wind tunnel experiments

1. Introduction The proportion of global renewable energy consumption, especially the wind energy, is increasing annually. More and more centralized wind farms have been built over the past decade in the emerging economies. The wind power technology is also developing rapidly. At present, the design power of single wind turbine has reached up to 10 MW [1]. Micro horizontal axis wind turbine (HAWT), compared with large HAWT, have lower power coefficient and worse economical performance, so the research work for micro HAWT is relatively limited. But what cannot be ignored is that the micro HAWT are more suitable for the distributed energy supply areas like roofs, villages and roadsides than large HAWT, because of the advantages of flexible deployment, rapid installation, easy-to-use and low cost. So it is still necessary to carry out research on micro wind turbines to improve their performance [2,3]. The CRWT is a structure proven to increase the power coefficient of HAWT effectively. Newman [4,5] derived the onedimensional mathematical model of the double stage rotors and

* Corresponding author. E-mail address: [email protected] (W. Yunke). https://doi.org/10.1016/j.energy.2018.07.049 0360-5442/© 2018 Elsevier Ltd. All rights reserved.

infinite stage rotors of wind turbine. By numerical simulation he thought the theoretical power coefficient limit value was predicted as 64% for 2-stage HAWT, and 66.7% for the infinite-stage rotors HAWT, higher than the famous Betz limit of single rotor HAWT known as 59%. Based on wind tunnel experimental results, Ushiyama et al. [6] demonstrated the technical feasibility of CRWT structure and its potential to improve the HAWT performance. Appa, Debleser, Stephens and Li et al. [7e10] also proposed several structures of small or micro CRWT. These structures all adopted similar structure, which is the two stage rotors counter rotating scheme to improve the power coefficient of HAWT. Li et al. [11e13] studied the effect of design parameters of straight blade pitch angle and rotors’ distance on the performance of the small CRWT (Fig. 1a) by using wind tunnel experiment and numerical simulation method. Gregg et al. [14] tested the airfoil performance of the optimized CRWT with similar methodology (Fig. 1b). And Habash et al. [15] tried to observe the wake flow between the two counter rotating rotors (Fig. 1c). All these studies show that counter-rotating scheme will significantly enhance the power coefficient of wind power system. In order to further improve the practicability of CRWT, Mitulet et al. [16] designed a new type of synchronous electric generator with counter-rotating armatures (Fig. 2a), which can obtain the

940

L. Zhiqiang et al. / Energy 161 (2018) 939e954

(a) Model of Li et al [11~13]

(b) Model of Gregg et al [14]

(c) Model of Habash et al [15]

Fig. 1. Different models of CRWT.

(a) Synchronous electric generator

(b) The power characteris Mitulet et al

of CRWT

Fig. 2. Research Results of Mitulet et al. [16].

two-axle output by only one generator. And they obtained the complete performance curves (Fig. 2b) of the 1 KW level CRWT by wind tunnel test. Through developing the Newman’s one-dimensional mathematical model [4,5], Chan et al. [17] proposed a new design of the upstream rotor with a ‘bladeless’ center in CRWT system (Fig. 3a), and the theoretical power coefficient limit predicted by this model is claimed up to 81.4%. Based on this model, Jung et al. [18] designed a 30 KW level prototype with bevel and planetary gear train (Fig. 3b). The experimental results of the prototype showed that the power coefficient of this preliminary design is only up to 50%. Subsequently, No et al. [19] established a 1 MW demonstration device of this design for commercial purpose of CRWT system (Fig. 3c). In addition, compared with the single rotor HAWT, the VAWT and the ducted HAWT, the test results indicated that the CRWT design can gain the better sensitive yaw, smaller tower loads and higher rotation stabilization [19]. In summary, the ‘center-bladeless’ design of upstream rotor could be helpful for improving the performance of large CRWT. But for micro CRWT, both rotors are much smaller than those of the large CRWT, so the ‘center-bladeless’ scheme could not be applicable. And limited by the larger load increment and manufacturing cost of two rotary rotors, the counter rotating scheme is rather more suitable for the micro HAWT application. Therefore, it is still necessary to continue the study of micro CRWT. By corrected blade element and momentum theory (BEMT) using the vortex lattice method, Hwang et al. [20] established an optimization procedure for a micro CRWT and validated it by a

comparison with the power curve of CRWT test. Then, they [21] studied the hydrodynamic performance of a number of micro CRWT models with different combinations of pitch angles of upstream and downstream blades by the method of water tunnel experiment. It can be seen that most of the previous research works were the performance test. But the fluid dynamics mechanism of CRWT improving the power coefficient has not been clarified. For the research of the flow details near the counter-rotating rotors region, especially the fluid-structure interaction on the surface of upstream and downstream rotors, and the flow interaction mechanism between two rotors, the experimental observation work is relatively difficult. Hence, it generally relies on numerical simulation method. However, because of the counter-rotation of both rotors, numerical simulation, especially modeling and programming using common CFD methods like conventional finite volume method (FVM) or finite difference method (FDM) is considerable difficult. And the numerical procedure modules for rotating parts like compressor and turbine provided by general commercial CFD programs are still not sufficient enough when dealing with two forced motion rotors counter-rotating at unknown rotational speed. More advanced CFD methods, such as dynamic mesh technology, while applied to the counter-rotating rotors under passive rotation condition, have also difficulties in dealing with the mesh interface between two rotors, so it is hard to accurately solve the 6DOF fluid-flow coupling problem in rotational domain [22e28]. To understand the CRWT better, the performance, the interaction between upstream and downstream rotors and the wake flow from the blades of two rotors in CRWT model proposed by Li et al.

L. Zhiqiang et al. / Energy 161 (2018) 939e954

941

[17]

(c) Prototype and commercial type designed by Jung et al[18] and No et al[19]

(b) 30KW level prototype designed by Jung et al[18]

Fig. 3. Latest development of CRWT.

in Ref. [11e13] are further studied using wind tunnel experiment and lattice Boltzmann method (LBM) in this paper. 2. CRWT model As shown in Fig. 4, double rotors coaxial counter-rotating scheme of CRWT is adopted in this study, and each torque from rotor-axles of the upstream and downstream is outputted respectively through the pulleys. The torque & tachometer is connected with the other end of every pulley. To simplify the expression in this paper, the upstream and downstream rotors are called as the front and rear rotors. Lanzafame et al. [29] deemed that the performance gap of the micro HAWT with twisted and non-twisted blades was not significant, but more researches [30e35] for aerodynamic design of micro HAWT indicated that the twisted blades designed with the classical BEM theory could improve the wind power coefficient

efficaciously. Despite the existence of conflicts in micro HAWT design methods, the study presented in this paper is mainly focused on the basic performance and aerodynamic influence between two stage rotors, instead of exclusively discussing the optimal design of blade modeling. The equal chord length and nontwisted straight blades can be convenient for analyzing the separation and stall mechanism in different stream surfaces of blade, so more abundant rules can be obtained by analyzing some experiments and simulation results. Moreover, the BEM optimization design theory for the second stage rotor is still lack of systematic research before this study. Overall consideration, in this paper, six equal chord length and non-twisted straight blades are adopted in each rotor. In the previous study, Li et al. [11e13] proved the rotors with the blade pitch angles of 60
and 30
have great difference in both speed and torque, which is more favorable to study the effect of different counter-rotating schemes on the performance of CRWT. Therefore, in this study, four different research schemes with the blade pitch angles of b ¼ 30
and b ¼ 30
are adopted, as shown in Fig. 5. The NACA4418 airfoil which is widely used in small wind turbines is selected for all of the rotors, and the chord length of each blade is equal from tip to hub with the larger rotors of 50 mm and the smaller rotors of 40 mm. To find the effect of front and rear rotor sizes on the performance of the CRWT, three kinds of schemes with different ratios of the diameters of the front and rear rotors are set up respectively in Fig. 6, where d is defined as the ratio of the diameters of the front and rear rotors:

.

d ¼ Dfront Drear

Fig. 4. Micro CRWT performance measurement setup.

(1)

According to the classification, the research cases are shown in Table 1. Cases 13e16 are the control cases of single rotor with pitch angle of b ¼ 30º and b ¼ 60º respectively. The diameters of rotors in

942

L. Zhiqiang et al. / Energy 161 (2018) 939e954

(a)

front

= 30Ϩ

rear

(b)

= 30Ϩ

front

= 30Ϩ

rear

= 60Ϩ

(c)

front

= 60Ϩ

rear

= 30Ϩ

(d)

front

= 60Ϩ

rear

= 60Ϩ

Fig. 5. Blade Cascade Diagrams of Different Pitch angles.

speed in test section with accuracy higher than 97.5%. The CRWT model shown in Fig. 1 a, the torque & tachometer shown in Fig. 4 and the DC generator are integrated as a whole then placed into the main test section of the wind tunnel. The output signals of torque tachometers as well as other necessary data are all collected by data acquisition system, as shown in Fig. 7b. The experimental wind speed can be set manually by the PC program.

(a) = 1

(b) = 1.25

3.1.2. Measurement The dynamic measurement of torque and rotational speed is completed by the JC torque & tachometer shown in Fig. 8a. The measurement precision of speed and torque is±0.1% at rated torque of 20N  m. The two ends of the torque meter are connected by a rigid coupling with the output shaft of the pulley and the DC generator respectively. The measuring principle is shown as Fig. 8b. The basic principle of JC torque & tachometer is: through the elastic axis and two groups of magnetic signal generator, the measured torque and rotating speed are converted to a pair of AC signals with phase difference. The frequencies of these AC signals are the same and proportional to the rotating speed of the elastic axis, and the change of their phase difference is proportional to the measured torque. The AC signals are sent to the 44723A module by special shielded cable, and the exact value of torque and rotating speed can be obtained as follows:

(c) = 0.8

Fig. 6. Three different rotor diameter ratio schemes.

cases 13 and 14 is 1.25 times of those in cases 15 and 16. In addition, previous studies [11e13,18,29] have found that the distance between rotors in a certain range would have relatively weak influence on the CRWT performance. So in this paper, the distance is set to a constant 0.1 m as recommended in Ref. [13]. 3. Research methodology 3.1. Experimental method 3.1.1. Experiment setup The wind tunnel is employed to provide the perfect environment with stable wind speed. As shown in Fig. 7a, a honeycomb and screens covers the entrance of the wind tunnel, which make the flow conditioned. A 10:1 area contraction is designed behind the honeycomb to provide a uniform, low-turbulence-intensity approach velocity. The compressor is located at the end part of the tunnel to adjust the flow speed in the range of 10e70 m/s with adjustment accuracy of 0.1 m/s passing through the test section with parameters of 0.8 m  0.8 m  1.5 m. Those entire designs make the wind tunnel has the capacity to provide stable wind

TðN  mÞ ¼ TN 

nðrpmÞ ¼ 60 

f  f0 fN  f0

(2)

F Z

(3)

Where TðN  mÞ is the measured torque, nðrpmÞ is the measured rotating speed, TN is the full range torque value, fN is the full range frequency value, f0 is the 0 torque frequency value, f is the measured frequency, F is the rotating speed signal frequency, and Z is the number of signal gear teeth.

Table 1 Research cases. Cases

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

bfront/
brear/


30 30 500 1

30 30 400 0.8

30 30 500 1.25

30 60 500 1

30 60 400 0.8

30 60 500 1.25

60 60 500 1

60 60 400 0.8

60 60 500 1.25

60 30 500 1

60 30 400 0.8

60 30 500 1.25

30 Null 500 Null

60 Null 500 Null

30 Null 400 Null

60 Null 400 Null

Dfront/mm

d

L. Zhiqiang et al. / Energy 161 (2018) 939e954

943

(a)

Fig. 7. Wind tunnel experimental setup.

(a) the JC Torque & tachometer

(b) Schematic diagram of measurement principle

Fig. 8. Measurement of Torque and Rotating speed.

3.2. Numerical method 3.2.1. A brief introduction to LBM In order to further obtain the flow details of the CRWT models based on the aforementioned wind tunnel experiments, especially the wake flow and the aerodynamics interaction of the two rotors, the lattice Boltzmann method (LBM) is used for numerical simulation. The lattice Boltzmann method is a modern approach in Computational Fluid Dynamics which is often used to solve the incompressible, time-dependent Navier-Stokes equations numerically [36]. Compared with the molecular dynamics equation which directly describes the microscopic molecular motion and the N-S equation which describes the macroscopic continuum motion, LBM is a kind of meso-scopic method based on the Boltzmann equation (BE) and is valid over a wider range of flow physics than the N-S equations [37]. For this reason, it is an invaluable tool in fundamental research, as it keeps the cycle between the elaboration of a theory and the formulation of a corresponding numerical model short and its strength lies however in the ability to easily represent

complex flow-solid interactions between the geometry and surrounding fluid. At the same time, it has proven to be an efficient and convenient alternative to traditional solvers for a large variety of industrial problems, many open-source or commercial solvers based on LBM have been developed during past decade, such as Palabos, OpenLB and XFlow. And in this paper, the following equations were solved by XFlow 2016 edition with its adaptive Cartesian moving grid generator. 3.2.2. The Lattice Boltzmann equation The Boltzmann kinetic equation can be written as the following form:

  eq fs;i ðx; tÞ  fi ns ðx; tÞ; uðx; tÞ ts  ts Fs ðx; tÞ þ ms ns

fs*;i ðx; tÞ ¼ fs;i ðx; tÞ 

h

(4)

944

L. Zhiqiang et al. / Energy 161 (2018) 939e954

Where the pre-collisional populations fs;i and post-collisional

populations fs*;i are related through the streaming step:

fs;i ðx; tÞ ¼ fs*;i ðx  hci ; t  hÞ

(5)

The Eq. (4) is mainly used to describe the collision procedure of micro molecular and for numerical solving the macro fluid motion, based on it, the lattice gas (LG) automata was proposed [38e40], which can be considered a discrete particle kinetics utilizing a discrete lattice and discrete time. The evolution equation of the LG automata is written as follows:

ni ðX þ ei ，t þ 1Þ ¼ ni ðX，tÞ þ Ui ðnðX，tÞÞ;

ði ¼ 0，1…，MÞ (6)

Where: ei is the local particle velocity. The LBM is a further derivative of the LG method, in which the LBE is obtained from a discrete kinetic equation for the particle distribution function (PDF) [41]. Similar to the kinetic equation in the LG automata in Eq. (6), the LBE can be written as:

fi ðX þ ei Dx，t þ DtÞ ¼ fi ðX，tÞ þ Ui ðf ðX，tÞÞ;

ði

¼ 0，1…，MÞ

(7) Fig. 9. The velocities ei of the D3Q19 lattice.

Where fi is the particle velocity distribution function along the ith direction; Ui ¼ Ui ðf ðX，tÞÞ is the collision operator which represents the rate of change of fi resulting from collision. Dt and Dx are time and space increments, respectively. When Dx=Dt ¼ jei j, Eqs. (6) and (7) have the same discretization. Ui depends solely on the local distribution function. In the LBM, space is discretized in a way that is consistent with the kinetic equation. The coordinates of the nearest neighbor points around X is X þ ei . Frischet al [39,40] employed the Chapman-Enskog expansion to derive the macroscopic hydrodynamic equation from Eq. (7). Based on the BGK collision models [42], Qian [43,44] and Chen et al. [45,46] developed the LBM to solve the second-order N-S equations, and finally the approximate N-S equation described with LBE is obtained [47]:



r

vua þ Vb  ua ub vt



  ¼ Va p þ nVb  Va rub þ Vb rua

(8) Fig. 10. The computational domain.

3.2.3. Adaptive LBM for CRWT Wood et al. [48] and Deiterding et al. [49] used the D3Q19 model of LBM to simulate the velocity distribution of single rotor HAWT and the wind farms, which demonstrated the potential of LBM in wind turbine numerical simulation. The LBM-D3Q19 model recommended by them is also applied for the CRWT aerodynamics simulation in this paper. As shown in Fig. 9, the lattice velocities in Eq. (9) are defined as:

8 i ¼ 0; < 0; ei ¼ ð±1; 0; 0Þc; ð0; ±1; 0Þc; ð0; 0; ±1Þc; i ¼ 1; …; 6; : ð±1; ±1; 0Þc; ð±1; 0; ±1Þc; ð0; ±1; ±1Þc; i ¼ 7; …; 18

(9)

3.2.4. Computational domain and boundary conditions The computational domain in the same parameters of wind tunnel test section is shown in Fig. 10. The fluid material is defined as air with molecular weight of 28.996 g/mole, temperature of 288.15 K (the same with experiment condition), dynamic viscosity

of 1.7894  105 Pa s and density of 1.225kg/m3. The inlet outlet condition is defined as Velocity Inlet with X-direction wind speed of 10 m/s and Pressure Outlet with average static pressure of 0.1 MPa. The 6-DOF solver applied in this study uses the rotors’ forces and moments in order to compute the translational and angular motion of the center of gravity of the two aluminum rotors with density material of 2.7  103 kg/m3. The governing equation for the translational motion of the center of gravity is solved for in the inertial coordinate system:

1 X! ! f G vG¼ m

(10)

Where ! v G is the translational motion of the center of gravity, m is ! the mass, and f G is the force vector due to gravity. The angular motion of the rotor, ! u B , is more easily computed

L. Zhiqiang et al. / Energy 161 (2018) 939e954

using body coordinates:

X!  ! MB  ! u B ¼ L1 u B  L! uB

(11)

! Where L is the inertia tensor, M B is the moment vector of the body, ! and u B is the rigid body angular velocity vector. The moments are transformed from inertial to body coordinates using:

! ! M B ¼ RM G

(12)

Cq Cj 4 S4 Sq C j  C 4 Sj C 4 Sq C j þ S4 Sj

C q Sj S4 Sq Sj þ C 4 C j C 4 Sq Sj  S4 C j

(13)

Where, in generic terms, Cc ¼ cosðcÞ and Sc ¼ sinðcÞ. The angles 4, q, and j are Euler angles that represent the rotors rotation about the X axis. After the angular and the translational accelerations are computed from Eq. (8), the rates are derived by numerical integration [50]. The angular and translational velocities are used in the dynamic mesh calculations to update the rigid body position.

3.2.5. Turbulence modeling The LBM is proven as an accurate and efficient method for direct numerical simulations (DNS) of simple turbulent flows (LBE-DNS) [43,44], but compared with the traditional CFD methods, the efficiency of LBE-DNS is relatively low when considering the huge scale of engineering CFD problems. To improve the practicability of LBM, Yu et al. [37,51,52] developed LBM by introducing the Smagorinsky large eddy simulation (LES) into LBE-LES model on the basis of Eggels’ work [53]. By extending the Smagorinsky model to the D3Q19 multiplerelaxation-time (MRT) model [54], Yu et al. [52] overcame the deficiencies of the numerical instability and inaccuracy under implementing boundary conditions. The D3Q19 MRT-LBE with Smagorinsky model is described in the following form proposed by Yu et al. in Ref. [52] as:

fi ðX þ ei Dx，t þ DtÞ  fi ðX，tÞ ¼ Ui ðf ðX，tÞÞ

i h ¼ M1  b S  m  mðeqÞ (14)

where M is a Q  Q matrix which linearly transforms the distribution functions f 2V to the velocity moments m2M:

m ¼ M  f; f ¼ M1  m

(15)

The discrete velocities in Eq. (14) are the same with ones defined in Eq. (7). The r0 is the mean density in the system, which is usually set to 1. The relaxation matrix b S is diagonal in the space M:

b S ¼ diagð0;  s1 ; s2 ; 0; s4 ; 0; s4 ; 0; s4 ; s9 ; s2 ; s9 ; s2 ; s9 ; s9 ; s9 ; s16 ; s16 ; s 16 Þ ¼ diag 0; se ; sε ; 0; sq ; 0; sq ; 0; sq ; sv ; sp ; sv ; sp ; sv ; sv ; sv ; sm ; sm ; sm (16) The viscosity is:

(17)

For LES, n ¼ y0 þ yt , where y0 and yt are the molecular viscosity and eddy viscosity, respectively. In the Smagorinsky model [55], the eddy viscosity yt from the filtered strain rate tensor:

.  2 Sab ¼ va ub þ vb ua

(18)

A filter length scale Dx and the Smagorinsky constant CS are written as:

yt ¼ ðCS Dx Þ2 S; S :¼ 2S : S

3

Sq S4 C q 5 C4 Cq

  1 1 1 2 c dt  3 sy 2

pffiffiffiffiffiffiffiffiffiffiffiffi

Where, R represents the following transformation matrix:

2

n¼

945

(19)

Martínez-Tossas et al. [56] used this LES method to simulate the wind turbine wakes with LBM. The research work showed that the computational results coincide with the dynamic model when the value of the Smagorinsky coefficient CS is appropriate, but it would be hard to predict an exact CS value for each case. In XFlow solver used in this paper, the dynamic Smagorinsky model is applied. The dynamic Smagorinsky model is a modification of the previous model where CS with initial value of 0.12 is dynamically computed based on the information provided by the resolved scales of motion, and thus CS may vary in space and time. This could avoid the problem of manual setting CS uncertainty.

3.2.6. Dynamic mesh adaption For the CRWT blade surface boundary condition, Xu et al. [57] gave a detailed algorithm for the fluid-solid coupling process of the moving boundary in LBM. And for the applying of 6-DOF model, just as real situations, the blade surface of CRWT will accelerate the rotation along its axis when it is under the aerodynamic force. When the resistance moment of the rotating axis is constant, the rotor will gradually achieve the maximum speed determined by the design conditions of the blade. During this process, the surface and surrounding meshes of the blade will also rotate following the blade surface, that is, every new step of calculation will change the mesh structure of the last step. This requires adaptive dynamic mesh technology. When using conventional FVM or FDM to solve the continuum N-S equations in complex geometric structure, it is quite difficult to apply adaptive dynamic mesh technology. However, for the Cartesian grid based LBM which applies particle distribution functions (PDFs) and particle collision model, because the PDFs will be updated in every calculation step, it is very convenient to use adaptive dynamic mesh technology to deal with moving boundary. By adopting the Cartesian adaptive dynamic mesh refinement (SAMR) proposed by Berger et al. [58], a local dynamic mesh adaptation method for the forced rotation HAWT rotor was recommended by Wood et al. [48] and Deiterding et al. [49]. Neumann et al. [59] presented a new adaptive Lattice Boltzmann implementation within the Peano frame work which was not only adaptive to the boundary surface, but also adaptive to the structure of vortex. In this paper, the LBM solver and mesh generator XFlow is used to generate the adaptive wake (vortex) refinement mesh and solve the D3Q19 LBE-LES. The local dynamic mesh structure is shown in Fig. 11. Based on verifying mesh independent solutions, the number of initial (at computer step 0), refined (at computer step 1) and final meshes are shown in Table 2 and Fig. 12. The Layer Number in Table 2 are tags of different mesh scales and the minimum mesh scale near the rotor blade surface is about 0.05 mm, which can adjust the geometry surface well.

946

L. Zhiqiang et al. / Energy 161 (2018) 939e954

(a) X-Y Plane Mesh View

(b) Z-Y Plane Mesh View

(c) 3D Domain Mesh View

(d) Boundary Layer Mesh Structure at the Blade Root

(e) Dynamic Mesh Adaption Processing with the rotation of rotor Fig. 11. Dynamic mesh (lattice) structure.

Table 2 Mesh refinement result. Layer No.

1

2

3

4

5

6

7

8

Total

Initial Refined Final

130352 1106536 7301760

22362 75763 636296

9256 15000 177389

1823 924 62551

192 0 14992

16 0 1544

0 0 128

0 0 0

164001 1198223 8194660

(a) Initial Mesh Scale

(b) Refined Mesh Scale

(c) Final Mesh Scale

Fig. 12. Mesh scale.

4. Results and discussion 4.1. Numerical method validation The numerical results are validated by X-FOIL and experimental data. The X-FOIL is an interactive program for the design and analysis of subsonic isolated airfoils, which was first developed by Mark Drela at MIT as a design tool for the MIT Daedalus project in the 1980s. It was further developed in collaboration with Harold Youngren. The current version is 6.99, released in December 2013. Despite of its vintage, X-FOIL is able to provide the airfoil lift coefficient (CL) and drag coefficient (CD) prediction close to the experimental data in wide ranges of Re numbers [60,61], so it is still one of the best choices to validate the CFD result. And in this paper, the X-FOIL data of NACA4418 airfoil in chord length of 50 mm at

Re ¼ 33000 and attack angle (alpha) of 0e20
which have been verified by Fei-Bin Hsiao’s wind tunnel experiment [34] is adopted. As shown in the Fig. 13a, the CL and CD curves predicted by LBM are agreement with the X-FOIL data, the maximum error point of CL curves between X-FOIL and LBM data is at alpha ¼ 15
and the relative error is lower than 7.5%. The prediction trend of X-FOIL and LBM on the lift-drag ratio (CL/CD) is basically consistent in Fig. 13b, and the maximum error point of CL/CD is at alpha ¼ 15
and the relative error is lower than 25%, which is mainly caused by the error transfer. The most important is the maximum lift-drag ratio points predicted by X-FOIL and LBM are close, this indicates that the LBM’s prediction of the separation point in this study is accurate. For the case 1 and case 7 in Table 1, the comparison between the calculated and experimental values of the rotational speed of the rotors with 0-loads at different wind speeds is shown in Fig. 13c. It

L. Zhiqiang et al. / Energy 161 (2018) 939e954 The Comparison of CL&CD@Re=33000 Predicted by X-foil and LBM

1.0

CL X-foil CL LBM CD LBM

0.4

CL/CD

CL&CD

6

CD X-foil

0.6

X-foil LBM

5 4

400

200 100

5

10

15 alpha

20

25

30

(a)The comparison of CL&CD predicted by X-foil and LBM

0

-100 -200 -300

2

0

No.1 Front Sim No.1 Rear Sim No.1 Front Exp No.1 Rear Exp No.7 Front Sim No.7 Rear Sim No.7 Front Exp No.7 Rear Exp

The comparison of Experimental and LBM results

300

3

0.2 0.0

The Comparison of CL/CD@Re=33000 Predicted by X-foil and LBM

7

0.8

500

Rotational Speed (Rad/s)

1.2

947

-400 0

5

10

15

20

25

30

10

15

Wind speed (m/s)

20

alpha

(b)The comparison of CL/CD predicted by X-foil and LBM

(c)The comparison of Experimental and LBM results

Fig. 13. Numerical method validation.

can be seen that the calculated results are in good agreement with the experimental results, and the relative error is lower than 5%. Hence, the numerical method for CRWT simulation in this paper is reliable. 4.2. Performance 4.2.1. Power coefficient Based on the transient rotational speed and torque measured by the experiment, the power coefficient curves of each Case can be obtained according to the definition of shaft power:

P ¼T u

(20)

And the definition of power coefficient:

Cp ¼ P

 1 3 rv s 2

(21)

Where T is the measured torque, u is the measured rotational speed, Cp is the power coefficient, r is the air density in experimental environment, v is experimental wind speed, s is the windward area of the CRWT. Notice, the CRWT were installed in the wind tunnel experiment section, and the ratio of the CRWT’s windward area to the wind tunnel experiment section is:

sCRWT p  r2 ¼ 37:1% ¼ sEXP section LL

(22)

Because the wall of the wind tunnel can be seen as an external duct to the CRWT, the blocking effect in the experimental section can’t be ignored. Probably, wind tunnel experiment of wind turbine may lead to the enhancement of measured power coefficient as reported in Ref. [62e67]. In this paper, the area of wind tunnel experiment section is alternative for the CRWT’s windward area. The power coefficient ðCp Þ curves of the front and rear rotors in Cases 1e16 under the wind speed of 10 m/s are shown in Fig. 14. 1) The Cp of front rotors As shown in Fig. 13a, when Dfront > Drear, the maximum values of Cp ðCp:max Þ of Cases 1,3,4,6 and Control Case 13 (marked with the gray box) with bfront ¼ 30º are much larger than those of Cases 7,9,10,12 and Control Case 14 (marked with the black box) with bfront ¼ 60º. And when Dfront < Drear, the maximum values of Cp ðCp:max Þ of Cases 2, 5 and Control Case 15 (marked with the green

box) with bfront ¼ 30º are also far larger than ones of Cases 7, 11 and Control Case 16 with bfront ¼ 60º. In addition, when Dfront > Drear, the Cp:max values in Cases 1, 3, 4 and 6 are all significantly higher than one of Control Case 13. However, when Dfront < Drear, the Cp:max values in Cases 2, 5 are slightly lower than one of Control Case 15. This phenomenon indicates that the rear rotor can enhance the power output of the front rotor in CRWT with Dfront > Drear, and the higher the rear rotor Cp:max , the lower increase of the front rotor Cp:max. Meanwhile, the rear rotor will lead to the decrease of front rotor Cp:max, in CRWT with Dfront < Drear, and the higher the rear rotor Cp:max, the lower the decline of the front rotor Cp:max. Also, the Cp of the front rotor with 30
is higher than the one with 60
in all Cases, so the pitch angle of the front rotor should not be too large. According to Fig. 14a, the front rotor Cp:max in Cases 1, 3, 4 and 6 are:

Cp:max:Case6 > Cp:max:Case3 zCp:max:Casep4 > Cp:max:Case1

2) The Cp of rear rotors Compared with the power output of the front rotor, it can be seen from the Fig. 14 (b) that the law of rear rotor Cp is a little more complex. When Dfront < Drear, bfront ¼ 60
and brear ¼ 30
, the front rotor Cp of Case 11 is the lowest, but the rear rotor Cp of it is the highest among the Cases 1e16. And the rear rotor Cp:max of Case 11 is much higher than that of Control Case 13, and close to the front rotor Cp:max of Case 6. When Dfront ¼ Drear, the Cp curve of Case 10 is nearly overlapping with the Cp curve of Control Case 13. This phenomenon indicates that the front rotor has little effect on the Cp of rear rotor. As shown in Fig. 14b, the effect of pitch angle on the rear rotor Cp is the same as that of the front rotor, and the Cp of the rear rotor with 30
is also higher than the one with 60
in all Cases. In summary, the higher the front rotor Cp:max, the lower the rear rotor Cp:max. And what’s more special is that the front rotor in Case 11 with smaller Dfront and bfront ¼ 60º will remarkably improve the rear rotor Cp:max. This is a phenomenon worth paying attention to. 3) The total Cp of CRWT system In this paper, the relative tip speed ratio (TSR) of the two-stage rotors is defined as the reference TSR of the whole CRWT system’s power coefficient:

L. Zhiqiang et al. / Energy 161 (2018) 939e954

Power Coefficient (Cp) of Front Rotor

0.30

Max Cp point

0.25

Cases of large rotors with pitch angle of 30 degrees

0.20

Cp

No.1 No.2 No.3 No.4 No.5 No.6 No.7 No.8 No.9 No.10 No.11 No.12 No.13 No.14 No.15 No.16

under wind speed 10m/s

0.15 0.10 Cases of small rotors with pitch angle of 30 degrees

0.05 0.00

0

Cases of rotors with pitch angle of 60 degrees

2

4 6 8 10 Tip Speed Ratio (TSR)

No.1 No.2 No.3 No.4 No.5 No.6 No.7 No.8 No.9 No.10 No.11 No.12 No.13 No.14 No.15 No.16

Power Coefficient (Cp) of Rear Rotors

0.30

under wind speed10m/s

0.25 Cases of large rotors with pitch angle of 30 degrees

0.20

Cp

948

0.15 0.10 Cases of small rotors with pitch angle of 30 degrees

0.05 0.00

12

0

2

Cases of rotors with pitch angle of 60 degrees

4 6 8 10 Tip Speed Ratio(TSR)

12

(a)Power coefficient ( ) curves of the front rotors. The curves

(b)Power coefficient ( ) curves of the rear rotors. The curves

marked with the gray box belong to Cases 1, 3, 4, 6 and Control

marked with the blue box belong to Cases 1, 2, 10

Case 13. The curves marked with the green box belong to Cases

Control Case 13. The curves marked with green box belong to

2, 5 and Control Case 15. The curves marked with the black box

Cases 3, 12 and Control Case 15. The curves marked with

belong to Cases 7~12 and Control cases 14 and 16.

black box belong to Cases 4~9 and Control Cases 14 and 16.

11 and

Fig. 14. C p curves of the front and rear rotors under the wind speed of 10 m/s.

Power Coefficient (Cp) of CRWT under wind speed=10m/s

0.35

Cp.Max=34.2%

0.30

Cp

0.25 0.20 0.15 0.10

No.1 No.2 No.3 No.4 No.5 No.6 No.7 No.8 No.9 No.10 No.11 No.12 No.13 No.14 No.15 No.16

0.05 0.00

0

2

4

6

8

10 12 14 16 18 20 22

Tip Speed Ratio (TSR) Cp curves of CRWT under the wind speed of 10m/s Fig. 15. Cp curves of the CRWT system.

TSRrel ¼

ufront Rfront  urear Rrear v

(23)

Due to the counter rotating of the two-stage rotors, the relative TSR is the sum of the absolute values of the TSRs of the two-stage rotors. The whole CRWT system’s power coefficient ðCp Þ curves under wind speed of 10 m/s is shown in Fig. 15. From the Fig. 15, when the pitch angle b of the large diameter rotor is 30
, the Cp:max of the Cases 1, 2, 3, 4, 6, and 11 are all higher

than the Cp:max of the Control Case 13. The Cp:max of the Case 1 with Dfront ¼ Drear and bfront ¼ brear ¼ 30
is the highest, reaching 34.2%, which is about 54.8% higher than that of the Control Case 13 with single rotor. Compared with Fig. 14, although the Cp:max of the front and rear rotors in Case 1are both not the highest, the Cp:max of the CRWT system is the highest. According to Fig. 14, the Cp:max in Cases 1, 2, 3, 4, 6, and 11 are:

Cp:max:Case1 > Cp:max:Case11 > Cp:max:Case2 zCp:max:Case6 > Cp:max:Case3 > Cp:max:Case4 > Cp:max:Case13 zCp:max:Case10 Based on the pitch angles of front and rear rotors and the total power coefficient, the working principles of CRWT studied in this paper can be categorized into three types: Cases 1, 2 and 3 belong to the first type, which is named “MainMain type” in this paper. Their common character is that the pitch angles of two-stage rotors are both 30
. Two-stage rotors both play an important role in power output of CRWT, and the differences of the total power coefficient mainly depend on the diameter ratio of two stage rotors. Cases 4, 5 and 6 belong to the second type, which is named “Main-Vice type” in this paper. The only difference between the first type and the second one is the pitch angle of the rear rotor. From the Fig. 14b, it can be seen that the rear rotor with brear ¼ 60º hardly make direct contribution to the total power coefficient. However, the rotation of the rear rotor can enhance the power output of the front rotor, and indirectly lead to the increase of total power coefficient. Compared with d  1 in Case 4 or Case 5, the smaller rear rotor diameter (d > 1) in Case 6 is more favorable to the power output of the front rotor. Hence, in this type, big and small rotors are called the main and vice rotors respectively. Cases 10, 11 and 12 belong to the third type, which is named “Vice-Main type” in this paper. The difference between this type

L. Zhiqiang et al. / Energy 161 (2018) 939e954

and the other types is the pitch angle of the front rotor. As shown in Fig. 15, the front rotor with bfront ¼ 60º is almost no power output in CRWT. Compared with d  1 in Case 10 and Case 12, the smaller front rotor diameter (d < 1) in Case 11 is more favorable to the power output of the rear rotor. In this type, big and small rotors are still called the main and vice rotors. In summary, the CRWT can enhance the total power coefficient in a certain range of wind speed by two schemes. One scheme is to use the second stage rotor to further utilize the wind energy wasted by the first stage rotor. Experimental results have shown that the equal diameter counter rotating rotors with certain pitch angle can achieve a maximum power coefficient in a wider range of speed, just like Case 1 in Fig. 15. The other scheme is to use the vice rotor to enhance the power output of the main rotor, whether the vice rotor is the front rotor or rear one. The previous CRWT theoretical model proposed by Newman [5] is based on the assumption that the front rotor won’t be influenced by the rear rotor. This study has proved that the rear rotor can enhance the power output of the front rotor, just like Case 6 in Fig. 14. The mechanism of the power coefficient enhanced by different schemes of CRWT will be discussed in detail according to the LBM simulation results in aerodynamics research part of this paper.

4.3. Aerodynamics 4.3.1. The flow mechanism of CRWT As mentioned above, Cases 1, Case 4&6 and Case 10&11 represent three different schemes to improve the power coefficient of CRWT. Considering that their Cp curves shown in Fig. 15 are almost similar under different wind speeds, their flow mechanism is studied in typical working conditions of 10 m/s wind speed. Fig. 16 shows the LBM simulated transient velocity contours in five cylinder sections along the flow direction of Case 1, 4, 6, 10, 11 and the Control Case 13 with main rotor (case 1 the front rotor) rotational speed of 105e108 rad/s. The section A and D are for presenting the rotors’ inflow conditions, the section B and E are for presenting the rotors’ outflow conditions, and the section C is for presenting the flow characteristics between two rotors. Taking Case 1 as an example, it can be seen that the low velocity zone grows larger gradually during the five sections which indicates that the wind speed is slowed down by the two rotors stage by stage. Compared with the sections B, C and D, the low velocity zone is divided into two different parts in Section E, which are the center area and the blade tip area. Obviously, the center low velocity zone is partly originated from the wake flow shed from the front rotor’s hub and partly generated by the blade root. While the low velocity area near the blade tip is totally caused by the “recycling” wind energy procedure of the rear rotors. Also, by comparing the central areas of Sections D and E, although the blending of upstream wake slows down the center area velocity, there’s still surplus wind energy that can be gained by the central part of the rear rotor. This indicates the central part of the rotor blade is still necessary for the full utilization of wind energy, which argued with the theoretical model proposed by Chantharasenawong et al. [17] and Jung et al. [18]. As shown in Fig. 15, the transient velocity contours in Sections A, B and C of Case 1, Case 4, Case 6 are basically consistent and slightly lower than those of the Control Case 13. Due to the rotation of the rear rotors, the center velocity contours in section D of Cases 4 and 6 looks more “chaotic” than one of Case 13, and the chaos area in section E of Cases 4 and 6 grows larger than one of Case 1 by the agitation of rear rotors. However, the wind speed from the middle to the tip of the rear rotor blades almost no change, the reason is that the success design of little Cp rotor with b ¼ 60
makes the

949

rotor more like an agitator than a wind turbine. Previous studies have shown that the Cp of main rotor can be improved with the existence of this “agitator”. Here, we speculate that the agitation effect in the wake flow probably reduces the pressure of downstream flow and thus accelerates the upstream flow, just like a diffuser. For Cases 10 and 11, the transient velocity contours are unable to provide very convincing information, so it is necessary to make further research. To further analyze the aerodynamic mechanism of these three types of CRWTs, three important flow parameters time averaged after 1.5s time steps (0.001s/time step) at above five sections are processed by surface integral average using following method:

Z Z 0

S

4ðni ÞdS S

(24)

Where the 4ðni Þ are the flow parameters at each lattice node ni in the section, the S is the area of the section. And the processed flow parameters are Vx (Velocity in X-axis direction), P* (Total pressure), and TI (Turbulence intensity). All these parameters are normalized with initial values at section A, and Fig. 17 shows the distribution of these parameters in each section along X direction. As shown in Fig. 17a, the parameter Vx curves shows the integral mean velocity loss along the X-axis. Compared with other Cases, the curve slope of Case 1 is the greatest, which presents the maximum velocity loss rate. The velocity loss rates of Cases 4&6 drop faster than the one of Control Case 13, but slower than that of Case 1 after section C. This phenomenon indicates that the rotation of rear rotors will lead to velocity loss more or less when the wind blows through the same front rotor with the same design. This further proves the rear rotor can influence the front rotor which was mistakenly neglected in Newman’s assumption [5] and many other relative studies. The velocity loss rates of Cases 10 &11 drop slower than that of Control Case 13 alone the first three sections, and accelerate decline along the last two sections. This phenomenon indicates large pitch angle of the front rotor is not beneficial to the utilization of wind energy. The total pressure P* reflects the kinetic energy contained in fluid and is usually used as a useful parameter to measure the potential of fluid working. Therefore, the P* loss of the rotor can be used to estimate the efficiency of wind energy utilization, and the parameter of exhaust expansion ratio ðpfront =prear Þ which is commonly used in gas turbine design can also be introduced to evaluate the performance of the wind turbine. As shown in Fig. 17b, the expansion ratios of different Cases have been normalized and can be easily compared by observing the points’ ordinate. It can be seen that the pressure drop of the front rotor with 30
pitch angle is larger than that of 60
pitch angle, and the front rotor in Case 6 has the maximum expansion ratio and can collect most wind energy among all the front rotors in compared cases, which coincide with the foregoing conclusion that the front rotor in Case 6 has the best front rotor performance. In Cases 10 &11, the pressure drops of the front rotors are negligible, compared with the significant pressure drops of the rear rotors, which also coincide with the previous conclusion that the power outputs of the rear rotors are enhanced by the rotation of the front rotors. For Case 1, the p* in section E is the lowest, so the total Cp of Case 1 is the highest. A noteworthy interesting phenomenon is that the pressure drop of the front or rear rotors of Case 1 is slightly less

950

L. Zhiqiang et al. / Energy 161 (2018) 939e954

(a)Section Locations Section A

Section B

(b)Uniformed Velocity Magnitude in Contours Section C

Section D

Section E

Case

Case1

Case4

Case6

Case 10

Case 11

Case 13

Fig. 16. The Velocity Contour of different schemes CRWT at five X-axis sections.

L. Zhiqiang et al. / Energy 161 (2018) 939e954

front rotor

B

0.95

C

1.0

D

Vx/U

0.90 No.13 No.1 No.4 No.6 No.10 No.11

0.80 0.75 0.70 0.0

0.2

rear rotor 0.4 0.6 x/L

0.8

B

C

No.13 No.1 No.4 No.6 No.10 No.11

D

0.8

E

0.85

A

0.6 0.4

6 5 4 3

front rotor

rear rotor

E

0.2

2

0.0

1

1.0

(a)Normalized Vx

No.13 No.1 No.4 No.6 No.10 No.11

TI

A

p*/P*

1.00

951

0.0

0.2

0.4

x/L

0.6

0.8

1.0

(b)Normalized Total Pressure

rear rotor

E

front rotor

D C

B A 0.0

0.2

0.4

x/L

0.6

0.8

1.0

(c)Normalized Turbulence Intensity

Fig. 17. The Normalized Integral Averaging Flow Parameters of different schemes CRWT at five X-axis sections.

than that of Case 13 with the same pitch angle, but the total pressure drop of Case 1 is higher than that of case 13. This phenomenon further confirms that the CRWT can increase the power coefficient of HAWT effectively. As shown in Fig. 17 (c), due to the same front rotor design, the turbulence intensity of Case 1, Case 4, Case 6 and Case 13 obviously decline from Section B to C and remain almost unchanged from Section C to D, which seems to be consistent. A possible explanation is that the counter-rotation of front and rear rotors may lead to a certain constraint on the upstream flow which not only restrains the development of turbulence but also changes the outflow condition of the upstream rotors. The turbulence intensity of Case 10 and 11 slightly grows from Section A to C. From Section C to E, the turbulence intensity of Cases 10 and 11 rapidly increases, and this will effectively restrain the blade surface flow separation and improve the Cp of rear rotors in Cases 10 and 11. A detail discussion about this mechanism is proposed with the wake interaction

analysis in the following. 4.3.2. The wake influence In the previous numerical method validation part, the maximum lift-drag ratio point of airfoil NACA4418 under the test condition Re ¼ 33000 with alpha ¼ 6e7.5ºhas shown that the flow around NACA4418 airfoil blade surface will separate when the local attack angle alpha is over 7.5
. As shown in Fig. 18, under the rotation condition, the alpha angle can be calculated by the local velocity angle 4 and the pitch angle b:

a¼4b

And the local velocity angle 4 is synthesized from two velocity vectors of wind speed Vx and rotational speed u  R, which can be easily calculated by:

4 ¼ arctan

Fig. 18. The velocity vectors and related angles.

(25)

jVx j ju  Rj

(26)

Moreover, under the rotation condition, there’s a delayed stall phenomenon which was first proposed by Himmelskamp [68] and studied by many other scholars. In this paper, the delayed stall theory proposed by Wang [69] and Du [70] is applied and the blade section at R ¼ 0.15 m and u ¼ 74 76rad=s is chosen as Fig. 19 shown. In order to explain the reason why the output power of the front rotors is very small in Case 10 & 11, the special separation mechanism of Vice-Main rotor scheme CRWT is studied. During the simulation, the wind speed Vx is approximately equal to 10 m/s, and the local velocity angles and alpha angles are set as 4 ¼ 41:25
, a ¼ 11:25
in Case 13, 4 ¼ 41:53
, a ¼ 11:53
in Case 10, and 4 ¼ 41:82
, a ¼ 11:82
in Case 11. As shown in Fig. 19, the periodic large vortex sheds from the trailing edge of the blade in Case 13, which will generate the additional drag force known as vortex drag (Bet, 1925) or induced drag (Plandtl, 1918), and lead to the loss of blade efficiency. Compared with Case 13, the flow separation in Case 10 is weaker, and there is almost no flow separation in Case 11. By analyzing the difference among the three Cases, the most possible reason is that the wake flow of the front rotor can enhance the downstream turbulence intensity as shown in Fig. 17c and the increased turbulence intensity can improve the momentum and energy transfer, which will lead to the separation restrained in the blade surface. Furthermore, as shown in Fig. 17c, between sections D & E in Case 4 & 6 and between sections A & B in Case 10 & 11, the smaller rotor with b ¼ 60
can always generate stronger turbulence than

952

L. Zhiqiang et al. / Energy 161 (2018) 939e954

hand, the design size of the micro rotor also amplifies the negative effects of blade roughness and transmission loss, thus reducing the overall efficiency. 2. The design of CRWT rotors is more complex than the design of single rotor HAWT. The BEM method which performs well in the single rotor HAWT designing should be corrected in the CRWT design works, considering the interaction of the two stage rotors. The design of generator and mechanical structure of CRWT are also difficult, which limiting the commercial application of CRWT. And for Micro CRWT, the diameter of rotors is relative small, so at least one rotor should rotate at a high absolute rotational speed to reach the design TSR, which would be negative for the equipment life and noise control of Micro CRWT. 3. The rotation of front rotor will lead to strong fluctuating of the rear rotor’s torque output in almost all CRWT schemes, as shown in Fig. 20. This special phenomenon indicates that although the rotation of the front rotor would delay or inhibit the flow separation of the rotor blade surface, it will create an additional high frequency aerodynamic loads on the rear rotor. Obviously, this is unfavorable to the stability and power output of Micro CRWT. However, considering that Micro CRWT is mainly used in lowcost distributed generation, these shortcomings are still acceptable.

Fig. 19. The velocity vectors in the middle section of rotors.

the larger rotor with b ¼ 60
. This may be because the smaller rotors rotate faster than the larger rotors (usmaller z69rad=s , ularger z33rad=s, in above presented cases). 4.4. Limitations As mentioned above, the different CRWT schemes are all effective for improving power coefficient of wind turbines, but it is still necessary to pay attention to the following limitations of the Micro CRWT. 1. The maximum system power coefficient among all of the CRWTs studied in this paper is only 34.2%, which is lower than that of the current commercial wind turbines. One reason is the use of non-torsional and non-optimized straight blades, on the other

wind speed=10m/s

1.5 1.0

Torque (N*m)

1# front torque 1# rear torque 2# front torque 2# rear torque 3# front torque 3# rear torque 13# torque 13#contra torque 15# torque 15# contra torque

Fluctuating region

0.5 0.0 -0.5 -1.0 -1.5

0

5

10

15

20

25

30

35

Time(s) Fig. 20. The fluctuating of the rear rotor torque output.

40

5. Conclusion To understand the CRWT better, the performance and the interaction between two-stage rotors in several micro CRWTs are studied using wind tunnel experiment and lattice Boltzmann method (LBM) in this paper. The fluid dynamics mechanism of CRWT improving the power coefficient is discussed, some conclusions are obtained as follows: 1. As a direct solving N-S equations method, LBM is high-precision and reliable. With the implement of LES and the Cartesian adaptive dynamic mesh, it can be well applied to a series of fluid-solid coupling and rigid body forced motion researches, such as the CRWT CFD studies in this paper. While in practice, it is also found that its solution scale is very large and its computing efficiency is relatively low, so it still needs further improvement for engineering CFD applying. 2. The counter rotating scheme is proven again for improving the wind power coefficient of HAWT in this paper, and the micro CRWT can enhance the power coefficient by two schemes. One scheme is to use the second stage rotor to further utilize the wind energy wasted by the first stage rotor. This scheme has been studied in many previous researches, and the experimental results presented in this paper give further evidence for its effect that the equal diameter counter rotating rotors with certain pitch angle can achieve a higher power coefficient in a wider range of speed than single rotor HAWT. 3. The other scheme is to use the vice rotor to enhance the power output of the main rotor, whether the vice rotor is the front rotor or rear one. For the Main-Vice type, experimental and LBM numerical results have shown that the power coefficient of upstream Main rotor can be enhanced by the downstream Vice rotor. This is because of the Vice rotor’s “agitation effect” decreasing the total pressure in downstream flow, which will lead to the total pressure drop in Main-Vice type CRWT higher than the single rotor HAWT. For the Vice-Main type, the Main rotor power coefficient improvement mechanism is relatively more complex: The front Vice rotor rotation increases the turbulence intensity in the downstream wake flow, which will lead

L. Zhiqiang et al. / Energy 161 (2018) 939e954

to the enhancement of fluid momentum and energy transfer, and restrain the blade surface separation. Thus the blade’s liftdrag ratio increase and the Main rotor power coefficient at certain rotational speed is improved. 4. In addition, previous CRWT researches have always considered that the rear rotor have no effect on front rotor. In this paper, the rear rotor influence on front rotor has been found. For the MainMain type, the rear rotor will have a negative effect on the front rotor by decreasing the total pressure drop between the up and downstream of front rotor. This will lead to the front rotor power coefficient loss at certain rotational speed. On the contrary, the rear rotor improve the front rotor power coefficient in MainVice type. The importance of this discovery is that, on the one hand, it points out the improper assumptions of the Newman’s CRWT theoretical model [5]. On the other hand, the coupling effect of the rear rotor and front rotor should be fully considered in the process of the CRWT design in the future. However, the quantitative rule of this effect has not yet been studied in detail by this paper works, meanwhile, under the restriction of experimental conditions and time, we did not test the performance of Micro CRWT under natural random wind speed conditions. So it remain to be further explored and studied in future research. Acknowledgements This work was supported by the Chinese National Natural Science Fund [grant numbers: 51276008]. References [1] Polinder H, Bang D, Van Rooij RPJOM, McDonald AS, Mueller MA. 10 MW wind turbine direct-drive generator design with pitch or active speed stall control. In: Electric machines & drives conference, 2007. IEMDC’07. IEEE international, vol. 2. IEEE; 2007, May. p. 1390e5. [2] Mithraratne N. Roof-top wind turbines for microgeneration in urban houses in New Zealand. Energy Build 2009;41(10):1013e8. [3] Peacock AD, Jenkins D, Ahadzi M, Berry A, Turan S. Micro wind turbines in the UK domestic sector. Energy Build 2008;40(7):1324e33. [4] Newman BG. Actuator-disc theory for vertical-axis wind turbines. J Wind Eng Ind Aerod 1983;15(1e3):347e55. [5] Newman BG. Multiple actuator-disc theory for wind turbines. J Wind Eng Ind Aerod 1986;24(3):215e25. [6] Ushiyama I, Shimota T, Miura Y. An experimental study of the two-staged wind turbines. Renew Energy 1996;9(1e4):909e12. [7] Appa, K. (2000). U.S. Patent No. 6,127,739. Washington, DC: U.S. Patent and Trademark Office. [8] Debleser, Y. (2003). U.S. Patent No. 6,504,260. Washington, DC: U.S. Patent and Trademark Office. [9] Stephens, B. K., & Price, D. (2010). U.S. Patent application No. 12/889,670. [10] Li Z, Wu Y. A micro wind turbine with the counter-rotating rotors. 2013. CN203098139U. [11] Li Z, Zhang Y. Numerical research of wind turbine with two contrary-turning impellers. In: Conference on China technological development of renewable energy source; 2010, December. [12] Zhiqiang Li, Yong Wang, Yusheng Zhang. Experimental study on interaction of the counter-rotating rotors in a wind turbine. In: Materials for renewable energy & environment (ICMREE), 2011 international conference on, vol. 1. IEEE; 2011, May. p. 570e4. [13] Li Z, Wang Y, Xiao H. Experimental study on structures of counter-rotating wind turbines. In: Materials for renewable energy and environment (ICMREE), 2013 international conference on, vol. 1. IEEE; 2014, August. p. 368e72. [14] Gregg JR, Merchant JS, Treuren KWV, Gravagne IA. Experimental analysis of a counter-rotating wind turbine. ASME 2009 international mechanical engineering congress and exposition, vol. 31; 2009. p. 1016. [15] Habash RWY, Groza V, Yang Y, Blouin C, Guillemette P. Performance of a contrarotating small wind energy converter. Int Sch Res Not 2014;2011(2). [16] Mitulet A, Oprina G, Chihaia RA, Nicolaie S, Nedelcu A, Popescu M. Wind tunnel testing for a new experimental model of counter-rotating wind turbine. In: Daaam international symposium on intelligent manufacturing and automation, Daaamvol. 100; 2015. p. 1141e9. [17] Chantharasenawong C, Suwantragul B, Ruangwiset A. Axial momentum theory for turbines with Co-axial counter rotating rotors. In: Commemorative international conference of the occasion of the 4th cycle anniversary of KMUTT; 2008, December.

953

[18] Jung SN, No TS, Ryu KW. Aerodynamic performance prediction of a 30kW counter-rotating wind turbine system. Renew Energy 2005;30(5):631e44. [19] No TS, Kim JE, Moon JH, Kim SJ. Modeling, control, and simulation of dual rotor wind turbine generator system. Renew Energy 2009;34(10):2124e32. [20] Hwang B, Lee S, Lee S. Optimization of a counter-rotating wind turbine using the blade element and momentum theory. J Renew Sustain Energy 2013;5(5), 052013. [21] Huang B, Shi J, Wei X. Model testing of a series of counter-rotating type Horizontal-Axis tidal turbines with 500 mm diameter. J Eng Gas Turbines Power 2017;139(10), 102602. [22] Haenler M, Ritschel U, Warnke I. Systematic modelling of wind turbine dynamics and earthquake loads on wind turbines. In: European wind energy conference and exhibition. Athens, Greece: European Wind Energy Association; 2006, February. p. 1e6. [23] Jonkman J, Butterfield S, Musial W, Scott G. Definition of a 5-MW reference wind turbine for offshore system development (No. NREL/TP-500-38060). Golden, CO: National Renewable Energy Laboratory (NREL); 2009. [24] Asr MT, Nezhad EZ, Mustapha F, Wiriadidjaja S. Study on start-up characteristics of H-Darrieus vertical axis wind turbines comprising NACA 4-digit series blade airfoils. Energy 2016;112:528e37. [25] Courbois A, Flamand O, Toularastel JL, Ferrant P, Rousset JM. Applying relevant wind generation techniques to the case of floating wind turbines. In: Sixth european and african conference on wind engineering. Nantes, France: EACWE; 2011. p. 7e13. July. [26] Cheboxarov VV, Cheboxarov VV. Research of behavior of large rotary pontoon of offshore wind turbine in waves. In: The ninth ISOPE Pacific/Asia offshore mechanics symposium. International Society of Offshore and Polar Engineers; 2010, January. [27] Makrygiannis GK, Dimas AA. Numerical simulation of fluid-structure interaction between sea waves and a spar-buoy wind turbine platform. In: 8th gracm international congress on computational mechanics Volos; 2015. 12 Julye15 July. [28] Asr MT, Osloob R, Mustapha F. Double-stage H-Darrieus wind turbine-rotor aerodynamics. Appl Mech Mater 2016:829. [29] Lanzafame R, Messina M. Design and performance of a double-pitch wind turbine with non-twisted blades. Renew Energy 2009;34(5):1413e20. ndez Juan, Greiner David. Wind blade chord and twist angle optimization [30] Me by using genetic algorithms. In: Topping B, Montero G, Montenegro R, editors. Proceedings of the fifth international conference on engineering computational technology, vol. 6. Las Palmas de Gran Canaria, Spain, Sept: Civil-Comp Press; 2006. [31] Song Q. Design, fabrication, and testing of a new small wind turbine blade. 2012. Doctoral dissertation. [32] Chaudhary MK, Roy A. Design & optimization of a small wind turbine blade for operation at low wind speed. World J Eng 2015;12(1):83e94. [33] Sirigu G, Cassaro M, Battipede M, Gili P, Frulla G. Wind generator innovative blade design: variable twist and start-up control. Int J Mech 2016;10:53e61. [34] Hsiao FB, Bai CJ, Chong WT. The performance test of three different horizontal axis wind turbine (HAWT) blade shapes using experimental and numerical methods. Energies 2013;6(6):2784e803. [35] Singh RK, Ahmed MR. Blade design and performance testing of a small wind turbine rotor for low wind speed applications. Renew Energy 2013;50:812e9. [36] Chen S, Doolen GD. Lattice Boltzmann method for fluid flows. Annu Rev Fluid Mech 1998;30(1):329e64. [37] Yu H, Girimaji SS, Luo LS. DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method. J Comput Phys 2005;209(2):599e616. [38] Hardy J, De Pazzis O, Pomeau Y. Molecular dynamics of a classical lattice gas: transport properties and time correlation functions. Phys Rev 1976;13(5): 1949. [39] Frisch U, Hasslacher B, Pomeau Y. Lattice-gas automata for the Navier-Stokes equation. Phys Rev Lett 1986;56(14):1505. [40] Frisch U, d’Humieres D, Hasslacher B, Lallemand P, Pomeau Y, Rivet JP. Lattice gas hydrodynamics in two and three dimensions (No. LA-UR-87e2524; CONF8610281-2). NM (USA); Observatoire de Nice, 06 (France); Ecole Normale Superieure, 75-Paris (France): Los Alamos National Lab.; 1986. [41] Wolfram S. Cellular automata and complexity: collected papers, vol. 1. Reading, MA: Addison-Wesley; 1994. [42] Bhatnagar PL, Gross EP, Krook M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys Rev 1954;94(3):511. [43] Qian YH. Lattice gas and lattice kinetic theory applied to the Navier-Stokes equations. Doktorarbeit. Paris: Universite Pierre et Marie Curie; 1990. res D, Lallemand P. Lattice BGK models for Navier-Stokes [44] Qian YH, d’Humie equation. EPL (Europhysics Letters) 1992;17(6):479. [45] Chen S, Chen H, Martnez D, Matthaeus W. Lattice Boltzmann model for simulation of magnetohydrodynamics. Phys Rev Lett 1991;67(27):3776. [46] Chen S, Wang Z, Shan X, Doolen GD. Lattice Boltzmann computational fluid dynamics in three dimensions. J Stat Phys 1992;68(3):379e400. [47] Qian YH, Orszag SA. Lattice BGK models for the Navier-Stokes equation: nonlinear deviation in compressible regimes. EPL (Europhysics Letters) 1993;21(3):255. [48] Wood Stephen L, Deiterding Ralf. A lattice Boltzmann method for horizontal axis wind turbine simulation at 14th International Conference on Wind Engineering. 2015. p. 18. Brazil. 21-26 Jun 2015.

954

L. Zhiqiang et al. / Energy 161 (2018) 939e954

[49] Deiterding R, Wood SL. Predictive wind turbine simulation with an adaptive lattice Boltzmann method for moving boundaries. In: Journal of Physics: Conference Series, vol. 753. IOP Publishing; 2016, September. p. 082005. No. 8. [50] Sokolichen A, Eigenberger G, Lapin A. Simulation of buoyancy driver bubbly flow: established simplifications and open questions. J Rev AIChE J January 2004;50(1):24e45. [51] Yu D, Mei R, Luo LS, Wei S. Viscous flow computations with the method of lattice Boltzmann equation. Prog Aero Sci 2003;39(5):329e67. [52] Yu H, Luo LS, Girimaji SS. LES of turbulent square jet flow using an MRT lattice Boltzmann model. Comput Fluids 2006;35(8e9):957e65. [53] Eggels JGM. Direct and large-eddy simulation of turbulent fluid flow using the lattice-Boltzmann scheme. Int J Heat Fluid Flow 1996;17(3):307e23. res D. Multipleerelaxationetime lattice Boltzmann models in three [54] D’Humie dimensions. Phil Trans 2002;360(1792):437. [55] Smagorinsky J. General circulation experiments with the primitive equations: I. The basic experiment. Mon Weather Rev 1963;91(3):99e164. [56] Martínez-Tossas LA, Churchfield MJ, Meneveau C. Large eddy simulation of wind turbine wakes: detailed comparisons of two codes focusing on effects of numerics and subgrid modeling. In: Journal of Physics: Conference Series, vol. 625. IOP Publishing; 2015. p. 012024. No. 1. [57] Xu L. A novel moving boundary condition based on Chapman-Enskog expansion with the lattice Boltzmann method. Doctoral dissertation. University of Pittsburgh; 2014. [58] Berger MJ, Colella P. Local adaptive mesh refinement for shock hydrodynamics. J Comput Phys 1989;82(1):64e84. [59] Neumann P, Neckel T. A dynamic mesh refinement technique for Lattice Boltzmann simulations on octree-like grids. Springer-Verlag New York, Inc; 2013. [60] Drela M. Xfoil: an analysis and design system for low Reynolds number

airfoils. Lect Notes Eng 1989;54:1e12. [61] Morgado J, Vizinho R, Silvestre MAR, P ascoa JC. Xfoil vs, cfd performance predictions for high lift low Reynolds number airfoils. Aero Sci Technol 2016;52:207e14. [62] Ohya Y, Karasudani T, Sakurai A. Development of high-performance wind turbine with brimmed diffuser. Nihon Koku Uchu Gakkai ronbunshu ¼ Journal of the Japan Society for Aeronautical and Space Sciences 2003;50(587): 477e82. [63] Ohya Yuji, Karasudani Takashi, Sakurai Akira, Inoue Masahiro. Development of high-performance wind turbine with a brimmed-diffuser part 2. J Jpn Soc Aeronaut Spaceences 2004;52(604):210e3. [64] Abe K, Nishida M, Sakurai A, Ohya Y, Kihara H, Wada E, et al. Experimental and numerical investigations of flow fields behind a small wind turbine with a flanged diffuser. J Wind Eng Ind Aerod 2005;93(12):951e70. [65] Ohya Y, Karasudani T, Sakurai A, Abe KI, Inoue M. Development of a shrouded wind turbine with a flanged diffuser. J Wind Eng Ind Aerod 2008;96(5): 524e39. [66] Kosasih B, Tondelli A. Experimental study of shrouded micro-wind turbine *. Procedia Eng 2012;49(29):92e8. borah Aline Tavares Dias do Rio Vaz, Mesquita ALA, Vaz JRP, Blanco CJC, [67] De Pinho JT. An extension of the blade element momentum method applied to diffuser augmented wind turbines. Energy Convers Manag 2014;87:1116e23. [68] (s1929). The rotol contra-rotating airscrew. Aircraft Eng Aero Technol. [69] Wang Q, Xu Y, Song JJ, Li CF, Ren PF, Xu JZ. 3d stall delay effect modeling and aerodynamic analysis of swept-blade wind turbine. J Renew Sustain Energy 2013;5(6):71e5. [70] Dumitrache A, Cardos V, Dumitrescu H. Delayed stall modeling of the rotating blades. Aip Conference, vol. 1281. American Institute of Physics; 2010. p. 1848e50. ̄ ̄

̄

̄