Numerical analysis of a tidal current generator with dual flapping wings

Accepted Manuscript Numerical analysis of a tidal current generator with dual flapping wings

Penglei Ma, Yong Wang, Yudong Xie, Zhipu Huo PII:

S0360-5442(18)30855-7

DOI:

10.1016/j.energy.2018.05.035

Reference:

EGY 12868

To appear in:

Energy

Received Date:

07 January 2018

Revised Date:

21 March 2018

Accepted Date:

05 May 2018

Please cite this article as: Penglei Ma, Yong Wang, Yudong Xie, Zhipu Huo, Numerical analysis of a tidal current generator with dual flapping wings, Energy (2018), doi: 10.1016/j.energy. 2018.05.035

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Numerical analysis of a tidal current generator with dual flapping wings

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

Penglei Ma1,2, Yong Wang1,2,*, Yudong Xie1,2,*, Zhipu Huo1,2 1. School of Mechanical Engineering, Shandong University, Jinan 250061, China 2. Key Laboratory of High-efficiency and Clean Mechanical Manufacture (Shandong University), Ministry of Education, Jinan 250061, China *Corresponding author China. [email protected] (Yong Wang); [email protected] (Yudong Xie) Abstract: Flapping wings, inspired by the mechanism of birds and fish, can act as generators to harvest energy from tidal currents. The hydraulic system is simplified as a spring-damper system to establish the coupling equations relating to the wing motion and the hydrodynamic forces. To provide guidance for design of a fully flow-induced flapping wings energy harvesting system, the behaviors of both system response and energy extraction performance are analyzed using two-dimensional numerical approach. Depending on the rotary actuator radius R, and the volume ratio β between the cylinder and rotary actuator, three distinguishable behaviors are observed in the system response and energy extraction performance. At larger R and smaller β, the dual wings tend to undergo a damped reduction flapping motion because the pitching motion consumes a significant amount of energy. Both decreasing R and increasing β can reduce the energy consumption of the pitching motion, and thus allow the dual wings to achieve a sustainable flapping motion. Although an irregular response can achieve a self-sustained flapping motion, it is unfavorable owing to its unstable power output. The regular response essential for stable energy harvesting is realized over a range of coupling parameters. The energy extraction performance of the system is closely associated with β but also slightly depends on R. Keywords: tidal current generator; dual flapping wings; energy harvesting; hydraulic coupling system; system response; flow-induced;

18

1. Introduction

19 20 21 22 23 24 25 26 27 28 29 30 31

Tidal energy has been considered as an ideal alternative energy source to reduce dependence on fossil fuels. Conventionally, turbines with rotating blades are used for harnessing tidal energy. However, these devices are usually gigantic and have many negative effects on the environment and economy. Inspired by birds and fish, applying flapping wings to generators has gradually attracted increasing attention [1, 2]. In their early stages, studies on flapping wings mostly focused on the propulsion mechanism. If a wing is free to flap in both the pitching and heaving degrees of freedom, power can be transferred from the air or water to the wing. The application of flapping wings for extracting energy from a uniform flow was first proposed by Mackinney and Delaurier [3]. With the increasing importance of renewable energy, this concept has attracted more interest in recent years. Following the study by Mackinney and Delaurier, Jones and Platzer studied the aerodynamics of a flapping wing using an unsteady panel code [4, 5]. Their results revealed that the wing could change from energy consumption to energy generation at a fixed flapping frequency and heaving amplitude 1

ACCEPTED MANUSCRIPT 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70

if the pitching amplitude was sufficiently large. Recently, Kinsey et al. have performed systematic investigations using both numerical approaches and prototype tests [6, 7]. The effects of viscosity (Reynolds number), geometric parameters (wing shape and location of pitching axis), and motion parameters (pitching amplitude, heaving amplitude and flapping frequency) on the energy extraction performance were investigated. Extensive research on flapping wings can be found in the studies conducted by Xiao et al. [8], Xie et al. [9, 10], and Tian et al. [11, 12]. Dual wings in a suitable configuration have better energy extraction performance than a single wing. Kinsey and Dumas studied the interaction between the wakes of the upstream wing and the downstream wing and determined that the relative position of the downstream wing is a critical parameters in energy extraction performance [13]. Ko et al. carried out numerical analyses and experimental investigations of a tidal current generator with a mirror-type configuration of front and rear swing flapping wings [14]. Xu et al. demonstrated that tandem flapping wings have a better energy extraction behavior for an optimal frequency and suitable separation distance [15]. Karbasian et al. studied the potential for power harvesting with flapping wings in tandem. Their results demonstrated that it is possible to improve the performance of a flapping wing system by suing multistage wings arranged in tandem [16, 17]. In the above studies, the wings move with a prescribed motion. Nevertheless, it is hard to simultaneously control the pitching and heaving motions of a wing while harvesting energy from a tidal current. Recently, some researchers have focused on a semi-activated system with an actuated pitching motion but a free heaving motion [18, 19]. Using a linearized thin plate model, Zhu et al. demonstrated that this mechanism has a positive power extraction over a large range of mechanical and operational parameters [20]. In addition, computations were carried out using a Navier-Stokes model, the results indicated that the performance of the semi-activated system is closely related to vorticity control mechanisms, particularly in the interaction between the wing and the leading edge vortex [21]. Using the same method, Peng et al. investigated a device consisting of a single flapping wing mounted on a damper and a rotational spring, which was shown to achieve self-induced oscillating motions. Many studies of flapping wing power generation concentrate on fully prescribed motion or prescribed pitching motion with flow-driven heaving motion. Although the harvesting energy of a flapping wing is greater than the power consumption of the control system, it is contrary to the objective that energy harvesting from tidal current makes use of an external control to achieve the wings flapping motion. To achieve wings with fully flow-induced flapping motion, a hydraulic system is employed to couple the pitching motion and heaving motion of dual wings. The primary objective is to observe the system response to different coupling parameters. The rest of the study is organized in the following structure. First, the physical problem is described and hydraulic coupling system is presented. The numerical method is also briefly described. Next, the responses of the coupling system and its energy harvesting performance are discussed. Finally, several conclusions are given.

71

2. Problem description and methodology

72

2.1 Problem description

73

In most existing studies, the wing follows a prescribed motion, which may be either a fully 2

ACCEPTED MANUSCRIPT 74 75 76 77 78 79 80 81 82 83 84

prescribed motion or a semi-activated motion. For a single wing, it is difficult to achieve a real sense of fully flow-induced and self-sustained flapping motion during energy harvesting. An external control system is usually used to achieve the flapping motion of the wings. This control system makes the apparatus more complicated. In addition, although the harvesting energy of a flapping wing is greater than the power consumption of the control system, it is counter to the objective of harvesting energy from tidal currents. In this study, a hydraulic system is used to couple the pitching motion and heaving motion of wings instead of an external control system [22]. Rather than actively governing the wing motion, stability of the flapping motion for dual wings will be completely induced by the flow. Dual wings with a parallel configuration are considered in this study because there is no vortex interaction between the dual wings in this configuration.

cylinder check valve

motor

wing 2

generator G

wing 1

rotary actuator

piston rod

power output system

85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100

Fig. 1 Schematic of hydraulic system to couple the dual wings Fig. 1 shows the hydraulic system used to couple the dual wings. This figure illustrates the hydraulic system setup and does not reflect the spatial configuration of the dual wings. For the purpose of illustration, components, such as control valves and accumulators, are not displayed in Fig. 1. In Fig. 1, Wing 1 is located at the lowest point of its heaving motion. The lift acting on the wing cannot sustain the heaving motion because of it has a zero pitching angle. Wing 2 is at the dead point of its pitching motion. The torque loaded on the wing cannot pitch it in an ideal direction. However, the lift acting on Wing 2 can propel it in a downward direction. Hydraulic power is transmitted from the heaving motion of Wing 2 to the pitching motion of Wing 1. The lift acting on Wing 1 will increase with increasing pitching angle and then cross the dead point of its heaving motion. Similarly, the heaving motion of Wing 1 will help Wing 2 overcome the dead point of its pitching motion. Part of heaving motion of one wing is used to propel the pitching motion of the other and the rest of the heaving motion is used to generate power. The non-rod piston chamber of the cylinder comprises a hydraulic circuit with a power generation system, as illustrated in Fig. 1. When moving the piston, the fluid in the chamber will rotate the motor, and thus produce power. The check valve provides a consistent direction of fluid flow to the motor. 3

ACCEPTED MANUSCRIPT 101

2.2 Motion kinematics

102 103 104

Before building the coupling equations, the parameters related to a wing are identified, as illustrated in Fig. 2. A Cartesian coordinate system is built in which x is the direction of flow and y is the direction of the wing heaving motion.

y

α

x

M

U∞

L

D

y

F

θ

105 106 107 108 109 110 111

Fig. 2 Sketch of wing hydrodynamics and pivot location In Fig. 2, L and D are the chord length and wing thickness, respectively. For most flapping wings in hydrodynamic applications, the pitching center is located in the one third length of the chord from the leading edge [2]. F and M are the lift force and pitching torque loaded on the wing, respectively. The drag force is not considered because it is irrelevant to the system response. The corresponding force and torque coefficients (CF and CM, respectively) are defined as follows: 𝐶𝐹 = 𝐶𝑀 =

112 113 114 115 116 117 118 119 120 121 122

𝑀

𝐹

(1)

1 2 2 2𝜌𝑈∞𝐿 𝑆,

(2)

where ρ is the fluid density, and S is the wing-span. Here, a wing-span equal to one unit length is used for the 2D model. The effective angle of attack α is dependent on the heaving-induced pitching angle (Fig. 2), as follows α = arctan ( ‒ 𝑦 𝑈∞) ‒ 𝜃,

(3)

where 𝑦 is the heaving velocity, and 𝜃 is the pitching angle. To build the equation relating to the hydrodynamics and motion, the hydraulic system is simplified as a spring-damper system. Following the study by Zhu et al., the cylinder is regarded as a damper [21, 23]. The spring of the cylinder and rotary actuator are negligible in this study. The damping of the damper is denoted as c. Before establishing the mathematical formulation, some parameters are further defined to unify the dimension. ∗

𝑐 = 123 124

1 2 2𝜌𝑈∞𝐿𝑆

𝑐

1 2 2𝜌𝑈∞𝐿𝑆

.

(4)

The gravity of the wing is considered, whereas the mass of the wing is insignificant compared to the added mass. The total mass, m, and total torque of inertial, I, are described as follows: ∗

𝑚 = 4

𝑚

1 2 2𝜌𝑈∞𝐿𝑆

(5)

ACCEPTED MANUSCRIPT 𝐼

∗

𝐼 = 125 126 127 128

1 2 2 2𝜌𝑈∞𝐿 𝑆.

(6)

As shown in Fig. 1, each wing can flap in two degrees of freedom: the pitching motion, 𝜃(t), and heaving motion, 𝑦(t). The system as whole also has two degrees of freedom because one wing’s heaving motion is coupled with the other wing’s pitching motion. Consequently, the coupling equations relating to hydrodynamics and wing motion can be given as follows [24]:

{ {

∗

∗

𝐶𝐹1 = 𝑚 𝑦1 + 𝑐 𝑦1 + 𝐶𝐹𝑀12

(7)

∗

𝐶𝑀1 = 𝐼 𝜃1 ‒ 𝐶𝐹𝑀21 ∗ 𝑅 ∗

∗

𝐶𝐹2 = 𝑚 𝑦2 + 𝑐 𝑦2 + 𝐶𝐹𝑀21 ∗

𝐶𝑀2 = 𝐼 𝜃2 + 𝐶𝐹𝑀12 ∗ 𝑅

,

(8)

129

where the subscripts 1 and 2 indicate the corresponding parameters for Wing 1 and Wing 2,

130

respectively; 𝜃 and 𝑦 are the angular acceleration and heaving acceleration, respectively; 𝐶𝐹𝑀12 is the

131

component force coefficient of Wing 1 imposing on Wing 2; 𝐶𝐹𝑀21 is component force coefficient

132 133 134 135 136

of wing 2 imposing on Wing 1; and R is the parameter associated with the radius of the rotary actuator, which has dimension of m. In the condition where there is no leakage of the hydraulic system, there is a linear relationship between one wing’s heaving motion and the other wing’s pitching motion, and which can be expressed as:

{

𝑦1 =‒ 𝛽𝜃2 𝑦2 = 𝛽𝜃1 ,

(9)

137 138

where 𝛽 is a factor associated with the volume ratio between the cylinder and the rotary actuator, expressed in m/rad.

139

The individual power coefficient, 𝐶𝑃, being harvested from the current is the sum of the heaving

140

contribution, 𝐶𝑃𝐹, and pitching contribution, 𝐶𝑃𝑀: 𝐶𝑃𝐹𝑖 = 𝐶𝑃𝑀𝑖 = 𝐶𝑃𝑖 =

2 3

𝜌𝑈∞𝐿𝑆 1

𝐶𝑃𝑖 = 𝐶𝑃𝐹𝑖 + 𝐶𝑃𝑀𝑖 = 𝑇 ∙ 141 142

2

(𝑡) ∙ 𝑦𝑖

𝐹 3 𝜌𝑈∞𝐿𝑆 𝑖

2

(10)

𝑀𝑖(𝑡) ∙ 𝜃𝑖

(11)

[𝐹𝑖(𝑡) ∙ 𝑦𝑖 + 𝑀𝑖 ∙ 𝜃𝑖].

(12)

3

𝜌𝑈∞𝐿𝑆

2 3 𝜌𝑈∞𝐿𝑆

[∫𝑇0𝐹𝑖(𝑡) ∙ 𝑦𝑖d𝑡 + ∫𝑇0𝑀𝑖 ∙ 𝜃𝑖d𝑡],

(13)

where T is the flapping period, and 𝐶𝑃 is the individual cycle-averaged power coefficient. The instantaneous overall power coefficient 𝐶𝑃𝑂 and cycle-averaged overall power coefficient,

5

ACCEPTED MANUSCRIPT 143

144 145

𝐶𝑃𝑂, are as follows: 𝐶𝑃𝑂 = 𝐶𝑃1 + 𝐶𝑃2

(14)

𝐶𝑃𝑂 = 𝐶𝑃1 + 𝐶𝑃2 .

(15)

The index of the energy harvesting performance is given by energy harvesting efficiency, η, which is expressed as: 1

1

𝜂 = 2𝜌𝑈∞𝐿𝑆 ∙ 𝐶𝑃𝑂 2𝜌𝑈∞𝑆𝐻 = 𝐶𝑃𝑂 ∙ 𝐿 𝐻 3

3

,

(16)

146 147 148

where H is the overall extent of the wing’s heaving motion, calculated as the distance between its highest position and lowest position. It depends primarily on the heaving amplitude but is also slightly influenced by the pitching amplitude and wing thickness.

149

3. Numerical approach

150

3.1 Governing equation and code setting

151

To compute the hydrodynamic forces, 𝐹𝑋, 𝐹𝑌, and torque, M, the solution of the Navier-Stokes

152 153

equations is more applicable at a large pitching angle. The governing equations for the unsteady incompressible flow around the wing are the 2D Navier-Stokes equations as follows: ∂𝑢𝑥 ∂𝑥 ∂𝑢𝑥 ∂𝑡 ∂𝑢𝑦 ∂𝑡

+

∂𝑢𝑦 ∂𝑦

(17)

=0

∂𝑢𝑥

∂𝑢𝑥

1∂𝑝

2

(18)

∂𝑢𝑦

∂𝑢𝑦

1∂𝑝

2

(19)

+ 𝑢𝑥 ∂𝑥 + 𝑢𝑦 ∂𝑦 = 𝐹𝑥 ‒ 𝜌∂𝑥 + 𝑣∇ 𝑢𝑥 + 𝑢𝑥 ∂𝑥 + 𝑢𝑦 ∂𝑦 = 𝐹𝑦 ‒ 𝜌∂𝑦 + 𝑣∇ 𝑢𝑦,

154

where 𝑢𝑥 and 𝑢𝑦 are the velocity in the x and y directions, respectively; 𝐹𝑥 and 𝐹𝑦 are the body force

155 156 157 158 159 160 161 162

in the x and y directions, respectively; p is the pressure, and v is the viscosity. The governing equations are solved using Fluent software based on the finite volume method [13, 15]. A grid model applied for computation is built with dynamic mesh and a sliding interface. The mesh model has been described in detail in previous studies [25-27]. Kinsey and Dumas found that the choice of turbulence affects the cycle-averaged performance indicators slightly [28]. Additionally, the Spalart-Allmaras turbulence model has a lower requirement for the mesh quality and lower time cost. Therefore, the selection of the one-equation Spalart-Allmaras turbulence model seems appropriate for computation in this study.

163

3.2 Iteration scheme

164 165 166

To transmit the force and torque, a user defined function is compiled to control the motion of the wings. An iteration scheme is performed to solve the coupling equations. The major steps are summarized as follows.

167

Step 1: In the initial stage, the initial pitching amplitude, 𝜃0, of Wing 1 and the initial heaving angle, 6

ACCEPTED MANUSCRIPT 168

𝑦0, of Wing 2 are provided. All other parameters are set to zero.

169 170 171 172

Step 2: The Navier-Stokes equation is solved using Fluent software. The force and torque exerted on the wing by the flow are calculated. Step 3: According to Eqs. (7), (8) and (9), the heaving acceleration and pitching acceleration can be expressed as: ∗ 𝑛

𝑛

173

∗ 𝑛

𝑛

(20) 𝑛

‒ 𝑅𝐶𝐹1 + 𝑅𝑐 𝑦1 + 𝐶𝑀2 𝑛 𝜃2 = ∗ ∗ 𝑅𝛽𝑚 + 𝐼

174

𝑛

𝑛

175

𝑦1 =‒ 𝛽𝜃2

176

𝑦2 = 𝛽𝜃1,

177 178 179 180

𝑛

𝑅𝐶𝐹2 ‒ 𝑅𝑐 𝑦2 + 𝐶𝑀1 𝑛 𝜃1 = ∗ ∗ 𝑅𝛽𝑚 + 𝐼

𝑛

𝑛

(21) (22) (23)

where the superscript n indicates corresponding parameters at time step n. To get the values of 𝜃 and 𝑦, the parameters on the right-hand side of the Eqs. (20) and (21) are substituted into the equations. Step 4: To update the position of wings, the classical Runge-Kutta method is used to calculate the displacement and moving velocity of the wings at time step n+1, as follows 𝜃

𝑛+1

𝜃 181

where

182

{

𝑛

𝑛

1

= 𝜃 + 6(𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4)∆𝑡 𝑛+1

𝑛

1

= 𝜃 + 6(𝑙1 + 2𝑙2 + 2𝑙3 + 𝑙4)∆𝑡,

(24) (25)

𝑛

𝑘1 = 𝜃 , 𝑙1 = 𝜃 𝑙 ∆𝑡 𝑛 𝑛 2 ,𝑙2 = 𝜃 + 𝑙1∆𝑡 2 𝑘2 = 𝜃 + 1 . 𝑙2∆𝑡 𝑛 𝑛 𝑙 ∆𝑡 2 ,𝑙3 = 𝜃 + 2 2 𝑘3 = 𝜃 + 𝑛

𝑛

𝑘4 = 𝜃 + 𝑙3∆𝑡 , 𝑙4 = 𝜃 + 𝑙3∆𝑡

183 184 185

The heaving velocity is calculated according to Eq. (9). The time step, ∆t, is sufficiently small to provide a high accuracy. Step 5: The above steps are repeated to complete the iteration.

186

3.3 Validation study

187 188 189 190 191

The grid number and time-step independence have a significant effect on the computational results. In a previous study, computations were carried out to verify the meshing reliability, and the results demonstrated that the grid quantity and time-step affect the computation results slightly [27]. Thus, the same mesh model is used for the computation in this study. This model has 105 cells (700 nodes on the wing) and a time step of 0.005s.

7

ACCEPTED MANUSCRIPT 4

Deng et al present study

3 2

CF

1 0 -1 -2 -3 -4 0.0

0.2

0.4

0.6

0.8

1.0

t/T

192 193 194 195 196 197 198 199

Fig. 3. Comparison of the lift coefficient from literature and in the present study. To ensure the reliability of the code settings and iteration scheme, a validation study was performed on this model in a previous study [27]. Here, further validation is carried out by comparison to a semi-activated motion model in which the wing experiences a prescribed pitching motion but a free heaving motion. The iteration scheme is similar to that described in 3.2. Deng et al. have performed numerous studies of the semi-activated motion model [29, 30]. The validation results are shown in Fig. 3. The figure illustrates that the results of the model in this study strongly agree with that of Deng et al [29], and it indicates the validity of the code parameter settings.

200

4. Results and discussions

201 202 203 204 205

A previous study focused on the effects of the pitching amplitude and damping coefficient on the response and energy extraction performance of the system [27]. The coupling parameters, R and β, were not considered. This study will investigate the system response at different coupling parameters. Some basic parameters are summarized in Table 1. Table 1 Values of basic parameters used in this study wing NACA 0015 Re ≈ 106 m* 1 I* 0.1 c* 5 θ0 70°

206

4.1 Types of system response

207 208 209 210 211 212 213

The dual wings coupling system has two degrees of freedom. If the flow is regarded as an external excitation, the whole system can be treated as a forced vibration system with two degrees of freedom. Peng et al. considered a flapping foil mounted on a damper and a rotational spring, and found that four different responses were observed, depending on the system configuration and mechanical parameters [31]. The computational results in this study indicate that the system exhibits three distinctive dynamic behaviors depending upon R and β, as shown in Fig. 4. These three response types can also be found in the study conducted by Peng on a single wing [31].

8

ACCEPTED MANUSCRIPT 1.8

R/L= 1, =0.25

R/L= 1, =0.25

R/L= 1, =0.25

y

0.9

0.0

-0.9

-1.8

214 215 216 217 218 219 220

0

20

40

t

(s)

60

80

100

Fig. 4 Time history of heaving motion in Wing 1 At large R and small β (e.g, R/L=1 and β=0.25), the system experiences a typical damped oscillation reduction, in which the heaving amplitude and flapping frequency gradually decrease (response I). The dual wings will eventually stop at the equilibrium position, with a zero pitching angle. When R/L =0.5 and β=0.25, the dual wings undergo an irregular flapping motion (response II). At other values of, R/L and β, the dual wings maintain a periodic pitching and heaving motion (response III). This is the desirable response, although it is not a perfect harmonic motion. t1

1.28

t3

t2

t4

CPF1 CPF2 CPM1 CPM2

0.96 0.64 0.32 0.00 -0.32 -0.64 -0.96 -1.28

0

2

4

t1

2.24

6

t

8

(s)

10

14

16

t4

t3

t2

12

CP1 CP2 CPO

1.68 1.12 0.56 0.00 -0.56 -1.12

221 222 223 224

0

2

4

6

t

8

(s)

10

12

14

16

Fig. 5 Time history of instantaneous power coefficients (response type I, R/L=1, β=0.25) Different values of R appear to result in entirely different response behaviors in responses I and II. However, it would be inappropriate to simply attribute this behavior to the parameter R. To clarify the reason for the observed behavior, the individual contributions of the heaving and pitching motion are examined. 9

ACCEPTED MANUSCRIPT 225

As shown Fig. 5, both 𝐶𝑃𝐹1 and 𝐶𝑃𝐹2 have positive values most of the time. This indicates that

226

the heaving motions are capable of harvesting energy. Conversely, both 𝐶𝑃𝑀1 and 𝐶𝑃𝑀2 are negative.

227 228 229 230 231 232 233 234 235 236 237 238 239

This signifies that the pitching motions consume energy. Interestingly, the power coefficient for the heaving contribution of one wing cannot simultaneously maintain a positive or negative value with the pitching contribution of the other wing. Moreover, the maximum negative value of the pitching contribution generally appears in the vicinity of the zero pitching angle, because a large amount of energy is used for propelling the wing’s pitching motion in that region, while little energy is used to overcome the damper. As a result, the dual wing flapping motion will eventually cease. To understand the reason that the pitching motion consumes a large amount of energy, Fig. 6 shows the near-body flow field around the wings at representative moment where the power coefficient of the pitching contribution is near its maximum. The figure shows the relative position but not the actual distance between the dual wings. Moments t1~t4 are marked in Fig. 5. At moments t1 and t3, the dual wings have the same pitching angle direction, and the direction of angular velocity of the dual wings are opposite. At moments t2 and t4, the directions of the pitching angles are different, while the dual wings have the same angular velocity direction.

240 241 242 243 244 245 246 247 248 249 250 251 252 253 254

(a) t=t1 (b) t=t2 (c) t=t3 (d) t=t4 Fig. 6 Pressure contours around the wings at representative moments At moment t1, Wing 1 has an anti-clockwise pitching motion, whereas Wing 2 is undergoing a clockwise pitching motion. The angular velocities of Wing 1 and 2 are 1.67rad/s and 1.08rad/s, respectively. When fluid flows through a wing which has a nonzero angle of attack, the boundary layer will be changed and form a vortex. A leading edge vortex will detach from the surface of Wing 1, and pitching motion will accelerate the vortex shedding. The vortex creates a deeply negative pressure zone in the wing head. Traditionally, this negative zone will cover the upper surface of Wing 1. A positive pressure zone is observed in the rear of the pitching axis owing to a pitching motion with high angular velocity. Therefore, a larger torque is loaded on the wing. The pitching motion consumes a significant amount of energy because of a larger and opposite torque and angular velocity. Although the pitching angle of Wing 2 is obviously larger than that of Wing 1, a relatively slight negative pressure appears in the wing head. A possible reason for this is that the pitching motion presses the flow and alleviates the boundary layer separation. Consequently, the torque acting on the pitching axis is relatively small, and the lift force loaded on the wing is large. At moment t1, the lift acting on Wing 2 is largely used to compensate for the energy consumption and to sustain the pitching 10

ACCEPTED MANUSCRIPT 255 256 257 258 259 260 261 262

motion of Wing 1. Analogous analyses can be applied to moments t2, t3 and t4. Almost without exception, the maximum energy consumption of the pitching motion occurs in the vicinity of the maximum pitching angle of the other wing. At large R, the lift force drives the pitching motion of the other wing more easily. Thus, the wing pitching motion is driven by the heaving motion of the other wing but does not depend on its own pitching torque. This leads to the opposite directions between the torque and angular velocity at most time. This is the primary reason that the pitching motion consumes a large amount of energy and the dual wings eventually stop at the equilibrium position. CPF1 CPM1

5 4

CPF2 CPM2

t1 t2

3 2 1 0 -1 -2 -3

0

10

CP1

7

20

CP2

t

(s)

30

40

50

40

50

t1 t2

CPO

6 5 4 3 2 1 0 -1 -2 -3

263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278

0

10

20

t

(s)

30

Fig. 7 Time history of the power coefficients (response type II, R/L=0.5, β=0.25) Fig. 7 plots the instantaneous power coefficients of response type II. Similar with the response type I, the instantaneous power coefficient of the pitching contribution has a negative value in most of the time. Energy is needed to sustain the wing pitching motion. Differ from the response type I, the energy consumption of the pitching motion is obviously lower than the energy harvesting of the heaving motion. Moreover, the power coefficient of the pitching contribution frequently appears a large positive value. Thus, the wings can flap continuously. In particular, when both the power coefficient of the heaving contribution and pitching contribution are less than zero simultaneously, the wing’s motion is completely propelled by another wing. At some moment, the energy of harvested by the dual wings cannot meet that required by the motion, which leads to a negative overall power coefficient. Under operating conditions without a control system, the dual wings will sustain their motion depending on their own inertia. Although response type II can achieve a fully flow-induced and self-sustained flapping motion, this condition should be avoided because of its instability and poor energy output performance. As Fig. 7 shows, the power coefficients (both the individual contributions of Wings 1 and 2 and the overall power coefficient) are extremely irregular. These fluctuations are possibly associated with 11

ACCEPTED MANUSCRIPT 279 280 281 282 283 284

a complicated wake consisting of a group vortex near the wing. The force and torque acting on the wing are closely associated with vortex formation, evolution and shedding. If there is no vortex, the pressure center around the wing is near the pitching axis and the torque acting on the wing is relatively small. When a vortex is generated, the value of the torque will change with the vortex evolution, which has an obvious effect on the pressure center around the wing. 3000 0 -3000

Pressure (Pa)

-6000 -9000 -12000

wing 1

-15000

upper surface lower surface

-18000

pitching axis

-21000 -24000 -27000

0.0

0.2

0.4

wing 2 upper surface lower surface

0.6

0.8

1.0

Chord

(a) t=t1 3600 0 -3600

Pressure (Pa)

-7200 -10800

wing 1 -14400

upper surface lower surface

-18000

pitching axis

-21600

wing 2

upper surface lower surface

-25200 -28800

0.0

0.2

0.4

0.6

0.8

1.0

Chord

285 286

(b) t=t2 Fig. 8 Vortex contours and pressure distribution around the wings Table 2 Response parameters corresponding to moments t1 and t2 in Fig. 8 t

287 288 289 290 291 292 293 294

CM1

𝜃1 (rad/s)

CM2

𝜃2 (rad/s)

t1 -0.278 -0.226 -2.070 1.121 t2 1.001 0.437 -2.181 -0.728 Fig. 8 shows the vortex fields and pressure distribution around the wings at moments t1 and t2 in Fig. 7; the corresponding numerical results are summarized in Table 2. For Wing 1, a leading vortex has already detached and is shedding from the wing surface at moment t1. It interferes with the trailing edge of the wing and generates a negative pressure zone. The pressure gradient in the head of the wing’s lower surface returns to a normal level when the vortex detaches. A negative torque acts on the wing. However, the direction of the torque is same as that of the angular velocity. At moment t2, the leading-edge vortex is not formed at the upper surface owing to the large angular velocity. A uniform negative pressure gradient is distributed on the upper surface. 12

ACCEPTED MANUSCRIPT 295 296 297 298 299 300 301 302 303 304 305

Boundary layer separation can be observed in the head of the lower surface because the wing presses the flow. This slightly affects the pressure distribution. There is a large positive torque which has an opposite direction from the angular velocity. Consequently, CPM1 changes sharply from moment t1 to moment t2. Wing 2 has a larger pitching angle and smaller angular velocity. A leading edge is observed, and it creates a deeply negative pressure gradient in the wing head, as illustrated in Fig. 8(a). This negative pressure zone has a slight effect on the torque because its area is very small. As pitching motion continues, the pitching angle approximately reaches its the maximum at moment t2. This reduces the pressure gradient but enlarges the area of the negative pressure zone. Consequently, the torque is improved. The value of CPM2 does not experience a sudden change because the direction of the torque is opposite to that of the angular velocity at moments t1 and t2. CPF1 CPM1

1.25

CPF2 CPM2

1.00 0.75 0.50 0.25 0.00 -0.25

306

0

10

20

30

40

t

50

60

(s)

2.0

CP1

70

CP2

80

90

100

CPO

1.5

1.0

0.5

0.0

-0.5

0

10

20

t

(s)

30

40

50

307 308 309 310 311 312 313 314 315 316

Fig. 9 Time history of the power coefficients (response type III, R/L=0.25, β=2) Fig. 9 shows the variation in the instantaneous power coefficient corresponding to response III. The power coefficient of the heaving contribution is greater than zero most of the time and less than zero for an extremely short time. This is the key that allows the dual wings flap continuously. In Fig. 9(b), the overall power coefficient always has a large positive value. There are no points at which the power coefficients of Wing 1 and Wing 2 are both negative simultaneously. This means that the system is capable of stably outputting power. Of course, the curve of the overall power coefficient is not a smooth flat line. Therefore, to output energy stably, an accumulator is considered in the design of the hydraulic system (the accumulator is not shown in Fig. 1).

317

4.2 Response regions within the parametric space

318 319

To provide guidance for the design of dual flapping wings energy harvesting system, systematic numerical simulations are carried out to determine the coupling system response within the parametric 13

ACCEPTED MANUSCRIPT 320 321

space.

322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339

Table 3 Response types for different parameter combinations β 0.25 0.5 1 2 4 R/L II III III III III 0.25 II III III III III 0.5 I II III III III 1 I II III III III 2 I II III III III 4 Here, the response types are not simply determined according to the heaving (or pitching) displacement of the dual wings but also considering the power coefficient curves. Although in some cases the responses seem to be regular, they are classified as response type II because the overall power coefficient has a negative value. At large R, the lift force of the same component can be converted to a larger torque acting on the pitching motion. Thus, the heaving motion of one wing has a more noticeable effect on the other wing’s pitching motion. The pitching motion of the wings will depend more on the heaving motion, which leads to the torque generally being in the opposite direction of the angular velocity. The pitching motion consumes a large amount of energy because of the higher angular velocity at smaller β values. Consequently, the dual wings experience a reduction in flapping motion at large R and small β values. According to Eq. (9), the angular velocity of the pitching motion will decrease with increasing β. The smaller angular velocity means there is little energy consumption by the pitching motion. Thus, the system has a better response, even at a larger R value. With a decrease in R, the effect of the lift force on the pitching motion of the other wing is weakened, and the effect of torque on the wing’s pitching motion is enhanced. The time when the torque and angular velocity have the same direction is increased. This reduces the energy consumption of the pitching motion. When β=0.25, the response type will vary between I and II at small R values.

340

4.3 Energy extraction performance

341 342 343

This section focuses on response III, in which the dual wings experience periodic heaving and pitching motion and power output is reliable and stable. Table 4 Response indicators for different combinations of R and β R/L β Y/L f* CPF1 CPM1 CPF2 CPM2 CP η (%) 0.5 0.5 0.649 0.096 0.838 0.145 0.838 0.144 1.965 33.78 0.5 1 0.833 0.058 0.710 0.063 0.711 0.063 1.546 24.08 0.5 4 2.408 0.013 0.244 0.001 0.246 0 0.654 4.12 0.25 2 1.456 0.032 0.434 0.015 0.432 0.015 0.896 8.99 1 2 1.196 0.028 0.312 -0.001 0.313 -0.004 0.620 7.65 2 2 1.216 0.028 0.316 -0.003 0.317 -0.004 0.626 7.61 Table 4 summarizes the response indicators for different combinations of R and β. Here, the reduced flapping frequency, f*, is defined as f*=L/TU∞. For response III, Wing 1 has the same flapping frequency and heaving amplitude as Wing 2. To reduce the impact of the fluctuation in the response, all given response indicators represent the averaged value of several flapping periods. From Table 4, little difference in the energy extraction performance is observed between Wing 1 and Wing

344 345 346 347 348

14

ACCEPTED MANUSCRIPT 349 350 351 352 353 354 355 356 357

2 with the same parameters. For β=2, a slight difference is observed in the response indicators for R/L=1 and R/L=2. Although its value is very small, the pitching motion has a negative contribution to the overall power coefficient. Thus, energy is needed to sustain the pitching motion. For R/L=0.25, the heaving amplitude is clearly larger than that for R/L=1 and R/L=2. As discussed above, the lift force of one wing slightly affects the pitching motion of the other wing at small R. The pitching motion depends more on its own torque most of the time. While lift is mainly used for sustaining the heaving motion. Consequently, the system response has a larger heaving amplitude and higher flapping frequency, and the power coefficient of the heaving contribution is greater.  = 0.5

4

=1

=2

=4

3 2

CF1

1 0 -1 -2 -3 -4

50

60

70

t

 = 0.5

1.00

(s)

80

=1

90

=2

100

=4

0.75 0.50

y

1

0.25 0.00

-0.25 -0.50 -0.75 -1.00

358 359 360 361 362 363 364 365 366 367 368 369 370 371 372

50

60

70

t

(s)

80

90

100

Fig. 10 Time history of the lift coefficient and heaving velocity at different β values At a fixed R (e.g., R/L=0.5 in Table 4), an obvious difference can be observed in the system response. The heaving amplitude of the system response increases with increasing β, while the flapping frequency decreases with increasing β. Thus, the system frequency is actually determined by the wings’ pitching motion. The system has a higher flapping frequency at smaller β because of the smaller pitching amplitude. Although the portion of the power coefficient of the pitching contribution is slightly higher at small R, the heaving motion is the main contributor to the overall power coefficient in most cases. Eq. (10) shows that three factors can affect the power coefficient: the magnitude of the heaving velocity, magnitude of the lift force and synchronization between the lift force and heaving velocity. The magnitude of the lift force depends largely on the angle of attack or pitching angle, whereas the heaving velocity is decided by the heaving amplitude and flapping frequency. At small β, although the wings have a small heaving amplitude, the heaving velocity is large owing to the high flapping frequency. The pitching amplitude is also large, resulting in a larger angle of attack during the wings’ movement. As a result, the magnitude of the lift force is large. Moreover, 15

ACCEPTED MANUSCRIPT 373 374 375 376

there is a synchronization between the lift force and heaving velocity, as shown in Fig. 10, which results in a higher power coefficient of the heaving contribution. Increasing β will decrease the magnitude of both the lift and heaving velocity. The power coefficient of the heaving contribution is thus decreased. t1 t2 t3

1.0

 = 0.5

=1

=2

=4

0.8 0.6 0.4

CM1

0.2 0.0 -0.2 -0.4 -0.6 -0.8

50

60

70

t1 t2 t3

1.8

t

(s)

 = 0.5

80

=1

90

=2

100

=4

1.2



0.6 0.0

-0.6 -1.2 -1.8

377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392

50

60

70

t

(s)

80

90

100

Fig. 11 Time history of the indicators related to pitching motion at different β values As shown in Table 4, at small β (e.g., β=0.5 and 1), the pitching motion has a positive contribution to the overall power coefficient. This contribution is too notable to neglect. Similar to the power coefficient of the heaving contribution, the pitching contribution also depends on three factors: the magnitude of torque, magnitude of angular velocity and synchronization of the torque and pitching angular velocity. At large β, the pitching amplitude is relatively small, and boundary layer separation and vortex shedding phenomenon will not occur. The torque is small and changes smoothly. Moreover, the magnitude of the angular velocity is small. As a result, the pitching motion has a limited influence on the overall power coefficient. The pitching amplitude increases quickly with decreasing β, and there is a larger angle of attack during the wings movement. The vortex shedding phenomenon occurs when the angle of attack exceeds the dynamic stall angle, and it has a noticeable effect on the pressure distribution and causes a larger fluctuation of the torque. Additionally, the pitching angular velocity is large at small β. Fig. 11 shows the synchronization of the torque and angular velocity. As a result, the power coefficient of the pitching contribution has a relatively large positive value most of the time.

16

ACCEPTED MANUSCRIPT

393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419

(a) t=t1 (b) t=t2 (c) t=t3 Fig. 12 Pressure contours around the wings at representative moments As discussed in section 4.1, when β=0.25, the pitching motion has a negative contribution to the overall power coefficient. It consumes a large amount of energy and leads to an unfavorable response (response types I and II). Fig. 12 shows a near-body field at representative moments (moments t1, t2 and t3 are marked in Fig. 11) for β=0.5. At moment t1, Wing 1 is located in the vicinity of the maximum pitching angle. A vortex has already detached and is shedding from the wing head, creating a negative pressure zone. The torque is small because the vortex is close to the pitching axis. Although the pressure gradient of the negative pressure zone decreases with further evolution of the vortex, the torque increases because the negative pressure zone moves toward to trailing edge. Moreover, the torque and the angular velocity are in the same direction, as shown in Fig. 11. At moment t3, the vortex generated at moment t1 interacts with the trailing edge. Meanwhile, a new trailing vortex has already formed and attaches closely to the wing surface, creating a deeply negative pressure zone. A large torque acts on the wing because the negative pressure zone is far from the pitching axis. Moreover, the angular velocity is close to its maximum value and is in the same direction as the torque. Consequently, the power coefficient of the pitching contribution has a large positive value at this moment. For Wing 2, a vortex can be found attaching to the wing surface at moment t1. However, the vortex has a slight effect on the torque owing to a large angular velocity. With increasing pitching angle and decreasing angular velocity, the gradient of the negative pressure zone increases. The torque is relatively small owing to the small negative pressure area. In addition, the angular velocity is close to zero. Consequently, the power coefficient of the pitching contribution has a small value. Table 4 indicates a greater energy harvesting performance at small R and β values. For example, when R/L=0.5 and β=0.5, the energy harvesting efficiency is 33.78%. However, it should be pointed out that a lower efficiency does not necessary imply less energy harvested. The efficiency for R/L=0.5 and β=4 is obviously lower than that for R/L=1 and β=2. Nevertheless, the overall power coefficient is larger at R/L=1 and β=2. This observation should be considered for practical applications of this turbine.

420

5. Conclusion

421

Through a hydraulic system coupling dual wings, the completely flow-induced flapping motion 17

ACCEPTED MANUSCRIPT 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440

is achieved. Systematic numerical simulations were performed to determine the response and energy harvesting performance of the coupling system in the parametric space. Depending on the rotary actuator radius R and the volume ratio β between the cylinder and the rotary actuator, three distinguishable behaviors are observed. At large R and small β, the dual wings tend to undergo a damped reduction flapping motion because a large amount of energy is consumed to sustain the pitching motion. The dual wings will eventually stop with a zero pitching angle. Decreasing R or increasing β will decrease the energy consumption of the pitching motion. The system response type will change from the reduction flapping motion to a sustainable flapping motion. Although an irregular response can achieve a fully flow-induced flapping motion, it is unfavorable owing to its unstable energy output. A regular response, essential for stable energy harvesting, is realized over a range of coupling parameters. The heaving amplitude, flapping frequency, and energy harvesting performance indicators are predicted for this response type. Both the heaving amplitude and flapping frequency are closely associated with β, but also depend slightly on R. The heaving amplitude increases while the flapping frequency decreases with increasing β. At large R, the heaving motion of one wing has a more noticeable effect on the pitching motion of the other wing. This results in the angular velocity and torque having opposite directions most of the time. Thus, the pitching motion has a small contribution to the overall power coefficient. A relatively larger portion of power coefficient from the pitching contribution is observed at small R.

441

6. Acknowledgments

442 443 444

The authors gratefully acknowledge the support of the National Natural Science Foundation of China (No.51475270) and the Scientific and Technological Development Project of Shandong Province (2014GGX104014).

445

References

446 447 448 449 450 451 452 453 454 455 456 457 458 459 460

[1] Young J, Lai JCS, Platzer MF. A review of progress and challenges in flapping foil power generation. PROG AEROSP SC. 2014;67:2-28. [2] Xiao Q, Zhu Q. A review on flow energy harvesters based on flapping foils. J FLUID STRUCT. 2014;46:174-91. [3] McKinney W, DeLaurier J. Wingmill: an oscillating-wing windmill. J ENERGY. 1981;5(2):109-15. [4] Jones K, Platzer M. Numerical computation of flapping-wing propulsion and power extraction. AIAA paper. 1997;97:0826. [5] Jones KD, Lindsey K, Platzer M. An investigation of the fluid-structure interaction in an oscillating-wing micro-hydropower generator. WIT Transactions on The Built Environment. 2003;71. [6] Kinsey T, Dumas G. Parametric Study of an Oscillating Airfoil in a Power-Extraction Regime. AIAA J. 2008;46(6):1318-30. [7] Kinsey T, Dumas G, Lalande G, Ruel J, Méhut A, Viarouge P, et al. Prototype testing of a hydrokinetic turbine based on oscillating hydrofoils. RENEW ENERG. 2011;36(6):1710-8. [8] Xiao Q, Liao W. Computational Study of Oscillating Hydrofoil with Different Plunging/Pitching Frequency. Journal of Aero Aqua Bio-mechanisms. 2010;1(1):64-70. 18

ACCEPTED MANUSCRIPT 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504

[9] Lu K, Xie Y, Zhang D, Xie G. Systematic investigation of the flow evolution and energy extraction performance of a flapping-airfoil power generator. Energy. 2015;89:138-47. [10] Lu K, Xie Y, Zhang D. Nonsinusoidal motion effects on energy extraction performance of a flapping foil. RENEW ENERG. 2014;64:283-93. [11] Tian F-B, Young J, Lai JCS. Improving power-extraction efficiency of a flapping plate: From passive deformation to active control. J FLUID STRUCT. 2014;51(Supplement C):384-92. [12] Liu Z, Tian F-B, Young J, Lai JC. Flapping foil power generator performance enhanced with a spring-connected tail. PHYS FLUIDS. 2017;29(12):123601. [13] Kinsey T, Dumas G. Optimal Tandem Configuration for Oscillating-Foils Hydrokinetic Turbine. J FLUID ENG. 2012;134(3):031103-. [14] Kim J, Le TQ, Ko JH, Sitorus PE, Tambunan IH, Kang T. Experimental and numerical study of a dual configuration for a flapping tidal current generator. Bioinspiration & biomimetics. 2015;10(4):046015. [15] Xu J, Sun H, Tan S. Wake vortex interaction effects on energy extraction performance of tandem oscillating hydrofoils. Journal of Mechanical Science and Technology. 2016;30(9):4227-37. [16] Karbasian HR, Esfahani JA, Barati E. Simulation of power extraction from tidal currents by flapping foil hydrokinetic turbines in tandem formation. RENEW ENERG. 2015;81:816-24. [17] Karbasian HR, Esfahani JA, Barati E. The power extraction by flapping foil hydrokinetic turbine in swing arm mode. RENEW ENERG. 2016;88:130-42. [18] Zhan J, Xu B, Wu J, Wu J. Power Extraction Performance of a Semi-activated Flapping Foil in Gusty Flow. Journal of Bionic Engineering. 2017;14(1):99-110. [19] Wang Z, Du L, Zhao J, Sun X. Structural response and energy extraction of a fully passive flapping foil. J FLUID STRUCT. 2017;72(Supplement C):96-113. [20] Zhu Q, Peng Z. Mode coupling and flow energy harvesting by a flapping foil. PHYS FLUIDS. 2009;21(3):033601. [21] Zhu Q, Haase M, Wu CH. Modeling the capacity of a novel flow-energy harvester. APPL MATH MODEL. 2009;33(5):2207-17. [22] Dumas G, Kinsey T, Lemay J, Jean Y, Plourde-Campagna M-A, Lalande G. Oscillating Hydrofoil, Turbine, Propulsive System and Method for Transmitting Energy. U. S2011. [23] Zhu Q. Energy harvesting by a purely passive flapping foil from shear flows. J FLUID STRUCT. 2012;34:157-69. [24] Shimizu E, Isogai K, Obayashi S. Multiobjective design study of a flapping wing power generator. J FLUID ENG. 2008;130(2):021104. [25] Ma P, Yang Z, Wang Y, Liu H, Xie Y. Energy extraction and hydrodynamic behavior analysis by an oscillating hydrofoil device. RENEW ENERG. 2017;113:648-59. [26] Ma P, Wang Y, Xie Y, Huo Z. Effects of time-varying freestream velocity on energy harvesting using an oscillating foil. Ocean Engineering. 2018;153:353-62. [27] Ma P, Wang Y, Xie Y, Zhang J. Analysis of a hydraulic coupling system for dual oscillating foils with a parallel configuration. Energy. 2018;143:273-83. [28] Kinsey T, Dumas G. Computational Fluid Dynamics Analysis of a Hydrokinetic Turbine Based on Oscillating Hydrofoils. J FLUID ENG. 2012;134(2):021104-. [29] Teng L, Deng J, Pan D, Shao X. Effects of non-sinusoidal pitching motion on energy extraction performance of a semi-active flapping foil. RENEW ENERG. 2016;85:810-8. [30] Deng J, Teng L, Pan D, Shao X. Inertial effects of the semi-passive flapping foil on its 19

ACCEPTED MANUSCRIPT 505 506 507

energy extraction efficiency. PHYS FLUIDS. 2015;27(5):053103. [31] Peng Z, Zhu Q. Energy harvesting through flow-induced oscillations of a foil. PHYS FLUIDS. 2009;21(12):123602.

20

ACCEPTED MANUSCRIPT

Highlights 1. A hydraulic system is used to couple dual wings’ pitching motion and heaving motion. 2. The response and energy harvesting of a dual flapping wings system are analyzed. 3. A regular response, essential for stable energy harvesting, is realized. 4. The energy harvesting is closely associated with β but also slightly depends on R.