- Email: [email protected]
Energy extraction of two flapping foils with tandem configurations and vortex interactions
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Energy extraction of two flapping foils with tandem configurations and vortex interactions G.D. Xu∗, W.H. Xu College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China
a r t i c l e
i n f o
Keywords: Flapping foils Tandem configurations Energy extraction Global phase shift Efficiency
a b s t r a c t The energy extraction of two flapping foils with tandem configurations has been studied based on velocity potential theory. The interactions of the hind foil and the vortices of the fore foil with typical oscillatory motions are investigated. The global phase shift, which combines the longitudinal distance and the phase shift of the motions of the two foils, is introduced as the indicating parameter of the interaction modes. The effects of longitudinal distance and phase shift of the motions on the energy extraction efficiency have been analysed. The most favourable global phase shift and the corresponding vortex pattern have been investigated. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction Huge amount of kinematic energy can be harvested from wind and current [1]. In the industry, the utilisation of the renewable wind and current energy is mainly based on the horizontal axis and vertical axis turbines. However, inspired by the locomotion of birds and ocean cetaceans, the hydrokinetic energy conversion systems based on flapping foil would be a suitable alternative. As discussed in the reviews of Young et al. [2] and Xiao and Zhu [3], the flapping foil is capable of operating in current with lower velocity. The flapping foil turbine would find application to those areas of lower current speed. Following the work of McKinney and DeLaurier [4], a few novel prototypes based on oscillatory foil have been built ([5,6], Pulse Tidal Ltd. http://www.pulsetidal.com). However, the single foil turbine is less efficient than the horizontal and vertical axis turbines. To improve the efficiency of the flapping foil turbine, tandem-foil system will be a competing alternative. Progresses have been made on the energy extraction of single flapping foil. The most popular kinematic mode is the combined sinusoidal heave and pith motions. Parameters, including Strouhal number St, heave amplitude h0 , amplitude of rotational angle 𝜃 0 or the nominal amplitude of the effective attack angle 𝛼 0 and phase difference 𝜀 of the heave and pitch motions, determine the trajectory of the foil [7,8]. The effects of these parameters on the energy extraction efficiency have been presented by Jones and Platzer [9], Jones et al. [10], Dumas and Kinsey [11], Zhu [12], Simpson et al. [13] and Ashraf et al. [14]. Inspired by the locomotion of flying birds and swimming fish, non-sinusoidal motion has been introduced to prescribe the trajectory of the flapping foil [14–17]. The prescribed non-sinusoidal motions, which ∗
Corresponding author. E-mail address: [email protected] (G.D. Xu).
http://dx.doi.org/10.1016/j.enganabound.2017.06.010 Received 1 November 2016; Received in revised form 18 June 2017; Accepted 18 June 2017 0955-7997/© 2017 Elsevier Ltd. All rights reserved.
actually increase the mean effective attack angle during the upstroke or down-stroke motion, can generally improve the efficiency. Liu et al. [18] further simulated the energy harvesting of a deformable foil. Comparing with the rigid foil, the superimposed deformations of the leading and trailing edges enhance the energy extraction efficiency significantly when proper magnitudes of the deformations are adopted. The propulsive efficiency of flying birds in V formation would be improved due to the vortex interaction. The interactions of the downstream wings and the vortices might enhance the thrust or reduce the energy consumption [19]. In the context of energy harvesting of multifoils, hind foil will perform its oscillatory motion in the vortex wake of the fore foil. The average flow velocity around the hind foil would be reduced as the fore foil extracts the kinematic energy from the current. However, foils in tandem configuration might benefit from the interactions of dynamic vortex and the foil when proper position and motions are adopted. The tandem-foil system would harness more power and result higher efficiency. Platzer et al. [6], Ashraf et al. [14] and Kinsey and Dumas [20] have studied the two foils with tandem configuration numerically. As discussed in Kinsey and Dumas [20], the efficiency is sensitive to the phase difference of the motions and the distance of the upstream and downstream foils; in their simulation the maximum efficiency of the tandem foils is as high as 0.64 at a favourable global phase shift 𝜓. Jones et al. [10] carried out experimental studies on tandem foils; the efficiency of the tandem foils is much lower than the numerical results. The difference could be due to the loss of mechanical system. The global phase shift, which indicates the relative position of the downstream foil and the vortex of the upstream foil, can be introduced to study the coupling of the vortex motion and the foils. The corresponding vortex–foil interaction modes with respect to the phase shift of the
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Nomenclature 𝜌 C d U lc
L 𝜔 T = 2𝜋/𝜔 hi (t) h0 𝜃 i (t)
𝜃0 𝜀 𝜎 𝛼𝑖 (𝑡) = 𝜃𝑖 (𝑡) − arctan
ℎ̇ 𝑖 (𝑡) 𝑈
𝛼0
𝑃𝑀𝑖 =
1 𝑇 1 𝑇
of flapping foils with small to moderate effective attack angle at higher Reynolds number, i.e. Re > 1.0 × 105 . The results obtained from inviscid theory had been verified in the work of Anderson et al. [7], La Mantia and Dabnichki [26] and Triantafyllou et al. [27]. Present study adopts the unsteady panel method; only trailing edge vortex shedding is considered. The vortex interactions of the tandem foils have been simulated extensively. The theoretical power coefficients and efficiency of the foils are calculated based on the hydrodynamic forces. For prescribed typical motion parameters, it is found that the efficiency of the hind foil fluctuates against the global phase shift 𝜓 = 2𝜋 𝑈𝐿𝑇 + 𝜎. By increasing 𝜎 or L, the lowest efficiency 𝜂 = 0.05 is found at 𝜓/2𝜋 = 0.05; while the highest efficiency 𝜂 = 0.24 is found at 𝜓/2𝜋 = 0.45. The corresponding vortex interaction modes are investigated. 2. Mathematical formulas 2.1. The motions of the tandem foils
global phase shift horizontal and vertical forces and moment of each foil 𝑓1𝑖 = 2𝐹1𝑖 2𝐹 4𝐹 , 𝑓3𝑖 = 𝜌𝐶 𝑈3𝑖2 , 𝑓5𝑖 = 𝜌𝐶 2 5𝑈𝑖 2 𝜌𝐶 𝑈 2
∫ 𝐹3𝑖 (𝜏)ℎ̇ 𝑖 (𝜏)𝑑𝜏
power due to vertical motion
ℎ1 (𝑡) = ℎ0 sin(𝜔𝑡 + 𝜎), ℎ2 (𝑡) = ℎ0 sin 𝜔𝑡,
(1)
∫ 𝐹5𝑖 (𝜏)𝜃̇ 𝑖 (𝜏)𝑑𝜏
power due to rotational motion
𝜃1 (𝑡) = 𝜃0 sin(𝜔𝑡 + 𝜎 + 𝜀), 𝜃2 (𝑡) = 𝜃0 sin(𝜔𝑡 + 𝜀),
(2)
power due to vertical and rotational motions power coefficients due to vertical motion, rotational motion power consumption due to rotational motion power coefficient
where 𝜔 and 𝜀 are the circular frequency and phase difference respectively, the subscripts 1, 2 indicate the number of the foil. The hind foil will perform oscillatory motion in the vortex wake of the fore foil. There exists a global phase shift of motions between the fore and hind foils in terms of L and 𝜎 [20]. We have 𝐿 𝜓 = 2𝜋 +𝜎 (3) 𝑈𝑇
+𝜎
𝑡+𝑇
𝑡 𝑡+𝑇 𝑡
Pi = PLi + PMi 𝑐𝑃 𝐿𝑖 =
Fig. 1. Sketch of the tandem flapping foils in current.
Fig. 1 shows the tandem configuration of two foils in unbounded flow. The chord of the foil is denoted as C. The Cartesian coordinate system oxz locates at the mean position of the fore foil with x in the direction against the current and z pointing upwards. The hind foil locates after the fore foil. These two foils are denoted as foil 1 and foil 2, respectively. The motions of the foils in the uniform current are prescribed through sinusoidal heave and pitch motions. We have
𝜔ℎ0 𝜋𝑈 𝜓 = 2𝜋 𝑈𝐿𝑇 F1i ,F3i ,F5i
𝑆𝑡 =
𝑃𝐿𝑖 =
density of the fluid chord of the foil the maximum vertical extension of the foil velocity of the uniform current the distance from leading edge to rotational centre (xci , zci ) and lc /C = 1/3, i = 1, 2 indicates the fore and hind foils longitudinal distance of the two foils at their rotational centre circular frequency of oscillatory motion period of the oscillatory motion heave motion of the foil the heave amplitude pitch motion, the positive direction is defined as counterclockwise the pitch amplitude the phase difference of the heave and pitch motions the phase shift of the motions of the two foils effective attack angle at the rotational centre, dot over hi (t) means the derivative against time nominal amplitude of effective attack angle Strouhal number
𝑃𝐿𝑖 ; 𝑐𝑃 𝑀𝑖 0.5𝜌𝐶 𝑈 3
𝑐𝑃𝑖𝑛𝑀𝑖 = −𝑐𝑃 𝑀𝑖 𝑃𝑖 0.5𝜌𝐶 𝑈 3 𝑃𝑖 0.5𝜌𝑈 3 𝑑
𝑐𝑃 𝑖 =
𝜂𝑖 = 𝜂t = 𝜂1 + 𝜂2
=
𝑃𝑀𝑖 0.5𝜌𝐶 𝑈 3
power efficiency the total efficiency of the tandem foils cPL , cPM , cP and 𝜂 are the corresponding coefficient for single foil
Substituting 𝑆𝑡 = 𝜓 = 𝑆𝑡
𝜔ℎ0 𝜋𝑈
and 𝑇 =
2𝜋 𝜔
into Eq. (3), we have
𝐶 𝐿 +𝜎 ℎ0 2𝐶
(4)
The global phase shift is the function of Strouhal number St, heave amplitude h0 , longitudinal distance L and phase shift 𝜎. The vortex interactions of the foils depend on these parameters. We note that the vortex structure of the flapping foil is determined when the motion parameters St, h0 /C, 𝛼 0 and 𝜀 are prescribed. Different motion parameters will result in different vortex structure. For prescribed motions, the global phase shift, which unifies the L and 𝜎, will affect the efficiency of the hind foil significantly.
motion 𝜎 and the longitudinal distance L are less understood. The mechanism of the vortex interaction and the corresponding global phase shift require systematic investigations. The energy harvesting of the tandem foils in unbounded flow with typical oscillatory motions will be analysed in present analysis. For oscillatory foils, leading edge vortex shedding would be reduced or even eliminated even though larger attack angle is considered [21,22]. The leading edge vortex shedding would be delayed at higher Reynolds number and would result in the increase of energy extraction efficiency [10,23–25]. These justify the unsteady panel method for the simulation
2.2. Governing equations and the boundary value problem We consider the tandem foils in a uniform current with periodical heave and pitch motions. By assuming that the Reynolds number is high 203
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
enough and the effective attack angle ranges from small to medium, the boundary layer is thin and the vortex only sheds at the trailing edge. The inviscid flow around the foils is described by the velocity potential flow theory. The total velocity potential Φ in the presence of the current satisfies Laplace equation. We have ∇2 Φ = 0
The Kelvin theorem requires that the circulation of the fluid domain should be constant. The total circulation of the vortices and foils shall be zero. The summation of the shed vortex and the circulation over the foils shall be equal but opposite in sign. We note that the strength of each point vortex calculated through Eq. (16) will not change during the evolution of the vortex. Eqs. (14) and (16) guarantee the circulation of the fluid domain be zero and the Kelvin theorem be satisfied. The position of the point vortex can be updated through following equations
(5)
The boundary condition on the foil surface S0 can be written as ⇀ ⇀ 𝜕Φ = (−𝜃̇ 𝑖 𝑍, ℎ̇ 𝑖 + 𝜃̇ 𝑖 𝑋 𝑖 ) 𝑛 = −𝜃̇ 𝑖 𝑍𝑖 𝑛𝑥 + (ℎ̇ 𝑖 + 𝜃̇ 𝑖 𝑋𝑖 )𝑛𝑧 𝜕𝑛
(6)
⇀
where 𝑛 = (𝑛𝑥 , 𝑛𝑧 ) is the inward normal of the foil surface, 𝜃̇ 𝑖 and ℎ̇ 𝑖 are the rotational velocity and the vertical velocity respectively, and ⇀
𝑋 𝑖 = (𝑋𝑖 , 𝑍𝑖 ) = (𝑥 − 𝑥𝑐𝑖 , 𝑧 − 𝑧𝑐𝑖 ), (xci ,zci ) denotes the rotational centre of the corresponding fore or hind foil. In the far field ∇Φ = −Ui, 𝑥∕𝐶 >> 0 or |𝑧|∕𝐶 >> 0
where i indicates the horizontal unit vector. The velocity potential 𝜙 due to the disturbance the foil motion is introduced. We have (8)
Substituting Eq. (8) into Eqs. (5)–(7), we obtain the governing equations and boundary conditions ∇2 𝜙 = 0 in the f luid domain
(17)
𝑧𝜅(𝑗 ) (𝑡 + Δ𝑡) = 𝑧𝜅(𝑗 ) (𝑡) + 𝑤𝜅(𝑗 ) 𝑑𝑡
(18)
where u𝜅(j) and w𝜅(j) are the local horizontal and vertical velocities, which will be referred in the appendix. As the shed point vortex of the fore foil approaches the hind foil, it could penetrate the foil surface due to numerical error. A fence scheme is adopted to avoid the penetration [28]. The fence is an artificial surface that has an offset 𝛿 to the foil surface. When the distance of the point vortex to the foil surface is smaller than 𝛿, the point vortex will be placed on the nearest artificial fence. The offset 𝛿 over the foil surface is approximately the size of the element on the foil. To solve the boundary value problem, we have the integral equation based on Green’s third identity
(7)
Φ(𝑥, 𝑧) = 𝜙(𝑥, 𝑧) − 𝑈 𝑥 or 𝜙(𝑥, 𝑧) = Φ(𝑥, 𝑧) + 𝑈 𝑥
𝑥𝜅(𝑗 ) (𝑡 + Δ𝑡) = 𝑥𝜅(𝑗 ) (𝑡) + (−𝑈 + 𝑢𝜅(𝑗 ) )𝑑𝑡
(9)
Λ(𝑝)𝜙(𝑝) =
𝑆01 +𝑆02
𝜕𝜙 𝜕Φ = + 𝑈 𝑖 = (𝑈 − 𝜃̇ 𝑖 𝑍𝑖 )𝑛𝑥 + (ℎ̇ 𝑖 + 𝜃̇ 𝑖 𝑋𝑖 )𝑛𝑧 on foil surface 𝜕𝑛 𝜕𝑛
(10)
∇𝜙 = ∇Φ + 𝑈 𝑖 = 0 in the far f ield
(11)
+
∫
𝜕 𝐺(𝑝, 𝑞 ) 𝜕𝜙(𝑞) 𝜙(𝑞) − 𝐺(𝑝, 𝑞) ]𝑑𝑆 𝜕 𝑛𝑞 𝜕 𝑛𝑞
𝑆 𝑤 1 +𝑆 𝑤 2
The Kutta condition requires that the velocity is finite and continuous at the trailing edge. We have [8] 𝜕 𝜙+ || 𝜕 𝜙− || 𝜕𝜇 || − |⇀ = 𝜕𝑠 |⇀ 𝜕𝑠 ||⇀ 𝜕𝑠 | |𝑥𝑇 𝑥𝑇 𝑥𝑇
[
∫
+
𝑚 ∑ ( 𝑗=1
⇀ ⇀ ⇀ ⇀ 𝜕 𝐺(𝑝, 𝑞 ) 𝜇(𝑞) 𝑑𝑆 +𝜇(𝑥′𝑇 1 )𝐻(𝑝, 𝑥′𝑇 1 )+𝜇(𝑥′𝑇 2 )𝐻(𝑝, 𝑥′𝑇 2 ) 𝜕 𝑛𝑤
𝐻 (𝑝, 𝑞)𝜅𝑗1 + 𝐻 (𝑝, 𝑞)𝜅𝑗2
)
(19)
where Λ(p) is the solid angle at point p(x, z), G(p, q) = ln r, 𝑟 = √ (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2 , and q(𝜉, 𝜂) is the point on the integral boundary,
(12)
where 𝜙+ and 𝜙- are the potential on the upper and lower surface at the − trailing edge of the foil, 𝜇 = 𝜙+ 𝑤 − 𝜙𝑤 is the difference of the potential
𝐻(𝑝, 𝑞) = arctan
𝑥−𝜉 𝑧−𝜂
is obtained through the integral by parts.
⇀
on the upper and lower surface of the wake, 𝑥𝑇 = (𝑥𝑇 , 𝑧𝑇 ) denotes the 𝜕 trailing edge and 𝜕𝑠 denotes the partial derivative along the foil surface or vortex wake. The pressure p across the wake shall be continuous or p + = p − . In the form of Lagrangian, we have 𝑑𝜇 =0 𝑑𝑡 The above equation can be written in its discretised form as ⇀
⇀
𝜇(𝑥′𝑇 , 𝑡 + Δ𝑡) = 𝜇(𝑥𝑇 , 𝑡)
2.3. Calculation of the hydrodynamic forces Using Bernoulli equation, we have the pressure on each foil surface [ 𝑝𝑖 = −𝜌
(13)
(14)
where 𝑥′𝑇 = (𝑥′𝑇 , 𝑧′𝑇 ) is the end point of the dipole element connecting
𝐹𝑗𝑖 =
⇀
to the trailing edge and the strength of dipole 𝜇(𝑥𝑇 , 0) = 0 at first time step. In the time stepping scheme, the hind foil would penetrate the continuous vortex wake of the fore foil. To avoid the breaking of vortex dipole, the continuous vortex dipole is reduced into point vortex. Following the time stepping procedures in Xu and Wu [8], the dipole element can be converted into equivalent point vortex and a new dipole connecting to the trailing edge will develop at j + 1 time step. We define the vortex strength 𝛾 = 𝜕𝜇 . From Eq. (13) we have 𝜕𝑠
⇀
∫
𝑝𝑖 𝑛𝑗 𝑑𝑆 , 𝑗 = 1, 3, 5,
(21)
𝑠0
where (n1 ,n3 ,n5 ) = (nx ,nz ,−Xnz + Znx ). 3. Results and discussions 3.1. Typical vortex structure of single foil The vortex structure of a NACA 0012 foil performing sinusoidal heave and pitch motions is √ analysed. The uniform current velocity U = 4.43 and C = 2.0 or 𝑈 ∕ 𝑔𝐶 = 1.0. When solving the boundary integral equation, the linear distribution of the potential on each panel and the wake dipole is assumed. With finer element near the leading and trailing edges, numerical test shows that N = 60 elements over the foil surface can result convergent solutions. The time step interval Δt should make sure that the motion of point vortex at each time step is
𝑑𝛾 =0 (15) 𝑑𝑡 The strength of newly developed point vortex at time step j can be obtained in its discretised form ⇀
(20)
𝑑 where 𝑑𝑡 is the derivative fixed on the foil. Consequently, the force and moment can be calculated through the integration of pressure, we have
⇀
𝜅𝑗 (𝑡 + Δ𝑡) = 𝜇(𝑥𝑇 , 𝑡) − 𝜇(𝑥′𝑇 , 𝑡)
] 𝑑𝜙 1 − (𝑈 − 𝜃̇ 𝑖 𝑍𝑖 , ℎ̇ 𝑖 + 𝜃̇ 𝑖 𝑋𝑖 ) ⋅ ∇𝜙 + ∇𝜙∇𝜙 𝑑𝑡 2
(16)
204
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Fig. 2. (a) The typical vortex structure of a single foil and (b) the induced velocity by the vortices at t/T = 4.0 with St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, 𝜀 = 90°.
⇀
less than the average size of the mesh on the foil, or Δ𝑡 < 2𝐶∕(𝑁 𝑈 ),
Fig. 3. (a) The time histories of the vertical forces and (b) the comparison with the results of Kinsey and Dumas [20] which is referred as K&D in the figure, tandem NACA 0015 foils with L/C = 5.4, 𝜎 = 𝜋, St = 0.28, h0 /C = 1.0 and 𝜃 0 = 70° are simulated.
⇀
where 𝑈 = |𝑈 + 𝑢, 𝑤| is the velocity of the local flow. As shown in Fig. 2(a), the vortex sheds from the trailing edge successively. The structure of these vortices is a typical reversed Karman vortex street. The arrows in Fig. 2(b) indicate the vortex induced flow as the uniform current speed is excluded. When the hind foil moves in the vortex street, the average local fluid velocity will decrease and the vortex flow will change its effective attack angle significantly. The vortex structure in Fig. 2 suggests that the phase shift 𝜎 and the distance L would determine the relative position of the hind foil and the vortices of the fore foil. In the following section the tandem NACA 0012 foils with same motion parameters in Fig. 2 will be investigated.
3.3. Effects of phase shift 𝜎 The phase shift 𝜎 is the phase difference of the heave motions; it determines the relative vertical position of the hind foil and the vortices of the fore foil. We fix the motion of the two foils at St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, and 𝜀 = 90°. Simulations are carried out with L/UT = 0.5 and 𝜎 ranging from − 𝜋 to 𝜋. The energy efficiencies of the two foils at various 𝜎 are shown in Fig. 4(a). The horizontal axis on top of the figure is the corresponding global phase shift, which is calculated through Eq. (3). As shown in Fig. 4(a), there is no significant variation of the efficiency of the fore foil 𝜂 1 ; while the efficiency of the hind foil 𝜂 2 experiences one oscillation as 𝜎 increases from − 𝜋 to 𝜋. The efficiency of the hind foil is a periodical function of the phase shift 𝜎. The lowest efficiency of hind foil is found at 𝜎/2𝜋 = −0.45 or 𝜓/2𝜋 = 0.05; the highest efficiency is achieved at 𝜎/2𝜋 = −0.05 or 𝜓/2𝜋 = 0.45. To further investigate the power extraction of the hind foil, the power due to lifting force cPL2 and the consumption power due to the rotational motion 𝑐𝑃𝑖𝑛𝑀2 are presented in Fig. 4(b). Both cPL2 and 𝑐𝑃𝑖𝑛𝑀2 experience one oscillation. Power coefficient due to lifting force cPL2 is generally lower than that of a single foil, as shown in Fig. 4(b). At 𝜎/2𝜋 = −0.05, the curve of cPL2 finds its crest while the curve of 𝑐𝑃𝑖𝑛𝑀2 finds its trough. These result in the highest power coefficient and efficiency. Fig. 5 shows the typical time histories of the hind foil at 𝜎/2𝜋 = −0.5 (𝜓/2𝜋 = 0.0) and 𝜎/2𝜋 = −0.05 (𝜓/2𝜋 = 0.45). The amplitudes of horizontal force f1 and vertical force f3 at 𝜎/2𝜋 = −0.05 are larger than those of 𝜎/2𝜋 = −0.5, as shown in Fig. 5(a) and (b). The moment f5 at 𝜎/2𝜋 = −0.05 is smaller than that of 𝜎/2𝜋 = −0.5. The time histories of the forces of a single foil performing the same motion are presented in Fig. 5. Although the force amplitudes are larger than those of the tandem foils and therefore result in higher power coefficient cPL , the energy consumption due to rotational moment 𝑐𝑃𝑖𝑛𝑀 compensates the extracted power cPL . Smaller moment amplitude means the reduction of the consumed power. Consequently, the efficiency of hind foil at 𝜎/2𝜋 = −0.05 finds its peak.
3.2. Validation When tandem foils are considered, the hind foil will perform its motion in the vortex wake of the fore foil. To validate the numerical code, the interaction of tandem NACA 0015 foils with L/C = 5.4, 𝜀 = 𝜋, St = 0.28, h0 /C = 1.0 and 𝜃 0 = 70° is simulated firstly. Convergence study is carried out before extensive simulations. We used 60, 90 and 140 elements on the foil surface; the discrepancies of the forces are negligible. In the following simulations, 140 elements on each foil are adopted. The simulations stop after eight oscillations. The time history of vertical force of first five periods is shown in Fig. 3(a). The force of the hind foil is not affected in the first half oscillation. The effects of the vortex interaction on the hind foil become stable after the first oscillation. In the following figures, only one stable period will be shown. To compare the hydrodynamic forces of the fore and hind foils, the relative time t′ = t − 𝜎 for the fore foil is introduced. As shown in Fig. 3(b), the curves of present simulation are generally in good agreement with those of Kinsey and Dumas [20], whose results were obtained through solving two-dimensional Reynolds average Navier–Stokes equation at Reynolds number 5.0 × 105 . The snapshot of the flow field in Kinsey and Dumas [20] showed that there was no flow separation at the leading edge of the foil. The boundary layer was thin comparing with the foil. It can be treated as the high Reynolds number problem. The amplitude of present curves is slightly larger since the viscosity of the fluid has been ignored. 205
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Fig. 4. (a) The efficiency of the tandem foils and (b) the cPL2 and cPM2 of the hind foil at various 𝜎 at L/UT = 0.5, with St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, 𝜀 = 90° the dash line and dash dot line present cPL and cPM of single foil, respectively.
Fig. 5. The time history of (a)f12 , (b)f32 , (c)f52 of the hind foil at 𝜎/2𝜋 = −0.5 and 𝜎/2𝜋 = −0.05.
Typical vortex flow of the tandem foils is further investigated. Fig. 6 shows the flow velocity induced by the vortices at 𝜎/2𝜋 = −0.5. The vortices of the fore foil coincide with the vortices of the hind foil; stronger vortices are developed. The pattern of the vortex is similar with that of a single foil. When the hind foil performs the upstroke motion, it encounters the positive vortices of the fore foil; the hind foil encounters the negative vortices during its downstroke motion. The local flow velocity has been reduced and the effective attack angle decreases. Therefore the extracted power decreases significantly. The vortex pattern in Fig. 6 is named as synchronised mode. Fig. 7 gives the vortex flow at 𝜎/2𝜋 = −0.05. The hind foil encounters the negative vortices in the upstroke motion and it encounters positive vortices in the downstroke motion. The local flow velocity of the hind foil decreases due to the vortices. However, the effective attack angle increases. Therefore a rel-
atively higher efficiency has been achieved. The vortex of the tandem foils develops into two columns of vortices or vortex pairs, with positive vortices above the mean position and negative vortices below the mean position. The vortex pattern of the tandem foils is named as the interlaced mode. 3.4. The effects of longitudinal distance L In terms of the global phase shift, the longitudinal distance L has equivalent effects as 𝜎, as indicated in Eq. (4). The distance L determines the relative longitudinal position of the hind foil and the vortex of the fore foil. Its effects on the performance of the tandem foil are investigated. Again, the parameters of the motion are fixed at St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, and 𝜀 = 90°. Simulations are carried out 206
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Fig. 6. The flow induced by the vortices at 𝜓/2𝜋 = 0, (a) t/T = 0, (b) t/T = 1/8, (c) t/T = 1/4, (d) t/T = 3/8.
Fig. 7. Flow induced by the vortices at 𝜓/2𝜋 = 0.45, (a) t/T = 0, (b) t/T = 1/8, (c) t/T = 1/4, (d) t/T = 3/8.
with 𝜎 = −𝜋 and L/UT ranging from 0.5 to 1.0. The efficiency of the two foils at various L is shown in Fig. 8. The efficiency of the fore foil is not affected by the distance L obviously; while the efficiency of the hind foil fluctuates with L. The peak of the efficiency of the hind foil is as high as that of the fore foil at L/UT = 0.9 (𝜓/2𝜋 = 0.4); the lowest efficiency is found at L/UT = 0.55 (𝜓/2𝜋 = 0.05). Comparing Fig. 8(a) and Fig. 4(a), the curves in these two figures are resemble. The highest and lowest efficiencies are found with similar global phase shift 𝜓. This suggests that the effects of distance L have similar effects as the phase shift 𝜎. The curves of cPL2 and 𝑐 𝑖𝑛 are further shown in Fig. 8(b). These 𝑃 𝑀2 two waving curves of cPL2 and 𝑐 𝑖𝑛 find the peak and trough at L/UT 𝑃 𝑀2 ≈ 0.9 respectively, which result in the highest efficiency of the hind foil. The typical time histories of the forces of hind foil at L/UT = 0.5 and L/UT = 0.95 are further presented in Fig. 9. The horizontal and vertical forces of L/UT = 0.95 in Fig. 9(a) and (b) have larger amplitudes than those of L/UT = 0.5; the amplitude of moment at L/UT = 0.95 is smaller
than that of L/UT = 0.5. As a result, the efficiency with L/UT = 0.95 is higher than that of L/UT = 0.5. Comparing with the curves of a single foil, the amplitudes of the forces of tandem foils are smaller. However, higher efficiency is achieved due to the decrease of the energy consumption of the pitch moment at L/UT = 0.95. It is expected that the vortex interactions would be the same as those in Figs. 6 and 7. The synchronised mode and interlaced mode will form at L/UT = 0.5 (𝜓/2𝜋 = 0) and L/UT = 0.95 (𝜓/2𝜋 = 0.45) respectively. 4. Concluding remarks The energy harvesting of flapping foils with tandem configurations in uniform current has been simulated through the boundary elements. The hind foil is moving in the vortex wake of the fore foil. The effects of vortex–foil interaction on the energy extraction performance are analysed. There are strong interactions of foils and vortices. The efficiency of the fore foil is not affected significantly; however the power 207
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Fig. 8. (a) The efficiency of the tandem foils and (b) the cPL2 and cPM2 of the hind foil against L at 𝜎 = −𝜋, with St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, 𝜀 = 90°.
Fig. 9. The time history of (a)f12 , (b)f32 , (c)f52 of the hind foil at L/UT = 0.5 and L/UT = 0.95 in Fig. 8.
efficiency of the hind foil fluctuates with the global phase shift 𝜓. The lowest efficiency of the hind foil is around 𝜂 2 = 0.05 and its highest efficiency is 𝜂 2 = 0.24. The efficiency curves against L and 𝜎 resemble each other. It is interesting to find that the peaks of efficiency are found at global phase 𝜓/2𝜋 ≈ 0.45, where the “interlaced mode” can be observed. The lowest efficiency appears when 𝜓/2𝜋 ≈ 0 or 𝜓/2𝜋 ≈ 1.0, where the “synchronised mode” is found. The global phase shift 𝜓 indicates the vortex interaction modes and 𝜎 and L have equivalent effects on the efficiency. Present simulations are limited to one set of typical motion parameters. The effects of the motion parameters, including St, h0 , 𝛼 0 and 𝜀, have been ignored. We noted that the change of these parameters will cause the variation of the vortex structure. Even though the same global phase shift is adopted, it will not necessarily cause the increase or decrease of the efficiency due to the complexity of the foil–vortex interactions. For instance, when small Strouhal number is adopted, the longitudinal distance of the vortices will increase and the strength of
the vortices will decrease; the variation of the effective attack angle will become less significant. When large Strouhal number is concerned, the distance of the vortices will decrease. The average effective attack angle will not increase or decrease significantly if the vertical and longitudinal distances of the vortices are smaller than the size of the foil. The effects of different motion parameters, including St, h0 , 𝛼 0 and 𝜀, require further investigations. Present study would be highly important for the design and construction of hydro-kinematic turbines based on flapping tandem foils. A proper distance and phase shift of the motions can be selected to assure the efficiency of the tandem foil system. For fixed longitudinal distance, the foils can adapt the phase shift of the motions to obtain the most favourable global phase shift. When the phase shift of the motions is fixed, we could adjust the longitudinal distance. The tandem configuration of the foils can be implemented through two crankshaft structures like those in Xu et al. [29], with the shaft of two sets of foil-crankshaft mechanisms being linked through chain wheels.
208
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Acknowledgement
Reference
The research is under the support of National Science Foundation of China (Grant nos. 11602067 and 51479044) and State Key laboratory of Ocean Engineering (Shanghai Jiao Tong University, Grant No. 1413).
[1] Schiermeier Q, Tollefson J, Scully T, Witze A, Morton O. Electricity without carbon. Nature 2008;454(7206):816–23. [2] Young J, Lai J, Platzer M. 2014 A review of progress and challenges in flapping foil power generation Prog Aerosp Sci 67, 2–28. [3] Xiao Q, Zhu Q. A review on flow energy harvesters based on flapping foils. J Fluids Struct 2014;46:174–91. [4] McKinney W, DeLaurier J. Wingmill: an oscillating-wing windmill. J Energy 1981;5(2):109–15. [5] Kinsey T, Dumas G, Lalande G, Ruel J, Mehut A, Viarouge P, et al. Prototype testing of a hydrokinetic turbine based on oscillating hydrofoils. Renew Energy 2011;36:1710–18. [6] Platzer M, Ashraf M, Young J, Lai J. Development of a new oscillating-wing wind and hydropower generator. In: Proceedings of the forty-seventh AIAA aerospace sciences meeting. Reston, Virginia: American Institute of Aeronautics and Astronautics; 2009. [7] Anderson JM, Streitlien K, Barrett DS, Triantafyllou MS. Oscillating foils of high propulsive efficiency. J Fluid Mech 1998;360:41–72. [8] Xu GD, Wu GX. Boundary element simulation of inviscid flow around an oscillatory foil with vortex sheet. Eng Anal Bound Elem 2013;37:825–35. [9] Jones KD, Platzer MF. Numerical computation of flapping-wing propulsion and power extraction. AIAA J 1997:97–0826 35th AIAA Aerospace Sciences Meeting, Reno, Nevada, Jan. [10] Jones KD, Lindsey K, Platzer MF. An investigation of the fluid–structure interaction in an oscillating-wing micro-hydropower generator. Fluid–structure interaction II. Southampton, UK: Wessex Institute of Technology Press; 2003. p. 73–82. [11] Dumas G, Kinsey T. Eulerian simulations of oscillating airfoils in power extraction regime. In: Rahman B, editor. Proceedings in advances in fluid mechanics VI. WIT Press; 2006. p. 245–54. [12] Zhu Q. Optimal frequency for flow energy harvesting of a flapping foil. J Fluid Mech 2011;675:495–517. [13] Simpson BJ, Licht S, Hover FS, Triantafyllou MS. Energy extraction through flapping foils. In: Proceedings of the twenty-seventh international conference on offshore mechanics and arctic engineering. OMAE; 2008. p. 2008–58043. [14] Ashraf MA, Young J, Lai JCS, Platzer MF. Numerical analysis of an oscillating-wing wind and hydropower generator. AIAA J 2011;49(7):1374–86. [15] Lu K, Xie Y, Zhang D. Nonsinusoidal motion effects on energy extraction performance of a flapping foil. Renew Energy 2014;64:283–93. [16] Platzer M, Ashraf M, Young J, Lai J. Extracting power in jet streams: pushing the performance of flapping wing technology. In: Proceedings of the twenty-seventh congress of the international council of the aeronautical sciences. Nice France: International Council of the Aeronautical Sciences, September 19–24; 2010 Paper 2010-2.9.1. [17] Xiao Q, Liao W, Yang SC, Peng Y. How motion trajectory affects energy extraction performance of a biomimic energy generator with an oscillating foil? Renew Energy 2012;37:61–75. [18] Liu WD, Xiao Q, Cheng F. A bio-inspired study on tidal energy extraction with flexible flapping wings. Bioinspiration Biomim 2013;8:1–16. [19] Lehmann F-O. Wing–wake interaction reduces power consumption in insect tandem wings. Exp Fluids 2009;46(5):765–75. [20] Kinsey T, Dumas G. Optimal tandem configuration for oscillating-foils hydrokinetic turbine. J Fluids Eng 2012;134(3):031103-1–031103-11. [21] Dickinson MH, Lehmann FO, Sane SP. Wing rotation and the aerodynamic basis insect flight. Science 1999;284:1954–60. [22] Ellington CP, Vanderberg C, Wilmott A, Thomas A. Leading edge vortices in insect flight. Nature 1996;384:626–30. [23] Campobasso MS, Piskopakis A, Drofelnik J, Jackson A. Turbulent Navier–Stokes analysis of an oscillating wing in a power-extraction regime using the shear stress transport turbulence model. Comput Fluids 2013;88:136–55. [24] Kinsey T, Dumas G. Parametric study of an oscillating airfoil in a power-extraction regime. AIAA J 2008;46:1318–30. [25] Kinsey T, Dumas G. Optimal operating parameters for an oscillating foil turbine at Reynolds number 500,000. AIAA J 2014;52(9):1885–95. [26] La Mantia M, Dabnichki P. Unsteady panel method for flapping foil. Eng Anal Bound Elem 2009;33:572–80. [27] Triantafyllou MS, Hover FS, Techet AH, Yue DKP. Review of hydrodynamic scaling laws in aquatic locomotion and fishlike swimming. J Am Soc Mech Eng 2005;58:226–37. [28] Pan YL, Dong XX, Zhu Q, Yue DKP. Boundary-element method for the prediction of performance of flapping foils with leading-edge separation. J Fluid Mech 2012;698:446–67. [29] Xu GD, Xu WH, Dai J. Numerical and experimental study of a flapping foil generator. Appl Ocean Res 2017;63:242–50. [30] Li S, Han R, Zhang AM, Wang QX. Analysis of pressure field generated by a collapsing bubble. Ocean Eng 2016;117:22–38.
Appendix Once the boundary value problem has been solved, the potential and normal velocity on the integral boundary are known. However the velocity of a point in the fluid domain cannot be obtained directly [30]. Here we shall derive a numerical scheme to calculate the local velocity. Combine the similar terms of the two foils in Eq. (19), we have [ Λ(𝑝)𝜙(𝑝) =
𝜕 ln 𝑟𝑝𝑞
∫
𝜕 𝑛𝑞
𝑆0 ⇀
𝜙(𝑞) − ln 𝑟𝑝𝑞
+𝜇(𝑥′𝑇 ) arctan
] 𝜕 ln 𝑟𝑝𝑞 𝜕𝜙(𝑞) 𝑑𝑆 + 𝜇(𝑞)𝑑𝑆 ∫ 𝜕 𝑛𝑤 𝜕 𝑛𝑞 𝑆𝑤
𝑥 − 𝑥′𝑇 𝑧 − 𝑧′𝑇
( ) 𝑥 − 𝑥𝑗 + arctan 𝜅𝑗 𝑧 − 𝑧𝑗 𝑗=1 𝑚 ∑
(22)
where S0 , Sw and (xj , zj ) denote the foil surface, wake panel and position of the point vortex, respectively. The solid angle Λ(p) = 2𝜋 if p is within the fluid domain. Applying the integral by part to the first term in the brackets and the second integral, we have 2𝜋𝜙(𝑝) = −
𝜕 𝜙(𝑞 ) 𝑥−𝜉 𝜕 𝜙(𝑞 ) arctan ln 𝑟𝑝𝑞 𝑑𝑆 𝑑𝑆 − ∫ 𝜕𝑠 𝑧−𝜂 𝜕𝑛
∫
𝑆0
+
𝑆0
) 𝑚 ( ∑ 𝑥 − 𝑥𝑗 𝜕 𝜇(𝑞 ) 𝑥−𝜉 𝑑𝑆 + arctan arctan 𝜅𝑗 ∫ 𝜕𝑠 𝑧−𝜂 𝑧 − 𝑧𝑗 𝑗=1
(23)
𝑆𝑤
Applying the derivative against x and z, we have 2𝜋
𝜕𝜙(𝑝) 𝜕 𝜙(𝑞 ) 𝑧−𝜂 =− 𝑑𝑆 ∫ 𝜕𝑥 𝜕𝑠 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2 𝑆0
−
∫
𝜕 𝜙(𝑞 ) 𝑥−𝜉 𝑑𝑆 𝜕𝑛 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2
∫
𝜕 𝜇(𝑞 ) 𝑧−𝜂 𝑑𝑆 𝜕𝑠 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2
𝑆0
+
𝑆𝑤
+
𝑚 ∑
(
(𝑥 − 𝑥𝑗 )2 + (𝑧 − 𝑧𝑗 )2
𝑗=1
2𝜋
𝑧 − 𝑧𝑗
) 𝜅𝑗
(24)
𝜕𝜙(𝑝) 𝜕 𝜙(𝑞 ) 𝑥−𝜉 = 𝑑𝑆 ∫ 𝜕𝑧 𝜕𝑠 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2 𝑆0
−
∫
𝜕 𝜙(𝑞 ) 𝑧−𝜂 𝑑𝑆 𝜕𝑛 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2
∫
𝜕 𝜇(𝑞 ) 𝑥−𝜉 𝑑𝑆 𝜕𝑠 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2
𝑆0
+
𝑆𝑤
−
𝑚 ∑ 𝑗=1
(
𝑥 − 𝑥𝑗 (𝑥 − 𝑥𝑗 )2 + (𝑧 − 𝑧𝑗 )2
) 𝜅𝑗
(25)
The velocity at point p in the fluid domain can be obtained through 𝑝) 𝜕𝜙(𝑝) (𝑢, 𝑤) = ( 𝜕𝜙( , 𝜕𝑧 ). 𝜕𝑥 209
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Energy extraction of two flapping foils with tandem configurations and vortex interactions G.D. Xu∗, W.H. Xu College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China
a r t i c l e
i n f o
Keywords: Flapping foils Tandem configurations Energy extraction Global phase shift Efficiency
a b s t r a c t The energy extraction of two flapping foils with tandem configurations has been studied based on velocity potential theory. The interactions of the hind foil and the vortices of the fore foil with typical oscillatory motions are investigated. The global phase shift, which combines the longitudinal distance and the phase shift of the motions of the two foils, is introduced as the indicating parameter of the interaction modes. The effects of longitudinal distance and phase shift of the motions on the energy extraction efficiency have been analysed. The most favourable global phase shift and the corresponding vortex pattern have been investigated. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction Huge amount of kinematic energy can be harvested from wind and current [1]. In the industry, the utilisation of the renewable wind and current energy is mainly based on the horizontal axis and vertical axis turbines. However, inspired by the locomotion of birds and ocean cetaceans, the hydrokinetic energy conversion systems based on flapping foil would be a suitable alternative. As discussed in the reviews of Young et al. [2] and Xiao and Zhu [3], the flapping foil is capable of operating in current with lower velocity. The flapping foil turbine would find application to those areas of lower current speed. Following the work of McKinney and DeLaurier [4], a few novel prototypes based on oscillatory foil have been built ([5,6], Pulse Tidal Ltd. http://www.pulsetidal.com). However, the single foil turbine is less efficient than the horizontal and vertical axis turbines. To improve the efficiency of the flapping foil turbine, tandem-foil system will be a competing alternative. Progresses have been made on the energy extraction of single flapping foil. The most popular kinematic mode is the combined sinusoidal heave and pith motions. Parameters, including Strouhal number St, heave amplitude h0 , amplitude of rotational angle 𝜃 0 or the nominal amplitude of the effective attack angle 𝛼 0 and phase difference 𝜀 of the heave and pitch motions, determine the trajectory of the foil [7,8]. The effects of these parameters on the energy extraction efficiency have been presented by Jones and Platzer [9], Jones et al. [10], Dumas and Kinsey [11], Zhu [12], Simpson et al. [13] and Ashraf et al. [14]. Inspired by the locomotion of flying birds and swimming fish, non-sinusoidal motion has been introduced to prescribe the trajectory of the flapping foil [14–17]. The prescribed non-sinusoidal motions, which ∗
Corresponding author. E-mail address: [email protected] (G.D. Xu).
http://dx.doi.org/10.1016/j.enganabound.2017.06.010 Received 1 November 2016; Received in revised form 18 June 2017; Accepted 18 June 2017 0955-7997/© 2017 Elsevier Ltd. All rights reserved.
actually increase the mean effective attack angle during the upstroke or down-stroke motion, can generally improve the efficiency. Liu et al. [18] further simulated the energy harvesting of a deformable foil. Comparing with the rigid foil, the superimposed deformations of the leading and trailing edges enhance the energy extraction efficiency significantly when proper magnitudes of the deformations are adopted. The propulsive efficiency of flying birds in V formation would be improved due to the vortex interaction. The interactions of the downstream wings and the vortices might enhance the thrust or reduce the energy consumption [19]. In the context of energy harvesting of multifoils, hind foil will perform its oscillatory motion in the vortex wake of the fore foil. The average flow velocity around the hind foil would be reduced as the fore foil extracts the kinematic energy from the current. However, foils in tandem configuration might benefit from the interactions of dynamic vortex and the foil when proper position and motions are adopted. The tandem-foil system would harness more power and result higher efficiency. Platzer et al. [6], Ashraf et al. [14] and Kinsey and Dumas [20] have studied the two foils with tandem configuration numerically. As discussed in Kinsey and Dumas [20], the efficiency is sensitive to the phase difference of the motions and the distance of the upstream and downstream foils; in their simulation the maximum efficiency of the tandem foils is as high as 0.64 at a favourable global phase shift 𝜓. Jones et al. [10] carried out experimental studies on tandem foils; the efficiency of the tandem foils is much lower than the numerical results. The difference could be due to the loss of mechanical system. The global phase shift, which indicates the relative position of the downstream foil and the vortex of the upstream foil, can be introduced to study the coupling of the vortex motion and the foils. The corresponding vortex–foil interaction modes with respect to the phase shift of the
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Nomenclature 𝜌 C d U lc
L 𝜔 T = 2𝜋/𝜔 hi (t) h0 𝜃 i (t)
𝜃0 𝜀 𝜎 𝛼𝑖 (𝑡) = 𝜃𝑖 (𝑡) − arctan
ℎ̇ 𝑖 (𝑡) 𝑈
𝛼0
𝑃𝑀𝑖 =
1 𝑇 1 𝑇
of flapping foils with small to moderate effective attack angle at higher Reynolds number, i.e. Re > 1.0 × 105 . The results obtained from inviscid theory had been verified in the work of Anderson et al. [7], La Mantia and Dabnichki [26] and Triantafyllou et al. [27]. Present study adopts the unsteady panel method; only trailing edge vortex shedding is considered. The vortex interactions of the tandem foils have been simulated extensively. The theoretical power coefficients and efficiency of the foils are calculated based on the hydrodynamic forces. For prescribed typical motion parameters, it is found that the efficiency of the hind foil fluctuates against the global phase shift 𝜓 = 2𝜋 𝑈𝐿𝑇 + 𝜎. By increasing 𝜎 or L, the lowest efficiency 𝜂 = 0.05 is found at 𝜓/2𝜋 = 0.05; while the highest efficiency 𝜂 = 0.24 is found at 𝜓/2𝜋 = 0.45. The corresponding vortex interaction modes are investigated. 2. Mathematical formulas 2.1. The motions of the tandem foils
global phase shift horizontal and vertical forces and moment of each foil 𝑓1𝑖 = 2𝐹1𝑖 2𝐹 4𝐹 , 𝑓3𝑖 = 𝜌𝐶 𝑈3𝑖2 , 𝑓5𝑖 = 𝜌𝐶 2 5𝑈𝑖 2 𝜌𝐶 𝑈 2
∫ 𝐹3𝑖 (𝜏)ℎ̇ 𝑖 (𝜏)𝑑𝜏
power due to vertical motion
ℎ1 (𝑡) = ℎ0 sin(𝜔𝑡 + 𝜎), ℎ2 (𝑡) = ℎ0 sin 𝜔𝑡,
(1)
∫ 𝐹5𝑖 (𝜏)𝜃̇ 𝑖 (𝜏)𝑑𝜏
power due to rotational motion
𝜃1 (𝑡) = 𝜃0 sin(𝜔𝑡 + 𝜎 + 𝜀), 𝜃2 (𝑡) = 𝜃0 sin(𝜔𝑡 + 𝜀),
(2)
power due to vertical and rotational motions power coefficients due to vertical motion, rotational motion power consumption due to rotational motion power coefficient
where 𝜔 and 𝜀 are the circular frequency and phase difference respectively, the subscripts 1, 2 indicate the number of the foil. The hind foil will perform oscillatory motion in the vortex wake of the fore foil. There exists a global phase shift of motions between the fore and hind foils in terms of L and 𝜎 [20]. We have 𝐿 𝜓 = 2𝜋 +𝜎 (3) 𝑈𝑇
+𝜎
𝑡+𝑇
𝑡 𝑡+𝑇 𝑡
Pi = PLi + PMi 𝑐𝑃 𝐿𝑖 =
Fig. 1. Sketch of the tandem flapping foils in current.
Fig. 1 shows the tandem configuration of two foils in unbounded flow. The chord of the foil is denoted as C. The Cartesian coordinate system oxz locates at the mean position of the fore foil with x in the direction against the current and z pointing upwards. The hind foil locates after the fore foil. These two foils are denoted as foil 1 and foil 2, respectively. The motions of the foils in the uniform current are prescribed through sinusoidal heave and pitch motions. We have
𝜔ℎ0 𝜋𝑈 𝜓 = 2𝜋 𝑈𝐿𝑇 F1i ,F3i ,F5i
𝑆𝑡 =
𝑃𝐿𝑖 =
density of the fluid chord of the foil the maximum vertical extension of the foil velocity of the uniform current the distance from leading edge to rotational centre (xci , zci ) and lc /C = 1/3, i = 1, 2 indicates the fore and hind foils longitudinal distance of the two foils at their rotational centre circular frequency of oscillatory motion period of the oscillatory motion heave motion of the foil the heave amplitude pitch motion, the positive direction is defined as counterclockwise the pitch amplitude the phase difference of the heave and pitch motions the phase shift of the motions of the two foils effective attack angle at the rotational centre, dot over hi (t) means the derivative against time nominal amplitude of effective attack angle Strouhal number
𝑃𝐿𝑖 ; 𝑐𝑃 𝑀𝑖 0.5𝜌𝐶 𝑈 3
𝑐𝑃𝑖𝑛𝑀𝑖 = −𝑐𝑃 𝑀𝑖 𝑃𝑖 0.5𝜌𝐶 𝑈 3 𝑃𝑖 0.5𝜌𝑈 3 𝑑
𝑐𝑃 𝑖 =
𝜂𝑖 = 𝜂t = 𝜂1 + 𝜂2
=
𝑃𝑀𝑖 0.5𝜌𝐶 𝑈 3
power efficiency the total efficiency of the tandem foils cPL , cPM , cP and 𝜂 are the corresponding coefficient for single foil
Substituting 𝑆𝑡 = 𝜓 = 𝑆𝑡
𝜔ℎ0 𝜋𝑈
and 𝑇 =
2𝜋 𝜔
into Eq. (3), we have
𝐶 𝐿 +𝜎 ℎ0 2𝐶
(4)
The global phase shift is the function of Strouhal number St, heave amplitude h0 , longitudinal distance L and phase shift 𝜎. The vortex interactions of the foils depend on these parameters. We note that the vortex structure of the flapping foil is determined when the motion parameters St, h0 /C, 𝛼 0 and 𝜀 are prescribed. Different motion parameters will result in different vortex structure. For prescribed motions, the global phase shift, which unifies the L and 𝜎, will affect the efficiency of the hind foil significantly.
motion 𝜎 and the longitudinal distance L are less understood. The mechanism of the vortex interaction and the corresponding global phase shift require systematic investigations. The energy harvesting of the tandem foils in unbounded flow with typical oscillatory motions will be analysed in present analysis. For oscillatory foils, leading edge vortex shedding would be reduced or even eliminated even though larger attack angle is considered [21,22]. The leading edge vortex shedding would be delayed at higher Reynolds number and would result in the increase of energy extraction efficiency [10,23–25]. These justify the unsteady panel method for the simulation
2.2. Governing equations and the boundary value problem We consider the tandem foils in a uniform current with periodical heave and pitch motions. By assuming that the Reynolds number is high 203
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
enough and the effective attack angle ranges from small to medium, the boundary layer is thin and the vortex only sheds at the trailing edge. The inviscid flow around the foils is described by the velocity potential flow theory. The total velocity potential Φ in the presence of the current satisfies Laplace equation. We have ∇2 Φ = 0
The Kelvin theorem requires that the circulation of the fluid domain should be constant. The total circulation of the vortices and foils shall be zero. The summation of the shed vortex and the circulation over the foils shall be equal but opposite in sign. We note that the strength of each point vortex calculated through Eq. (16) will not change during the evolution of the vortex. Eqs. (14) and (16) guarantee the circulation of the fluid domain be zero and the Kelvin theorem be satisfied. The position of the point vortex can be updated through following equations
(5)
The boundary condition on the foil surface S0 can be written as ⇀ ⇀ 𝜕Φ = (−𝜃̇ 𝑖 𝑍, ℎ̇ 𝑖 + 𝜃̇ 𝑖 𝑋 𝑖 ) 𝑛 = −𝜃̇ 𝑖 𝑍𝑖 𝑛𝑥 + (ℎ̇ 𝑖 + 𝜃̇ 𝑖 𝑋𝑖 )𝑛𝑧 𝜕𝑛
(6)
⇀
where 𝑛 = (𝑛𝑥 , 𝑛𝑧 ) is the inward normal of the foil surface, 𝜃̇ 𝑖 and ℎ̇ 𝑖 are the rotational velocity and the vertical velocity respectively, and ⇀
𝑋 𝑖 = (𝑋𝑖 , 𝑍𝑖 ) = (𝑥 − 𝑥𝑐𝑖 , 𝑧 − 𝑧𝑐𝑖 ), (xci ,zci ) denotes the rotational centre of the corresponding fore or hind foil. In the far field ∇Φ = −Ui, 𝑥∕𝐶 >> 0 or |𝑧|∕𝐶 >> 0
where i indicates the horizontal unit vector. The velocity potential 𝜙 due to the disturbance the foil motion is introduced. We have (8)
Substituting Eq. (8) into Eqs. (5)–(7), we obtain the governing equations and boundary conditions ∇2 𝜙 = 0 in the f luid domain
(17)
𝑧𝜅(𝑗 ) (𝑡 + Δ𝑡) = 𝑧𝜅(𝑗 ) (𝑡) + 𝑤𝜅(𝑗 ) 𝑑𝑡
(18)
where u𝜅(j) and w𝜅(j) are the local horizontal and vertical velocities, which will be referred in the appendix. As the shed point vortex of the fore foil approaches the hind foil, it could penetrate the foil surface due to numerical error. A fence scheme is adopted to avoid the penetration [28]. The fence is an artificial surface that has an offset 𝛿 to the foil surface. When the distance of the point vortex to the foil surface is smaller than 𝛿, the point vortex will be placed on the nearest artificial fence. The offset 𝛿 over the foil surface is approximately the size of the element on the foil. To solve the boundary value problem, we have the integral equation based on Green’s third identity
(7)
Φ(𝑥, 𝑧) = 𝜙(𝑥, 𝑧) − 𝑈 𝑥 or 𝜙(𝑥, 𝑧) = Φ(𝑥, 𝑧) + 𝑈 𝑥
𝑥𝜅(𝑗 ) (𝑡 + Δ𝑡) = 𝑥𝜅(𝑗 ) (𝑡) + (−𝑈 + 𝑢𝜅(𝑗 ) )𝑑𝑡
(9)
Λ(𝑝)𝜙(𝑝) =
𝑆01 +𝑆02
𝜕𝜙 𝜕Φ = + 𝑈 𝑖 = (𝑈 − 𝜃̇ 𝑖 𝑍𝑖 )𝑛𝑥 + (ℎ̇ 𝑖 + 𝜃̇ 𝑖 𝑋𝑖 )𝑛𝑧 on foil surface 𝜕𝑛 𝜕𝑛
(10)
∇𝜙 = ∇Φ + 𝑈 𝑖 = 0 in the far f ield
(11)
+
∫
𝜕 𝐺(𝑝, 𝑞 ) 𝜕𝜙(𝑞) 𝜙(𝑞) − 𝐺(𝑝, 𝑞) ]𝑑𝑆 𝜕 𝑛𝑞 𝜕 𝑛𝑞
𝑆 𝑤 1 +𝑆 𝑤 2
The Kutta condition requires that the velocity is finite and continuous at the trailing edge. We have [8] 𝜕 𝜙+ || 𝜕 𝜙− || 𝜕𝜇 || − |⇀ = 𝜕𝑠 |⇀ 𝜕𝑠 ||⇀ 𝜕𝑠 | |𝑥𝑇 𝑥𝑇 𝑥𝑇
[
∫
+
𝑚 ∑ ( 𝑗=1
⇀ ⇀ ⇀ ⇀ 𝜕 𝐺(𝑝, 𝑞 ) 𝜇(𝑞) 𝑑𝑆 +𝜇(𝑥′𝑇 1 )𝐻(𝑝, 𝑥′𝑇 1 )+𝜇(𝑥′𝑇 2 )𝐻(𝑝, 𝑥′𝑇 2 ) 𝜕 𝑛𝑤
𝐻 (𝑝, 𝑞)𝜅𝑗1 + 𝐻 (𝑝, 𝑞)𝜅𝑗2
)
(19)
where Λ(p) is the solid angle at point p(x, z), G(p, q) = ln r, 𝑟 = √ (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2 , and q(𝜉, 𝜂) is the point on the integral boundary,
(12)
where 𝜙+ and 𝜙- are the potential on the upper and lower surface at the − trailing edge of the foil, 𝜇 = 𝜙+ 𝑤 − 𝜙𝑤 is the difference of the potential
𝐻(𝑝, 𝑞) = arctan
𝑥−𝜉 𝑧−𝜂
is obtained through the integral by parts.
⇀
on the upper and lower surface of the wake, 𝑥𝑇 = (𝑥𝑇 , 𝑧𝑇 ) denotes the 𝜕 trailing edge and 𝜕𝑠 denotes the partial derivative along the foil surface or vortex wake. The pressure p across the wake shall be continuous or p + = p − . In the form of Lagrangian, we have 𝑑𝜇 =0 𝑑𝑡 The above equation can be written in its discretised form as ⇀
⇀
𝜇(𝑥′𝑇 , 𝑡 + Δ𝑡) = 𝜇(𝑥𝑇 , 𝑡)
2.3. Calculation of the hydrodynamic forces Using Bernoulli equation, we have the pressure on each foil surface [ 𝑝𝑖 = −𝜌
(13)
(14)
where 𝑥′𝑇 = (𝑥′𝑇 , 𝑧′𝑇 ) is the end point of the dipole element connecting
𝐹𝑗𝑖 =
⇀
to the trailing edge and the strength of dipole 𝜇(𝑥𝑇 , 0) = 0 at first time step. In the time stepping scheme, the hind foil would penetrate the continuous vortex wake of the fore foil. To avoid the breaking of vortex dipole, the continuous vortex dipole is reduced into point vortex. Following the time stepping procedures in Xu and Wu [8], the dipole element can be converted into equivalent point vortex and a new dipole connecting to the trailing edge will develop at j + 1 time step. We define the vortex strength 𝛾 = 𝜕𝜇 . From Eq. (13) we have 𝜕𝑠
⇀
∫
𝑝𝑖 𝑛𝑗 𝑑𝑆 , 𝑗 = 1, 3, 5,
(21)
𝑠0
where (n1 ,n3 ,n5 ) = (nx ,nz ,−Xnz + Znx ). 3. Results and discussions 3.1. Typical vortex structure of single foil The vortex structure of a NACA 0012 foil performing sinusoidal heave and pitch motions is √ analysed. The uniform current velocity U = 4.43 and C = 2.0 or 𝑈 ∕ 𝑔𝐶 = 1.0. When solving the boundary integral equation, the linear distribution of the potential on each panel and the wake dipole is assumed. With finer element near the leading and trailing edges, numerical test shows that N = 60 elements over the foil surface can result convergent solutions. The time step interval Δt should make sure that the motion of point vortex at each time step is
𝑑𝛾 =0 (15) 𝑑𝑡 The strength of newly developed point vortex at time step j can be obtained in its discretised form ⇀
(20)
𝑑 where 𝑑𝑡 is the derivative fixed on the foil. Consequently, the force and moment can be calculated through the integration of pressure, we have
⇀
𝜅𝑗 (𝑡 + Δ𝑡) = 𝜇(𝑥𝑇 , 𝑡) − 𝜇(𝑥′𝑇 , 𝑡)
] 𝑑𝜙 1 − (𝑈 − 𝜃̇ 𝑖 𝑍𝑖 , ℎ̇ 𝑖 + 𝜃̇ 𝑖 𝑋𝑖 ) ⋅ ∇𝜙 + ∇𝜙∇𝜙 𝑑𝑡 2
(16)
204
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Fig. 2. (a) The typical vortex structure of a single foil and (b) the induced velocity by the vortices at t/T = 4.0 with St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, 𝜀 = 90°.
⇀
less than the average size of the mesh on the foil, or Δ𝑡 < 2𝐶∕(𝑁 𝑈 ),
Fig. 3. (a) The time histories of the vertical forces and (b) the comparison with the results of Kinsey and Dumas [20] which is referred as K&D in the figure, tandem NACA 0015 foils with L/C = 5.4, 𝜎 = 𝜋, St = 0.28, h0 /C = 1.0 and 𝜃 0 = 70° are simulated.
⇀
where 𝑈 = |𝑈 + 𝑢, 𝑤| is the velocity of the local flow. As shown in Fig. 2(a), the vortex sheds from the trailing edge successively. The structure of these vortices is a typical reversed Karman vortex street. The arrows in Fig. 2(b) indicate the vortex induced flow as the uniform current speed is excluded. When the hind foil moves in the vortex street, the average local fluid velocity will decrease and the vortex flow will change its effective attack angle significantly. The vortex structure in Fig. 2 suggests that the phase shift 𝜎 and the distance L would determine the relative position of the hind foil and the vortices of the fore foil. In the following section the tandem NACA 0012 foils with same motion parameters in Fig. 2 will be investigated.
3.3. Effects of phase shift 𝜎 The phase shift 𝜎 is the phase difference of the heave motions; it determines the relative vertical position of the hind foil and the vortices of the fore foil. We fix the motion of the two foils at St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, and 𝜀 = 90°. Simulations are carried out with L/UT = 0.5 and 𝜎 ranging from − 𝜋 to 𝜋. The energy efficiencies of the two foils at various 𝜎 are shown in Fig. 4(a). The horizontal axis on top of the figure is the corresponding global phase shift, which is calculated through Eq. (3). As shown in Fig. 4(a), there is no significant variation of the efficiency of the fore foil 𝜂 1 ; while the efficiency of the hind foil 𝜂 2 experiences one oscillation as 𝜎 increases from − 𝜋 to 𝜋. The efficiency of the hind foil is a periodical function of the phase shift 𝜎. The lowest efficiency of hind foil is found at 𝜎/2𝜋 = −0.45 or 𝜓/2𝜋 = 0.05; the highest efficiency is achieved at 𝜎/2𝜋 = −0.05 or 𝜓/2𝜋 = 0.45. To further investigate the power extraction of the hind foil, the power due to lifting force cPL2 and the consumption power due to the rotational motion 𝑐𝑃𝑖𝑛𝑀2 are presented in Fig. 4(b). Both cPL2 and 𝑐𝑃𝑖𝑛𝑀2 experience one oscillation. Power coefficient due to lifting force cPL2 is generally lower than that of a single foil, as shown in Fig. 4(b). At 𝜎/2𝜋 = −0.05, the curve of cPL2 finds its crest while the curve of 𝑐𝑃𝑖𝑛𝑀2 finds its trough. These result in the highest power coefficient and efficiency. Fig. 5 shows the typical time histories of the hind foil at 𝜎/2𝜋 = −0.5 (𝜓/2𝜋 = 0.0) and 𝜎/2𝜋 = −0.05 (𝜓/2𝜋 = 0.45). The amplitudes of horizontal force f1 and vertical force f3 at 𝜎/2𝜋 = −0.05 are larger than those of 𝜎/2𝜋 = −0.5, as shown in Fig. 5(a) and (b). The moment f5 at 𝜎/2𝜋 = −0.05 is smaller than that of 𝜎/2𝜋 = −0.5. The time histories of the forces of a single foil performing the same motion are presented in Fig. 5. Although the force amplitudes are larger than those of the tandem foils and therefore result in higher power coefficient cPL , the energy consumption due to rotational moment 𝑐𝑃𝑖𝑛𝑀 compensates the extracted power cPL . Smaller moment amplitude means the reduction of the consumed power. Consequently, the efficiency of hind foil at 𝜎/2𝜋 = −0.05 finds its peak.
3.2. Validation When tandem foils are considered, the hind foil will perform its motion in the vortex wake of the fore foil. To validate the numerical code, the interaction of tandem NACA 0015 foils with L/C = 5.4, 𝜀 = 𝜋, St = 0.28, h0 /C = 1.0 and 𝜃 0 = 70° is simulated firstly. Convergence study is carried out before extensive simulations. We used 60, 90 and 140 elements on the foil surface; the discrepancies of the forces are negligible. In the following simulations, 140 elements on each foil are adopted. The simulations stop after eight oscillations. The time history of vertical force of first five periods is shown in Fig. 3(a). The force of the hind foil is not affected in the first half oscillation. The effects of the vortex interaction on the hind foil become stable after the first oscillation. In the following figures, only one stable period will be shown. To compare the hydrodynamic forces of the fore and hind foils, the relative time t′ = t − 𝜎 for the fore foil is introduced. As shown in Fig. 3(b), the curves of present simulation are generally in good agreement with those of Kinsey and Dumas [20], whose results were obtained through solving two-dimensional Reynolds average Navier–Stokes equation at Reynolds number 5.0 × 105 . The snapshot of the flow field in Kinsey and Dumas [20] showed that there was no flow separation at the leading edge of the foil. The boundary layer was thin comparing with the foil. It can be treated as the high Reynolds number problem. The amplitude of present curves is slightly larger since the viscosity of the fluid has been ignored. 205
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Fig. 4. (a) The efficiency of the tandem foils and (b) the cPL2 and cPM2 of the hind foil at various 𝜎 at L/UT = 0.5, with St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, 𝜀 = 90° the dash line and dash dot line present cPL and cPM of single foil, respectively.
Fig. 5. The time history of (a)f12 , (b)f32 , (c)f52 of the hind foil at 𝜎/2𝜋 = −0.5 and 𝜎/2𝜋 = −0.05.
Typical vortex flow of the tandem foils is further investigated. Fig. 6 shows the flow velocity induced by the vortices at 𝜎/2𝜋 = −0.5. The vortices of the fore foil coincide with the vortices of the hind foil; stronger vortices are developed. The pattern of the vortex is similar with that of a single foil. When the hind foil performs the upstroke motion, it encounters the positive vortices of the fore foil; the hind foil encounters the negative vortices during its downstroke motion. The local flow velocity has been reduced and the effective attack angle decreases. Therefore the extracted power decreases significantly. The vortex pattern in Fig. 6 is named as synchronised mode. Fig. 7 gives the vortex flow at 𝜎/2𝜋 = −0.05. The hind foil encounters the negative vortices in the upstroke motion and it encounters positive vortices in the downstroke motion. The local flow velocity of the hind foil decreases due to the vortices. However, the effective attack angle increases. Therefore a rel-
atively higher efficiency has been achieved. The vortex of the tandem foils develops into two columns of vortices or vortex pairs, with positive vortices above the mean position and negative vortices below the mean position. The vortex pattern of the tandem foils is named as the interlaced mode. 3.4. The effects of longitudinal distance L In terms of the global phase shift, the longitudinal distance L has equivalent effects as 𝜎, as indicated in Eq. (4). The distance L determines the relative longitudinal position of the hind foil and the vortex of the fore foil. Its effects on the performance of the tandem foil are investigated. Again, the parameters of the motion are fixed at St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, and 𝜀 = 90°. Simulations are carried out 206
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Fig. 6. The flow induced by the vortices at 𝜓/2𝜋 = 0, (a) t/T = 0, (b) t/T = 1/8, (c) t/T = 1/4, (d) t/T = 3/8.
Fig. 7. Flow induced by the vortices at 𝜓/2𝜋 = 0.45, (a) t/T = 0, (b) t/T = 1/8, (c) t/T = 1/4, (d) t/T = 3/8.
with 𝜎 = −𝜋 and L/UT ranging from 0.5 to 1.0. The efficiency of the two foils at various L is shown in Fig. 8. The efficiency of the fore foil is not affected by the distance L obviously; while the efficiency of the hind foil fluctuates with L. The peak of the efficiency of the hind foil is as high as that of the fore foil at L/UT = 0.9 (𝜓/2𝜋 = 0.4); the lowest efficiency is found at L/UT = 0.55 (𝜓/2𝜋 = 0.05). Comparing Fig. 8(a) and Fig. 4(a), the curves in these two figures are resemble. The highest and lowest efficiencies are found with similar global phase shift 𝜓. This suggests that the effects of distance L have similar effects as the phase shift 𝜎. The curves of cPL2 and 𝑐 𝑖𝑛 are further shown in Fig. 8(b). These 𝑃 𝑀2 two waving curves of cPL2 and 𝑐 𝑖𝑛 find the peak and trough at L/UT 𝑃 𝑀2 ≈ 0.9 respectively, which result in the highest efficiency of the hind foil. The typical time histories of the forces of hind foil at L/UT = 0.5 and L/UT = 0.95 are further presented in Fig. 9. The horizontal and vertical forces of L/UT = 0.95 in Fig. 9(a) and (b) have larger amplitudes than those of L/UT = 0.5; the amplitude of moment at L/UT = 0.95 is smaller
than that of L/UT = 0.5. As a result, the efficiency with L/UT = 0.95 is higher than that of L/UT = 0.5. Comparing with the curves of a single foil, the amplitudes of the forces of tandem foils are smaller. However, higher efficiency is achieved due to the decrease of the energy consumption of the pitch moment at L/UT = 0.95. It is expected that the vortex interactions would be the same as those in Figs. 6 and 7. The synchronised mode and interlaced mode will form at L/UT = 0.5 (𝜓/2𝜋 = 0) and L/UT = 0.95 (𝜓/2𝜋 = 0.45) respectively. 4. Concluding remarks The energy harvesting of flapping foils with tandem configurations in uniform current has been simulated through the boundary elements. The hind foil is moving in the vortex wake of the fore foil. The effects of vortex–foil interaction on the energy extraction performance are analysed. There are strong interactions of foils and vortices. The efficiency of the fore foil is not affected significantly; however the power 207
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Fig. 8. (a) The efficiency of the tandem foils and (b) the cPL2 and cPM2 of the hind foil against L at 𝜎 = −𝜋, with St = 0.3, h0 /C = 1.0, 𝛼 0 = 15°, 𝜀 = 90°.
Fig. 9. The time history of (a)f12 , (b)f32 , (c)f52 of the hind foil at L/UT = 0.5 and L/UT = 0.95 in Fig. 8.
efficiency of the hind foil fluctuates with the global phase shift 𝜓. The lowest efficiency of the hind foil is around 𝜂 2 = 0.05 and its highest efficiency is 𝜂 2 = 0.24. The efficiency curves against L and 𝜎 resemble each other. It is interesting to find that the peaks of efficiency are found at global phase 𝜓/2𝜋 ≈ 0.45, where the “interlaced mode” can be observed. The lowest efficiency appears when 𝜓/2𝜋 ≈ 0 or 𝜓/2𝜋 ≈ 1.0, where the “synchronised mode” is found. The global phase shift 𝜓 indicates the vortex interaction modes and 𝜎 and L have equivalent effects on the efficiency. Present simulations are limited to one set of typical motion parameters. The effects of the motion parameters, including St, h0 , 𝛼 0 and 𝜀, have been ignored. We noted that the change of these parameters will cause the variation of the vortex structure. Even though the same global phase shift is adopted, it will not necessarily cause the increase or decrease of the efficiency due to the complexity of the foil–vortex interactions. For instance, when small Strouhal number is adopted, the longitudinal distance of the vortices will increase and the strength of
the vortices will decrease; the variation of the effective attack angle will become less significant. When large Strouhal number is concerned, the distance of the vortices will decrease. The average effective attack angle will not increase or decrease significantly if the vertical and longitudinal distances of the vortices are smaller than the size of the foil. The effects of different motion parameters, including St, h0 , 𝛼 0 and 𝜀, require further investigations. Present study would be highly important for the design and construction of hydro-kinematic turbines based on flapping tandem foils. A proper distance and phase shift of the motions can be selected to assure the efficiency of the tandem foil system. For fixed longitudinal distance, the foils can adapt the phase shift of the motions to obtain the most favourable global phase shift. When the phase shift of the motions is fixed, we could adjust the longitudinal distance. The tandem configuration of the foils can be implemented through two crankshaft structures like those in Xu et al. [29], with the shaft of two sets of foil-crankshaft mechanisms being linked through chain wheels.
208
G.D. Xu, W.H. Xu
Engineering Analysis with Boundary Elements 82 (2017) 202–209
Acknowledgement
Reference
The research is under the support of National Science Foundation of China (Grant nos. 11602067 and 51479044) and State Key laboratory of Ocean Engineering (Shanghai Jiao Tong University, Grant No. 1413).
[1] Schiermeier Q, Tollefson J, Scully T, Witze A, Morton O. Electricity without carbon. Nature 2008;454(7206):816–23. [2] Young J, Lai J, Platzer M. 2014 A review of progress and challenges in flapping foil power generation Prog Aerosp Sci 67, 2–28. [3] Xiao Q, Zhu Q. A review on flow energy harvesters based on flapping foils. J Fluids Struct 2014;46:174–91. [4] McKinney W, DeLaurier J. Wingmill: an oscillating-wing windmill. J Energy 1981;5(2):109–15. [5] Kinsey T, Dumas G, Lalande G, Ruel J, Mehut A, Viarouge P, et al. Prototype testing of a hydrokinetic turbine based on oscillating hydrofoils. Renew Energy 2011;36:1710–18. [6] Platzer M, Ashraf M, Young J, Lai J. Development of a new oscillating-wing wind and hydropower generator. In: Proceedings of the forty-seventh AIAA aerospace sciences meeting. Reston, Virginia: American Institute of Aeronautics and Astronautics; 2009. [7] Anderson JM, Streitlien K, Barrett DS, Triantafyllou MS. Oscillating foils of high propulsive efficiency. J Fluid Mech 1998;360:41–72. [8] Xu GD, Wu GX. Boundary element simulation of inviscid flow around an oscillatory foil with vortex sheet. Eng Anal Bound Elem 2013;37:825–35. [9] Jones KD, Platzer MF. Numerical computation of flapping-wing propulsion and power extraction. AIAA J 1997:97–0826 35th AIAA Aerospace Sciences Meeting, Reno, Nevada, Jan. [10] Jones KD, Lindsey K, Platzer MF. An investigation of the fluid–structure interaction in an oscillating-wing micro-hydropower generator. Fluid–structure interaction II. Southampton, UK: Wessex Institute of Technology Press; 2003. p. 73–82. [11] Dumas G, Kinsey T. Eulerian simulations of oscillating airfoils in power extraction regime. In: Rahman B, editor. Proceedings in advances in fluid mechanics VI. WIT Press; 2006. p. 245–54. [12] Zhu Q. Optimal frequency for flow energy harvesting of a flapping foil. J Fluid Mech 2011;675:495–517. [13] Simpson BJ, Licht S, Hover FS, Triantafyllou MS. Energy extraction through flapping foils. In: Proceedings of the twenty-seventh international conference on offshore mechanics and arctic engineering. OMAE; 2008. p. 2008–58043. [14] Ashraf MA, Young J, Lai JCS, Platzer MF. Numerical analysis of an oscillating-wing wind and hydropower generator. AIAA J 2011;49(7):1374–86. [15] Lu K, Xie Y, Zhang D. Nonsinusoidal motion effects on energy extraction performance of a flapping foil. Renew Energy 2014;64:283–93. [16] Platzer M, Ashraf M, Young J, Lai J. Extracting power in jet streams: pushing the performance of flapping wing technology. In: Proceedings of the twenty-seventh congress of the international council of the aeronautical sciences. Nice France: International Council of the Aeronautical Sciences, September 19–24; 2010 Paper 2010-2.9.1. [17] Xiao Q, Liao W, Yang SC, Peng Y. How motion trajectory affects energy extraction performance of a biomimic energy generator with an oscillating foil? Renew Energy 2012;37:61–75. [18] Liu WD, Xiao Q, Cheng F. A bio-inspired study on tidal energy extraction with flexible flapping wings. Bioinspiration Biomim 2013;8:1–16. [19] Lehmann F-O. Wing–wake interaction reduces power consumption in insect tandem wings. Exp Fluids 2009;46(5):765–75. [20] Kinsey T, Dumas G. Optimal tandem configuration for oscillating-foils hydrokinetic turbine. J Fluids Eng 2012;134(3):031103-1–031103-11. [21] Dickinson MH, Lehmann FO, Sane SP. Wing rotation and the aerodynamic basis insect flight. Science 1999;284:1954–60. [22] Ellington CP, Vanderberg C, Wilmott A, Thomas A. Leading edge vortices in insect flight. Nature 1996;384:626–30. [23] Campobasso MS, Piskopakis A, Drofelnik J, Jackson A. Turbulent Navier–Stokes analysis of an oscillating wing in a power-extraction regime using the shear stress transport turbulence model. Comput Fluids 2013;88:136–55. [24] Kinsey T, Dumas G. Parametric study of an oscillating airfoil in a power-extraction regime. AIAA J 2008;46:1318–30. [25] Kinsey T, Dumas G. Optimal operating parameters for an oscillating foil turbine at Reynolds number 500,000. AIAA J 2014;52(9):1885–95. [26] La Mantia M, Dabnichki P. Unsteady panel method for flapping foil. Eng Anal Bound Elem 2009;33:572–80. [27] Triantafyllou MS, Hover FS, Techet AH, Yue DKP. Review of hydrodynamic scaling laws in aquatic locomotion and fishlike swimming. J Am Soc Mech Eng 2005;58:226–37. [28] Pan YL, Dong XX, Zhu Q, Yue DKP. Boundary-element method for the prediction of performance of flapping foils with leading-edge separation. J Fluid Mech 2012;698:446–67. [29] Xu GD, Xu WH, Dai J. Numerical and experimental study of a flapping foil generator. Appl Ocean Res 2017;63:242–50. [30] Li S, Han R, Zhang AM, Wang QX. Analysis of pressure field generated by a collapsing bubble. Ocean Eng 2016;117:22–38.
Appendix Once the boundary value problem has been solved, the potential and normal velocity on the integral boundary are known. However the velocity of a point in the fluid domain cannot be obtained directly [30]. Here we shall derive a numerical scheme to calculate the local velocity. Combine the similar terms of the two foils in Eq. (19), we have [ Λ(𝑝)𝜙(𝑝) =
𝜕 ln 𝑟𝑝𝑞
∫
𝜕 𝑛𝑞
𝑆0 ⇀
𝜙(𝑞) − ln 𝑟𝑝𝑞
+𝜇(𝑥′𝑇 ) arctan
] 𝜕 ln 𝑟𝑝𝑞 𝜕𝜙(𝑞) 𝑑𝑆 + 𝜇(𝑞)𝑑𝑆 ∫ 𝜕 𝑛𝑤 𝜕 𝑛𝑞 𝑆𝑤
𝑥 − 𝑥′𝑇 𝑧 − 𝑧′𝑇
( ) 𝑥 − 𝑥𝑗 + arctan 𝜅𝑗 𝑧 − 𝑧𝑗 𝑗=1 𝑚 ∑
(22)
where S0 , Sw and (xj , zj ) denote the foil surface, wake panel and position of the point vortex, respectively. The solid angle Λ(p) = 2𝜋 if p is within the fluid domain. Applying the integral by part to the first term in the brackets and the second integral, we have 2𝜋𝜙(𝑝) = −
𝜕 𝜙(𝑞 ) 𝑥−𝜉 𝜕 𝜙(𝑞 ) arctan ln 𝑟𝑝𝑞 𝑑𝑆 𝑑𝑆 − ∫ 𝜕𝑠 𝑧−𝜂 𝜕𝑛
∫
𝑆0
+
𝑆0
) 𝑚 ( ∑ 𝑥 − 𝑥𝑗 𝜕 𝜇(𝑞 ) 𝑥−𝜉 𝑑𝑆 + arctan arctan 𝜅𝑗 ∫ 𝜕𝑠 𝑧−𝜂 𝑧 − 𝑧𝑗 𝑗=1
(23)
𝑆𝑤
Applying the derivative against x and z, we have 2𝜋
𝜕𝜙(𝑝) 𝜕 𝜙(𝑞 ) 𝑧−𝜂 =− 𝑑𝑆 ∫ 𝜕𝑥 𝜕𝑠 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2 𝑆0
−
∫
𝜕 𝜙(𝑞 ) 𝑥−𝜉 𝑑𝑆 𝜕𝑛 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2
∫
𝜕 𝜇(𝑞 ) 𝑧−𝜂 𝑑𝑆 𝜕𝑠 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2
𝑆0
+
𝑆𝑤
+
𝑚 ∑
(
(𝑥 − 𝑥𝑗 )2 + (𝑧 − 𝑧𝑗 )2
𝑗=1
2𝜋
𝑧 − 𝑧𝑗
) 𝜅𝑗
(24)
𝜕𝜙(𝑝) 𝜕 𝜙(𝑞 ) 𝑥−𝜉 = 𝑑𝑆 ∫ 𝜕𝑧 𝜕𝑠 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2 𝑆0
−
∫
𝜕 𝜙(𝑞 ) 𝑧−𝜂 𝑑𝑆 𝜕𝑛 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2
∫
𝜕 𝜇(𝑞 ) 𝑥−𝜉 𝑑𝑆 𝜕𝑠 (𝑥 − 𝜉)2 + (𝑧 − 𝜂)2
𝑆0
+
𝑆𝑤
−
𝑚 ∑ 𝑗=1
(
𝑥 − 𝑥𝑗 (𝑥 − 𝑥𝑗 )2 + (𝑧 − 𝑧𝑗 )2
) 𝜅𝑗
(25)
The velocity at point p in the fluid domain can be obtained through 𝑝) 𝜕𝜙(𝑝) (𝑢, 𝑤) = ( 𝜕𝜙( , 𝜕𝑧 ). 𝜕𝑥 209