Numerical investigation of the effects of different parameters on the thrust performance of three dimensional flapping wings

JID:AESCTE AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.1 (1-13)

Aerospace Science and Technology ••• (••••) •••–•••

1

67

Contents lists available at ScienceDirect

68

2 3

Aerospace Science and Technology

4

69 70 71

5

72

6

www.elsevier.com/locate/aescte

7

73

8

74

9

75

10

76

12 13

Numerical investigation of the effects of different parameters on the thrust performance of three dimensional flapping wings

14 16 17 18 19

a

Gong Chunlin , Han Jiakun

b,c

c

a

, Yuan Zongjing , Fang Zhe , Chen Gang

a

Shaanxi Aerospace Flight Vehicle Design Key Laboratory, School of Astronautics, Northwestern Polytechnical University, China b State Key Laboratory of Strength and Vibration for Mechanic Structures, School of Aerospace, Xi’an Jiaotong University, China c Shaanxi Key Laboratory of Service Environment and Control for Flight Vehicles, Xi’an Jiaotong University, China

20

23 24 25 26 27 28 29 30 31 32 33

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 24 May 2018 Received in revised form 26 September 2018 Accepted 18 October 2018 Available online xxxx

After billions of years of evolution, many creatures employing flapping wing in nature tend to have excellent flight capabilities. To understand the bionic wings flow mechanism will be helpful to design high performance underwater vehicles and new conception aircrafts. The geometric parameters, kinematic parameters and flow parameters have great effects on the bionic wings thrust performance. Facing the diverse parameters, it’s very difficult to explore the three-dimensional (3D) bionic flapping wing flow mechanism with traditional numerical simulation method. In this paper, a general largescale parallel solver using Immersed Boundary-Lattice Boltzmann Method (IB-LBM) was developed. The evolution procedures of the 3D flapping wing leading edge vortex and wake flow vortex structures were analyzed in detail. Our study explained the 3D flapping wing thrust performance variation with different wing shapes, aspect ratios and pitch-bias angles of attack. Using Chinese supercomputer TianHeII presents a wide range of possibilities for the further development of parallel IB-LBM, employing tens of millions grids will help us to obtain more complete and accurate 3D flapping wing flow field information. It’s indicated that the obtuse wing has the best thrust performance compared with other sharp wing shapes. With the increase of the aspect ratio, the thrust coefficient of flapping wing increases firstly and then decreases, and with pitch-bias angles of attack increases, the thrust coefficient decreases quickly or even shown resistance phenomenon at the large pitch-bias angel of attack. The discussion of these parameters will provide a theoretical basis for improving flapping-like vehicles propulsive performance. © 2018 Published by Elsevier Masson SAS.

DP

22

Keywords: Lattice-Boltzmann method Immersed boundary method Flapping wing Thrust performance Vortex structure

TE

21

34 35

EC

36 37 38 39 40

CO RR

41 42 43 44

1. Introduction

45 47 48 49 50 51 52 53 54 55 56 57 58 59 60

With the rapid development of aviation technologies and related disciplines, the modern aircraft design is pursuing long distance flight, high maneuverability and low energy performances [1]. In nature, flying and swimming have been a source of inspiration for new vehicle. Fishes, birds and insects employing flapping wing are much more excellent than existing artificial flight and underwater vehicles. For example, the turning radius of anglerfish is only 6% of its body length [2], while the underwater vehicle has to reduce more than 50% of its speed to decrease the turning radius. The creatures’ amazing ability of flying or swimming inspires a few investigations about flapping wing thrust performance [3–5]. In particularly, the flapping wing benefits from using simple translational motions compared to existing conventional propulsion systems, which not only exhibits superior performance in unsteady flow but also is easier to manufacture. Therefore, a lot of new

UN

46

61 62 63 64 65 66

b,c

RO

15

OF

11

E-mail address: [email protected] (C. Gang). https://doi.org/10.1016/j.ast.2018.10.021 1270-9638/© 2018 Published by Elsevier Masson SAS.

77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109

concepts about bio-inspired vehicle have been investigated [6–8], which will help us to create high performance aircraft in the military and civil engineering fields. The actual circumstances of birds, insects and fishes employing flapping motion are very complicated, so that experimental observation, theoretical analysis and numerical simulation are all used together to reveal the mechanism of the bionic movement. In the early of 1930s, like-flapping fish swimming has been experimentally observed systematically [9], and later, Lighthill has proposed the large-amplitude slender theory for fishes movement in the theoretical analysis [10]. Theoretical analysis can give quantitative hydrodynamic results by means of simplified approximation and explains the inherent fluid mechanism of like-flapping fish swimming. But the theoretical analysis can only solve the bionic movement with very simple structures and simple movement parameters. Experimental investigation can obtain force and some flow field information to further reveals the nature of the bionic movements and provide basis test data for numerical method. Lin et al. designed a flapping mechanism to simulate the motion of a bird, the thrust and lift were measured employing wind tunnel

110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 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

OF

5

RO

4

the vorticity is considered. Therefore, the propulsion efficiency is overestimated due to underestimation of wake energy loss. However, for the parameters analysis of the hydrodynamic mechanism for three-dimensional flapping wings, it is very important and useful to capture much more details of the evolution procedure of the unsteady vortex structures in much larger scale simulations. In solving the three-dimensional flapping instability characteristics, the traditional numerical simulation method is used to overcome body-fitted mesh low efficiency and viscous dissipation, and it is urgent to develop on the basis of a large number of grid parallel solving methods. Different from the traditional CFD method using body-fitted mesh to solve the problem of moving boundary, the immersed boundary method (IBM) which was firstly proposed by Peskin [36], modulates the boundary motion into force rather than the boundary condition. Without the time-consuming mesh update procedure at each time step, the immersed boundary method is especially suitable for solving the unsteady flows with large deformation wings. It has been widely used in the fields of animals swimming, biological valve movement, transportation of fluid particle and swing of filament [37], such as the propulsion performance of two-dimensional fin fish [38]. The deformation of the wing involving the moving boundary can improve the performance of the flapping wing MAV. As an alternative of traditional NS equation, the Lattice Boltzmann Method (LBM) is a mesoscopic simulation technique which takes the transport equation as its govern equation. The LBM solvers have the same numerical accuracy as the NS solver under the unremitting efforts of many scholars, which has become an important tool for simulating complex flow in a broad application prospect [39–42]. Since IBM and LBM both use Cartesian grids, it is interesting to combine the two approaches to simulate complex unsteady flow problems. The IB-LBM method was first introduced by Feng and Michaelides where the restoring force caused by the boundary deformation was obtained by punishment force method [43,44], which has good stability for the unsteady flow problems with large deformation boundary and complex shape geometric boundary. It is widely used in bionic motion simulations [45–47]. In this paper, using Chinese supercomputer TianHe-II presents a wide range of possibilities for the further development of parallel IB-LBM, employing tens of millions grids will help us to obtain more complete and accurate 3D flapping wing flow field information at different geometry parameters. The in evolution procedures of the 3D flapping wing leading edge vortex and wake flow vortex structures were analyzed in detail. Our study explained the 3D flapping wing thrust performance variation with different wing shapes, aspect ratios and pitch-bias angles of attack. This paper is organized as follow: Firstly, the IB-LBM numerical method and solver validation is briefly introduced in Sec. 2; and independence verification is obtained in Sec. 3. our problems description and different parameters on thrust performance was numerical investigated in Sec. 4; and Sec. 5 is conclusion.

EC

3

tests [11]. The experimental visualization of the three-dimensional wake structure of flapping wing was achieved using PIV combined with traditional force measurements [12]. However, experiment research for bionic movement has also some disadvantages, such as poor repeatability, lack of enough detail of the flow field and also very expensive. It still is a great challenge to develop more high-resolution and cheaper experimental techniques. Therefore, how to accurately simulate three-dimensional flapping motion using numerical method becomes one of the key issues to biological and biomimetic flow simulation. Numerical simulations can obtain many flow detail information compared to experimental methods. In recent years, numerical simulations are becoming more and more attractive to investigate the parameter effects on the flow filed structure and thrust performance of pitching and heaving flapping wings. The vortex structure and hydrodynamic performance of a pitching and heaving ellipsoidal flapping was investigated by an Navier–Stokes solver [13]. The aerodynamic characteristics of a flapping wing with different pitching elasticity were calculated by solving unsteady Reynolds-averaged Navier–Stokes equations [14]. An immersed-boundary-method-based incompressible Navier–Stokes equation solver was also used to examine the variations in the wake structures and propulsive performance of a pitching-rolling motion plates over a range of major parameters such as Strouhal numbers, Reynolds numbers, and plate aspect ratios [15]. Three-dimensional instabilities in the wake of a plunging and pitching wing of infinite aspect ratio at low Reynolds number were considered by varying the mean pitch angle and the phase shift by two-dimensional and three-dimensional DNS simulation [16]. The effect of wing mass in the free flight of a two-dimensional flapping wing was numerical investigated by an immersed boundary-lattice Boltzmann method [17]. Dash et al. pointed out flapping wing experiences thrust deterioration when the flapping frequency exceeds a certain critical value basing on N–S method [18], which need choose new angle of attack profiles. As it can be seen as the above, different parameters have great effects on the thrust performance and the flow structures of rigid flapping wings, such as the flow parameters [15,16], kinematic parameters [19] and geometry parameters [13]. It is well known that changes in geometric parameters under threedimensional effects have a huge impact on the flapping wing. However, it is necessary to overcome the problem of mesh volume and moving boundary using traditional CFD method to solve three-dimensional problems. In the study of many geometric parameters, three-dimensional phenomena are often ignored. In the past research, two-dimensional flapping wings draws much more attentions [20–24]. Many numerical studies are carried out in which the spanwise variability of the geometry and flow field were ignored [25–27], and large aspect-ratio wings were used for focusing on the two-dimensional effect [28,29]. These twodimensional cases given many coarse explanations on the fluid mechanism and hydrodynamics force of flapping wings. For example, a two-dimensional flapping foil will produce reverse Karman vortex street under the appropriate combination of parameters. For the thrust flapping foils, the reverse Karman vortex street in the trailing area will appear when the Strouhal number is between 0.25 and 0.35. The influence of control parameters, such as the angle of attack and the asymmetric pair also have great effects on propulsion performance [30–32]. A large number of experiments and calculations show that three-dimensional vortex structures are far more complex than those of the two-dimensional foils. The vortex topology induced by three-dimensional flapping wings was obviously different from the reverse Karman vortex streets generated from two-dimensional flapping wings [33–35]. By considering the computational time cost, most of the 2D simulations are based on several millions of meshes. In fact, in the two-dimensional approximation, only the cross-flow component of

CO RR

2

[m5G; v1.246; Prn:24/10/2018; 13:32] P.2 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

UN

1

AID:4801 /FLA

DP

2

TE

JID:AESCTE

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119

2. Numerical method and validation

120 121

2.1. The lattice Boltzmann method

122 123

The lattice Boltzmann equation uses a square or cube mesh to discretize the fluid region into a series of grid points. The fluid is seen as a fluid micelle that can only move toward the surrounding grid point or stand still. The lattice Boltzmann equation describing the motion of fluid particles is as follows [48]:

124 125 126 127 128 129

f i (x + ei δt , t + δt ) − f i (x, t )  3 1 eq = − f i (x, t ) − f i (x, t ) + ωi f · ei δt τ 2

130 131

(1)

132

JID:AESCTE AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.3 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

6 7 8 9 10 11

The corresponding weighting coefficient is:

14 16 17 18 19

22

eq fi



9

25 26 27 28 29 30 31 32

35

40 41

2

ρ=



fi

(6)

47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

0,

u=

i

ei f i + 0.5fδt

(7)

ρ

p

= c 2s

(13)

F(s, t ) = α

ρ

(8)

√

where c s is the lattice sound velocity, c s is equal to 1/ 3 for the D3Q19 model. 2.2. The immersed boundary method

The idea of Immersed Boundary Method (IBM) and traditional Computational Fluid Dynamics (CFD) are different in dealing with moving boundary. IBM considered the boundary as internal force in fluid, and then took it into the Navier–Stokes equations or the Boltzmann equation with external force term. Compared to CFD methods using body-fitted meshes IBM does not need to regenerate new mesh at each time step when solving complex moving boundary problems, which greatly improves the computational efficiency. IBM has been applied in the fields of biomimetic motion and filament swing, and has become a research hotspot in computational fluid dynamics. In order to clearly introduce the concept of IBM, assume that a flexible filament (denoted by Γ ) is immersed in a two-dimensional incompressible viscous fluid (denoted by Ω ). The boundary of the flexible filament Γ is denoted by Lagrangian coordinate s, and the fluid area Ω is represented by Euler coordinates x. The coordinates on the flexible filament can be expressed as x = X(s, t ). The entire systems of the immersed boundary method combined with the incompressible NS equations are listed below [43]:



∂u ρ + u · ∇ u = −∇ p + μ∇ 2 u + f ∂t

74 75 76 77 78 79 80 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98

 









u X(s, t ), t 1 − u dt 1 + β u(X(s, t ), t − u]

99

(14)

where u is the surface velocity of the object. For a given moving object, its speed u(X(s, t ), t ) can be previously interpolated from the velocity information on the nearby Euler points. The feedback coefficients of α and β correspond to speed time integral and speed feedback respectively. For a detailed introduction to the speed correction IB-LBM algorithm, refer to Suzuki’s work [53]. The basic scheme of 3D IB-LBM parallel solver for flapping wing can be linked to calculation process and mathematical formulation through the following Table 1.

102 103 104 105 106 107 108 109 110 111 112

2.4. Solver validation

113

To validate the accuracy and reliability of the IB-LBM solver, several cases are tested. The first case is the steady flow around a sphere and the second case is the unsteady flow past a flapping wing. These validation cases have been explained in our previous articles [54]. Moreover, the unsteady flow past of biomimetic wing with IB-LBM has been studied in our previous works [55]. Here, the unsteady flow past a three-dimensional ellipsoidal flapping wing [13] was used to shown and validate the performance of the solver for bionic wings simulation. Pitching-and-heaving mode has been used widely in the simulation of fish and birds. The center of the wing heaves in the z-direction according to:

Z (t ) = A z sin(2π f )

100 101

0

The pressure of the fluid at the lattice point can be given as:

72

81

|r | ≤ 2x; |r | ≥ 2x.

The concentrated force can be constructed by different methods at a point in IBM, which can be divided into virtual boundary model [50], spring force model [51] and direct force model [52]. All three models can simulate the boundary conditions of the object, but they all have their own scope of application. The virtual boundary model is suitable for processing rigid objects and their active motion. The spring force model and the direct force model are suitable for handling fluid–structure coupled motion. The object of this paper is the numerical simulation of rigid active flapping wing, so the virtual boundary model is chosen. The concentrated force in the virtual boundary model is calculated by:

t

i



UN

46

τ and (5)

44 45

(4)

2

42 43



The velocity and mass of the fluid at the lattice point can be computed by:

37 39

3

υ = (2τ − 1)/6

36 38

1 (1 + cos( 2π|rx| )) 4

2.3. The IB-LBM

The relationship between dimensionless relaxation time fluid viscous coefficient υ is:

33 34

(3)

= ρωi 1 + 3ei · u + (ei · u) − u · u

23 24

δ(r ) =

2

71

(12)

73

The equilibrium state distribution function of D3Q19 model is:

20 21



⎧ ⎨ 1 /3 i = 0 ωi = 1/18 i = 1, 2, · · · , 6 ⎩ 1/36 i = 7, 8, · · · , 18

CO RR

15

 u(x, t )δ x − X(s, t ) dx

where p is fluid pressure, ρ is fluid density, u is fluid viscosity, F(s, t ) is concentrated force density at immersed boundary. f(x, t ) is fluid force density. δ(x − X(s, t )) is the Dirac delta interpolation function. The three-dimensional Dirac delta interpolation function in equation (12) is as follows [49]:

(2)

13

69 70



Ω

⎧ i=0 ⎨ (0, 0, 0) i = 1, 2, · · · , 6 ei = (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) ⎩ (±1, ±1, 0), (0, ±1, ±1), (±1, 0, ±1) i = 7, 8, · · · , 18

12

∂X = ∂t

Γ

OF

5

67 68

(11)

RO

4

(10)

DP

3

∇ ·u=0

  f(x, t ) = F(s, t )δ x − X(s, t ) ds

TE

2

eq

where f i is the particle density distribution function, f i is the corresponding equilibrium state distribution function, τ is the dimensionless relaxation time, δt is the time step, f is the fluid density, ei is the particle velocity, ωi is the weight coefficient, ei and ωi are related to discrete density model. In this paper, D3Q19 model is chosen. In D3Q19 model, the velocity is modeled as:

EC

1

3

(15)

At the same time, the wing pitches about its center according

114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130

to:

131

(9)

θ(t ) = α0 + A θ sin(2π f t + ϕ )

(16)

132

JID:AESCTE

AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.4 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

4

1

67

Table 1 Correspondence between inhouse code and mathematical formulation.

2 3

The solver code calculation process

4

Initialize the fluid field Calculate the concentrated force Disperse the concentrated force on the Euler grid point Get new velocity field

5 6 7 8

68 69

Mathematical formulation eq f i (x, t )] + 32

1

f i (x + ei δt , t + δt ) − f i (x, t ) = − τ [ f i (x, t ) − t F(s, t ) = α 0 [u(X(s, t ), t 1 ) − u]dt 1 + β[u(X(s, t ), t ] − u] f(x, t ) = Γ F(s, t )δ(x − X(s, t ))ds

ωi f · ei δ t

73

eq

f i (x + ei δt , t + δt ) − f i (x, t ) = − τ1 [ f i (x, t ) − f i (x, t )] +

3 2

ωi f · ei δ t

76

Table 2 The calculation parameters of the verification example.

12 13 14 15

77

The heaving amplitude A z

The pitching amplitude A θ

The pitch-bias angle α0

The phase shift ϕ

St number

Re number

0.5

π /6

0

π /2

0.6

200

16 17

RO

18 19 20 21 22 23

DP

24 25 26 27 28 29 31 32

Fig. 1. Comparison of verification results. The hydrodynamics coefficients of verification example and our results. (a) Lift coefficient with combined parameter of frequency and time. (b) Thrust coefficient with combined parameter of frequency and time.

35

39 40

where A z and A θ represent the heaving and pitching amplitudes respectively and f is the flapping frequency. ϕ is the phase shift between heaving and pitching. α0 is the pitch-bias angle. The Strouhal number (St) and Reynolds number (Re) are defined as, respectively:

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

Re = U ∞ L /υ

(17)

where Stream speed U ∞ is set as 1, υ is fluid viscous coefficient. The thrust and lift coefficients are defined as, respectively:

CT =

T 2 A 0.5ρ U ∞

CL =

L 2 A 0.5ρ U ∞

(18)

where T is the force in the wing span direction, and L is the force in the heaving direction. A  denotes twice the area of the ellipsoidal wing. The calculation parameters of the verification example are shown in Table 2. All of the above parameters are dimensionless in the current study. As shown in Fig. 1, the asterisk scatter chart shows the result of Dong; the solid line is the calculation result of this article, the numerical results of the thrust and lift coefficients were good agreement with Dong’s who used the finite volume NS equations with immersed boundary method [13]. Therefore, provides a powerful platform for solving the problem of flapping wing unsteady flow mechanism and is broadly applicable to a variety of parameter research.

UN

43

St = 2 A z f /U ∞

CO RR

41 42

62 63

3. Independence validation

64 65 66

three-dimensional NACA0012 wing at a constant inflow, the mesh and time step independence are verified in detail. Of course, each of the cases given in this paper is verified in their mesh and time step independence, but it’s not listed because of the paper length.

EC

38

TE

30

37

OF

11

36

74 75

10

34

71 72

9

33

70

To illustrate the independence of the selected mesh and time step, here we consider the pitching-heaving flapping motions of a

3.1. Mesh independence validation

The pitching-heaving flapping motions of a three-dimensional NACA0012 wing at a constant inflow U ∞ as shown in Fig. 2 (a). The direction of incoming fluid coincides with X axis and Y axis is spanwise direction of the wing. The length of the wing in y direction is 2 and the one in x direction is 1. The maximum thickness of the NACA0012 wing is 0.12. The specific movement mode can be described in equation (15) and (16). The focus of our study is to capture more flow field information using a large number of an overall-grid instead of local hybrid grid. The flow field calculation domain information is given in Fig. 2 (b), whose size is 18 × 10 × 10. This choice of domain size was based on our experience with running test simulations on a number of different domains. The total mesh of fluid domain is 540 × 300 × 300, which is 48.6 million. It was treated as nominal grid. The overall grid size for this refined gird is 630 × 350 × 350 which amounts to about 77.1 million grid points. The overall gird size for this roughed gird is 450 × 250 × 250 which amounts to about 28.1 million. The time step dt equals 0.001. Fig. 3 shows a comparison of key dynamic quantities for the grid independence study carried out for the St = 0.6, Re = 200, ϕ = π /2 and A θ = π /6 case. To give a clearer explanation of mesh independence, a quantitative analysis was introduced. As shown in

78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.5 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

5

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

OF

1

13 14

Fig. 2. (a) is NACA0012 flapping wing model, (b) is the flow field calculation domain. The flow field employed a standard Cartesian grid.

15 16 17

RO

18 19 20 21 22 23

DP

24 25 26 27 28 29 30

33 34 35 36 37 38 39 40

Fig. 3. The hydrodynamics coefficients of flapping wings at different grid numbers for St = 0.6, Re = 200, force coefficient.

Coarse Nominal Finer

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

3.2. Time independence validation

86 87 88 89 90

97 98 99 100 101 102

0.001

28.1 48.6 77.1

0.4497 0.4555 0.4594

0.86 – 1.3

The same calculation model and kinetic parameters as in the previous section is chosen. To verify the time independence, we use identical normal grid but different time steps is dt = 0.0005, dt = 0.001 and dt = 0.0015, named as compact time, normal time and loose time respectively. Fig. 5 shows the lift and thrust force coefficients under different time steps. From this figure, it can be seen that lift and thrust force coefficients are extremely similar with qualitative analysis. We also give the average thrust coefficient as a comparison to verify the time independence, as shown in Table 4. The relative error is given, and the independence of

85

96

The relative errors (%)

Table 3, the average thrust coefficient is compared to verify the mesh independence. Compared with nominal mesh, the relative errors between the refined and coarse meshes are about 0.86%, 1.3%, respectively. Meanwhile, it is shown that the lift and thrust force coefficients didn’t obvious change with the number of grid changed by comparing Fig. 3. Fig. 4 shows the comparison of flow field vorticity cloud diagrams under the three grid states. It is shown that nominal grid and finer grid are basically consistent on the vorticity, which indicated the solver has a good convergence in the grid selection. The selected computational domain and the computational grid are justified.

84

95

The thrust coefficient C t

UN

43

83

94

Number of overall-grid (×106 )

CO RR

42

82

93

Time step dt

41

81

92

Table 3 Mesh independence validation: the average thrust coefficient comparison. Mesh

80

91

ϕ = π /2 and A θ = π /6 case. (a) Lift force coefficient. (b) Trust

EC

32

TE

31

79

103 104 105 106 107 108

time step is further explained from the quantitative point of view. Fig. 6 shows the comparison of flow field vorticity cloud diagrams under the three time steps. We can find that the obtained flow field information is the same as the problem described by the kinetic curve, which is shown that the solver has a good convergence in the time selection.

109 110 111 112 113 114 115

4. Results and discussions

116 117

4.1. Problem statement

118 119

Our focus is on using the parallel solver that has been developed to explore the 3D flapping thrust performance and flow mechanism under different parameters. In this paper, the main parameters we consider include: the wing shape, the aspect ratio and pitch-bias angel of attack. The motion mode of the flapping wing adopts the above described conventional pitching and heaving. After the mesh and time step independence are verified, the same settings as in the previous section verification example are used. According to our experience, the calculation domain size still chooses 18 × 10 × 10. The parameters that need to be changed are given in Table 5. In this paper, we don’t consider the influence of kinematic parameters for the time being, so the other parameters are consistent with Table 3 in the solver verification case.

120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.6 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

6

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76 77

11

Fig. 4. The vorticity contour plot of NACA0012 flapping wing spread spanwise symmetry plane under three different grids. (a) Coarse grid, (b) Nominal grid, (c) Finer grid.

OF

12 13 14 15 16 17

RO

18 19 20 21 22 23

DP

24 25 26 27 28 29 31

Fig. 5. The hydrodynamics coefficients of flapping wings at different time steps for St = 0.6, Re = 200, coefficient.

32 34 35

Mesh

Time step dt

Nominal

0.0005 0.0010 0.0015

36 37 38

Number of overall-grid (×10 ) 6

48.6

39 40

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

4.2. The influence of wing shape on thrust performance

The momentum theorem is usually used to analysis the bionic wing’s performance [48]. As shown in Fig. 7, there are symmetry jets about the x– y plane in the wake regions by giving the one cycle average velocity profile. It’s indicated that the average lift coefficient of is zero, which is also the characteristic that all symmetrical wing should exhibit. The average velocity in the wake zone present a jet shape, which is consistent with the other researchers [56,57] have also shown the jet-like velocity profile is a public characteristics related to the flapping wing movement. The above results further verify that our solver performs well in solving

82 83 84 85 86 87 88 89 90

94 95 96 97 98 99 100

The relative errors (%) 1.6 – 0.08

Fig. 6. The vorticity contour plot of NACA0012 flapping wing spread spanwise symmetry plane under three different time steps. (a) Loose time, (b) Nominal time, (c) Compact time.

UN

42

0.4630 0.4555 0.4551

81

93

CO RR

41

The thrust coefficient C t

80

92

Table 4 Time step independence validation: the average thrust coefficient comparison.

EC

33

79

91

ϕ = π /2 and A θ = π /6 case. (a) Lift force coefficient, (b) Trust force

TE

30

78

the 3D flapping wing. According to momentum theorem for incompressible flow, thrust can be obtained by the integration of the velocity square distribution along the control surface times density. The computing domain we choose is the same. Because of the velocity magnitudes at the outlet surface for the three cases are very similar, the thrusts are very alike too. The jet region average velocity strength in Fig. 8 (a), (b) and (c) are obviously different. The mean velocity strength distribution of the cuboid wing shown as in Fig. 7 (c) is significantly smaller than those of the other two airfoils. Although the thrust variation between the three is very similar according to the momentum theorem, the thrust of the cuboid wing thrust should be smaller from the velocity strength.

101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.7 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3

Table 5 Different parameters of the flapping wing that need to be studied. The wing shape

The aspect ratio

4 5

The pitch-bias angle of attack

1. The influence of wing shape on thrust performance

6 7

Ellipsoidal wing

8 9 10

NACA0012 wing

AR = 2

α0 = 0◦

11

Cuboid wing

13 15 16 17 18 19 20 21 22 23

2. The influence of aspect ratio on thrust performance NACA0012 AR = 1 AR = 2 AR = 3 AR = 5 AR = 7

α0 = 0◦

RO

14

3. The influence of pitch-bias angle of attack on thrust performance α0 = 0◦ NACA0012 AR = 3 α0 = 10◦ α0 = 20◦ α0 = 30◦

DP

24

28 29 30 31 32 33 34 35 36 37

EC

27

Therefore, we give quantitative analysis from the perspective of the rising resistance coefficient. Fig. 8 (a) and (b) show the relationship between the lift and thrust coefficients of the above three flapping wings and the f t combination parameters. It can be seen from Fig. 8 (a) that the amplitude of the three case lift coefficients is different. However, they are all symmetric about the x-axis, which lead to the average lift is zero. It is consistent with the expected zero-lift force according to the symmetric average velocity profile shown in Fig. 7. However, the time-related thrust coefficient curve shown in Fig. 8 are asymmetric about the 0 axis, thus the periodic thrust is generated. The average thrust coefficients of the ellipsoidal wing, the NACA0012

38 39 40

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

Fig. 7. Average velocity contour in one flapping cycle at the spanwise symmetry plane under different wing shapes. (a) Ellipsoidal wing, (b) NACA0012 wing, (c) Cuboid wing. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

UN

43

CO RR

41 42

TE

25 26

wing and the cuboid wing are 0.4809, 0.4587 and 0.1125, respectively. The average thrust coefficient of the ellipsoidal wing is the biggest, followed by NACA0012 wing. The thrust coefficient of the cuboid wing is the smallest which should be related to the sharp leading edge and non-streamlined cross-sectional profile. It can be indicated from this that the obtuse wing has the best thrust performance compared with other sharp wing shapes. The values of the average thrust coefficients are consistent with the strength of the average velocity distribution in the wake regions in Fig. 7. The projection configuration of the NACA0012 wing and cuboid wing are the same rectangular. However, their average thrust coefficients are different and are smaller than that of the ellipsoidal wing. It indicates that the thrust coefficient of the flapping wing with ellipse projection shape is larger than those with rectangular projection shape. It is in fact that most of the plane projection shape of fish body and insect wings are similar to the oval, which shows that biological evolution obeyed some profound natural laws. The NACA0012 wing and the plate wing have the same projection shape, however their thrust coefficients are obvious different. This indicates that the streamlined flapping wings can improve their thrust performance as much as possible. This streamlined wing shape is easier to manufacture than the elliptical wing, under the premise of ensuring similar aerodynamic characteristics, the research on the influence of parameter changes on its thrust performance is more instructive for the design of new actual aircraft. From the above results, we can only obtain the aerodynamic characteristics of different wing shape, and why it is necessary to further analyze from the perspective of vorticity. In this paper, the leading edge vortex of different shape flapping wings are compared to reveal the mechanism of force produce. Fig. 9 shows the flapping wing span symmetrical plane vorticity diagram of at entire computational domain. We selected transient clouds at the same position and the same time to compare the differences between the three cases. The flapping motion produces a vortex street with the same jet velocity, but the cuboid wing vortex strength is obviously smaller. In order to explain the influence of different wing shapes on the flapping wing thrust performance

OF

12

7

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125

60

126

61

127

62

128

63

129 130

64 65 66

Fig. 8. The hydrodynamics coefficients of flapping wings at different at different wing shape. (a) Lift force coefficient, (b) Thrust force coefficient. Both represent the relationship between thrust coefficients and combined parameter of frequency and time.

131 132

JID:AESCTE

8

AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.8 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

12

Fig. 9. The vorticity contour plot of NACA0012 flapping wing spread wing span symmetry plane at different wing shapes. (a) Ellipsoidal wing, (b) NACA0012 wing, (c) Cuboid wing. The three states used the same color space.

OF

11 13 14 15 16 17

RO

18 19 20 21 22 23

DP

24 25 26 27

Fig. 10. Partial magnification of the vorticity contour in the wing span symmetry plane. Where, three-dimensional flapping wing is black. The three states used the same color space. (a) Ellipsoidal wing, (b) NACA0012 wing, (c) Cuboid wing.

28 29 30

TE

31 32 33 34 35

EC

36 37 38 39 40

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

Fig. 11. The vortices varieties of the NACA0012 wing movement in one cycle at different time.

from the perspective of flow mechanism, the role of the leading edge vortex is discussed. As shown in Fig. 10, the vorticity cloud map is partially enlarged to clearly observe the change of the leading edge vortex. It can be seen from Fig. 10 (a) that the leadingedge vortex V1 of the ellipsoidal wing is closer to the surface than the other wings. The ellipsoidal wing leading edge vortex is not elongated as quickly as the other two cases, and V1 is fuller and stronger. A leading edge vortex that is closer to the flapping wing surface and has stronger vortex strength results in greater thrust. From Fig. 10 (b) and (c), it can be seen that the leading-edge vortex V1 of NACA0012 wing is closer to the wing upper surface than the cuboid wing. The thrust coefficient of NACA0012 wing is also slightly larger than that of the cuboid wing. Therefore, thrust coefficient of the cuboid wing is smallest actually.

UN

42

CO RR

41

These results indicate that the leading edge vortex V1 is an important determinant mechanism of the thrust performance of bionic flapping wings. The blunt leading edge of the cuboid wing ultimately affects the formation of the leading edge vortex, which will determine the average velocity field behind the airfoil. However, the leading edge of cuboid wing and NACA0012 wing are blunt, causing similar vortex and velocity distribution. Their thrust coefficients are also close to each other. In summary, the leading edge shape is one of the factors affecting the thrust magnitude. It suggests that in bionic flapping wing design, the blunt leading edge wing should be avoided in order to obtain a higher propulsion performance. To further explore the flow mechanism of the flapping wing in a motion cycle, Fig. 11 shows the NACA0012 wing vorticity evolution in one cycle at five different moments. When the wing begins to do the pitching and heaving movement

77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.9 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

9

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76 77

11 13

Fig. 12. The vortex structure of the flapping wing identified by Q-criterion (Q = 0.1) in the wake zone of flapping wings at different wing shapes. Where, the picture marked the same order vortex structure with thin red lines.

OF

12 14 15 16 17

RO

18 19 20 21 22 23

DP

24 25 26 27 28 29 30

TE

31 32 33 35

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

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

(Fig. 11 (a)), the leading edge vortex begins to move across the wing to the trailing edge and meets the vortex shedding in the previous cycle. When the wing is moved up to the intermediate position, the vortex begins to adhere to the lower surface as shown in Fig. 11(b). The strength of the vortex on the wing lower surface is significantly greater than that of wing upper surface vortex, which will help lift the wing. When the wing reaches the highest point as shown in Figure (c), the strength of the vortex on the wing lower surface is significantly reduced, although the vortex attached to the lower surface will affect its movement, it can be seen that the upper surface will produce a stronger vortex is shown in Figure (d). The strength of the vortices on the upper and lower surfaces produced is not the same when the wing is flapping, which will help to generate thrust. From Fig. 11 (d) to (e), the leading edge vortex is always attached to the upper surface during the wing down, resulting in a downward thrust that causes the wing to move downward. As the wing moves downward, the strength of the lower surface vortex becomes smaller, and the upper surface rapidly produces a stronger vortex. The lower surface vortex that hinders the downward movement gradually disappears, which is also the cause of the cyclical change in thrust. Furthermore, the wake flow structure can also indicate many flow mechanism of the 3D flapping wing. Fig. 12 shows the vortex structure identified by Q-criterion (Q = 0.1) in the wake of the three different shape flapping wings. As can be seen from Fig. 12 (a), the vortex rings of the ellipsoid wing gradually decrease compared to the other wings. However, their shapes keep complete vortex wing when they are evolving downstream. They can induce stronger velocity so the thrust performance of the ellipsoidal wing rank top among the three flapping wings. Fig. 12 (c) shows the

CO RR

39

81

101

UN

38

80

Fig. 13. One-cycle average velocity contours of NACA0012 wing with different aspect ratios.

36 37

79

EC

34

78

vortex rings intertwine and affect each other when they leave the cuboid flapping wing. Then the vortex rings fracture, resulting in a rapid reduction of velocity in the wake of cuboid flapping wing, which leads to the worst thrust performance among the four flapping wing. In order to see this phenomenon more clearly, we draw the same order of vortex structure with red lines in the figure. It can be seen from the figure that, the elliptical wing is similar to the NACA0012 wing vortex structure at the red line position, but the vortex structure of the cuboid wing will begin to rupture, and its induction speed will be weak at the same position. As a result, the horizontal velocity component becomes smaller, and the thrust is of course the smallest of the three. We can also see that the vortex structure of NACA0012 and the vortex structure of the elliptical wing are very similar, which inspires us to choose obtuse leading edge wing shape instead of a circle considering the process difficulty in the actual processing and manufacturing. This is also the reason why people are beginning to pay attention to the NACA0012 flapping wing.

103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121

4.3. The influence of aspect ratio on thrust performance

122 123

In this section, the vortex structures and hydrodynamic coefficients of flapping wings with different aspect ratio were investigated. The geometric model is rectangle wing with NACA0012 wing. The length of the flapping wing along x-axis is constant 1 and the lengths along y-axis are 1, 2, 3, 5 and 7, which means the aspect ratio AR is also 1, 2, 3, 5, 7, respectively. The other computational parameters were as same as those in Table 2. The one-cycle average velocity contour of the five flapping NACA0012 wings with different aspect ratios are shown in Fig. 13. With the variation of

124 125 126 127 128 129 130 131 132

JID:AESCTE

10

AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.10 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

OF

1

13 14 15 16 18

Fig. 14. The hydrodynamics coefficients of flapping wings at different at different aspect ratios. (a) Lift force coefficient. (b) Thrust force coefficient. Both represent the relationship between thrust coefficients and combined parameter of frequency and time.

RO

17 19 20 21 22 23

DP

24 25 26 27 28 29 30

TE

31 32 33 34 35

EC

36 37 38 39 40

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

Fig. 15. Partial magnification of the vorticity contour in the wing span symmetry plane at different aspect ratios. The five states used the same color space.

the aspect ratio, the average velocity contours also change obviously. With the increase of the aspect ratio, the velocity contours in Fig. 13 (a) and Fig. 13 (b) are similar with bifurcated jets, the bifurcated jets begin to blend as a main jet in Fig. 13 (c), and then the wake jets disappear as shown in Fig. 13 (d) and Fig. 13 (e). The reason for this phenomenon is that the three-dimensional effect is not obvious because of the increase in the aspect ratio. It can be inferred that the thrust in those cases is small. Although the jet structures are similar for the AR1 wing and AR2 wing, the velocity intensity of AR2 wing is stronger than that of AR1 wing. The thrust coefficient of the AR2 wing is expected to be bigger than that of AR1 wing. The thrust is judged by the average velocity. AR2 should generate the maximum thrust because the average velocity in the wake is higher than the other four cases. However, from the average velocity shape, the main jet of AR3 will also affect the magnitude of the thrust. The velocity intensity in both side of the main jet of the AR3 wing is higher than those near the center line. For the AR5 and AR7 wing, the wake flows still appear as a main jet, however the angle and the velocity intensity of the jets are less than that of AR3 flapping wing. However, further quantitative analysis is still necessary to identify the rank of the thrust coefficients.

UN

42

CO RR

41

Fig. 14 gives the hydrodynamic coefficients at different aspect ratios. The average lift coefficients of five wings with different aspect ratios are 0. The average thrust coefficients of five flapping wings with the aspect ratios of 1, 2, 3, 5 and 7 are 0.3426, 0.4587, 0.5098, 0.3852 and 0.1997 respectively. When the aspect ratio changes from 1 to 3, the average thrust coefficient increases. This is consistent with strongest average velocity contour when the aspect ratio is 3. As the aspect ratio changes from 3 to 7, the thrust coefficients decreases. The thrust coefficient quantitative analysis can be found that the average velocity field prediction results are basically consistent with the actual thrust distribution. By comparing the average velocity field distribution, we can find that smaller angle of the jet in wake flow will produce small velocity, which also induces the decreased thrust performance Furthermore, we analyze the vortex structure of several special circumstances to reveal the mechanism of force produced. Fig. 15 shows the vorticity contours on the wing span symmetry plane of different NACA0012 flapping wings. It is found that when the aspect ratio is 1, the leading-edge vortex is very small and far from the wing surface. As the aspect ratio increases, the shape of the leading edge vortex becomes more and more full and the strength is also enhanced. When the aspect ratio continues to 5 and 7,

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.11 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

11

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

OF

1

13 14 15 16 17

RO

18 19 20 21 22

Fig. 16. The vortex structure of the flapping wing identified by Q-criterion (Q = 0.1) in the wake zone of flapping wings at different aspect ratios.

23

DP

24 25 26 27 28 29

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

the strength of the leading edge vortex begins to become smaller, which is also the cause of the thrust change. As for the tailing-edge vortex, the tailing-edge vortex is relatively slender and far away from the AR1 wing. When the aspect ratio is 2, the tailing-edge vortex is approximately similar to the that of AR1 wing. When the aspect ratio continues to 5 and 7, the strength of the tailing-edge vortex produced on the upper surfaces are stronger. In summary, comparing the thrust coefficients, we can find that different aspect ratios will produce different leading edge vortices and trailing-edge vortices. The strength of the vortex is the main factor affecting the thrust generation. The reason why the thrust of the AR3 is the largest is not only because the leading edge vortex is closer to the wing surface wing and the vortex shape is full, but also because of the stronger leading edge vortex and the weaker trailing edge vortex. Fig. 16 shows the vortex structure identified by Q-criterion when the NACA0012 flapping wings are at the lowest position in their flapping trajectory. It can be seen that the spanning distance of the vortex structure in wing surface is close to the length of the wing. This means that the spanwise length of the vortex structure on the wing surface become larger with the aspect ratio increases. However, the vortices induced by the flapping wings at different aspect ratios are basically the same in the direction of the airflow direction. During the evolution of the vortex ring to the downstream, the span distance is reduced. When the aspect ratio is 1 and 2, the main vortex induced by the wing surfaces is narrow and also the strength is relatively small. Influenced by the leadingedge vortex, the vortex structure becomes vortex wing along the downstream. When the aspect ratio is 3, 5 and 7, the main vortex become longer and longer. Under the same Q value, it can be seen from the fullness of the vortex that the strength of the leading edge vortex is weaker than the strength of the main vortex, which fails enclose of vortex structures. In the downstream wake flow,

the vortex takes on intertwined vortex structures. The wingtip vortex will induce disturbed velocity along its normal direction, which makes vortex structure go away from each other. This phenomenon coincides with those in Fig. 19 where the angle of wake vortex decrease as the aspect ratio increases.

EC

33

Fig. 17. Velocity contour in the spanwise symmetry plane at different pitch-bias angles.

CO RR

32

UN

31

TE

30

4.4. The influence of pitch-bias angle of attack on thrust performance

In the previous cases, the flapping wing pitch-bias angle of attack was assumed to be 0◦ . In fact, the wing pitch-bias angle of attack was also crucial for the propulsion performance of the 3D flapping wings. In this section, the effects of the pitch-bias angles on hydrodynamic coefficients and vortex structure are investigated for NACA0012 flapping wing. Fig. 17 shows the one-cycle average velocity. It can be seen that with the increase of the pitchbias angle of attack, the average speed of the intensity gradually decreased. According to the momentum theorem, the thrust coefficient will decrease as the angle increase. At the same time, for the angle of attack of 0◦ , the average velocity is symmetrical about the x– y plane. However, the average velocity contours become non-symmetrical bifurcated jets, which makes the wing produce a certain lift. It can be seen from Fig. 18 (a) that the lift coefficient curves at four different pitch-bias angles of attack change greatly. The mean lift coefficients for 0◦ , 10◦ , 20◦ , and 30◦ are 0, −0.7007, −1.2762, −1.6797, respectively. This is consistent with the average velocity distribution in the wake flow as shown in Fig. ??. Fig. 18 (b) shows the thrust coefficient time response. For the pitch-bias angle of 0◦ and 10◦ , the forward amplitude is significantly larger than the backward amplitude which lead to the positive thrust force. While for the pitch-bias angle of 20◦ and 30◦ , the negative amplitude is larger than their positive counterparts, so the drag force is generated. The average thrust coefficients for each pitch-bias angle

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:4801 /FLA

12

[m5G; v1.246; Prn:24/10/2018; 13:32] P.12 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

OF

1

13 14 15 16 18

Fig. 18. The hydrodynamics coefficients of flapping wings at different at different pitch-bias angles. (a) Lift force coefficient. (b) Thrust force coefficient. Both represent the relationship between thrust coefficients and combined parameter of frequency and time.

RO

17 19 20 21 22 23

DP

24 25 26 27 28 29

Fig. 19. Vorticity contour in wing span symmetry plane at different pitch-bias angles.

30

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

TE

34

EC

33

tex is full and close to the wing surface. The vortex structures of NACA0012 wing and ellipsoidal wing are very similar, which causes also the similar thrust performance. The vortex structure of the cuboid wing is small and far from the wing which leads to the worst thrust performance. With the increase of aspect ratio, the average thrust coefficient of NACA0012 flapping wings increases firstly and then decreases. The two-dimensional characteristics of the vorticity distribution is obviously exhibited with the increase of the aspect ratio. It means that for the small aspect ratio flapping wing, its three-dimensional effects in thrust performance must be taking into accounting for. With the pitch-bias angle of attack increases, the thrust coefficient of NACA0012 flapping wing decrease, or even shown resistance phenomenon at the large pitch-bias angel of attack. The evolution procedures of the leading edge vortex and wake flow vortex structures are the possible main reasons which could be used to explain the thrust performance changes with parameter variation of the flapping wing. In future work, much more fluid parameters such as Reynolds numbers, Strouhal number, and more different motion patterns such as amplitude, frequency, would be further investigated. At the same time, we will further explore the influence of 3D thickness and leading edge curvature on thrust performance. The effect of the structural flexibility is also very important and obvious for true fish, insects and birds. Taking into account for the fluid-structural interaction effects based on IB-LBM method is also expected in next step, although it is a challenge task in large scale simulation.

5. Conclusions

CO RR

32

can be obtained by averaging the thrust curves in various cases: 0.5098, 0.3714, −0.0104, −0.5055. It can be seen that for these flapping wings, the thrust coefficient will decrease or even become to resistance phenomenon with the pitch-bias angle of attack increases. In order to make the flow mechanism more clearly, the vorticity contours on the symmetry plane of the flapping wings in different pitch-bias angle were presented in Fig. ??. With the increase of the pitch-bias angle of attack, the strength of the leading edge vortex hardly changes much, but the strength of the trailing-edge vortex increases significantly, but the cover distance above the upper wing surface is gradually reduced as shown in Fig. 19 (a) and Fig. 19 (b). Furthermore, with the increase of the pitch-bias angle of attack, the tail vortex is no longer keeping the elongated shape, but gradually moves close to the wing surface. The strength of the trailing edge vortex increases, which leads to the smaller thrust force. The selection of the pitch-bias angle of attack is very important for the bionic flapping wing movement design which has significant effects on the lift and thrust performance.

The influence of geometrical and motion parameters on the flow mechanism and thrust performance of bionic flapping wing is one of the key issues in the design of wing shape and motion control law for bionic flapping aircraft and underwater vehicles. To capture the flow mechanism and fine vortex structures of three-dimensional complex flapping wings is a time-consuming challenge task, especially for the parameter influence investigations in large-scale calculation. In this paper, a three-dimensional parallel IB-LBM solver suitable for running on TianHe-II Supercomputer was developed and validated. The wing shape, aspect ratio and pitch-bias angle of attack have great influences on the thrust performance of the flapping wings. The numerical results show the ellipsoidal wing has the best thrust performance in the three flapping wings, whose leading-edge vor-

UN

31

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123

Conflict of interest statement

124 125

None declared.

126 127

Acknowledgements

128 129

This work was supported by The National Defense Fundamental Research (Grant No. JCKY2016204B102); Fundamental Research Funds for the Central Universities (Grant No. G2016KY0302);

130 131 132

JID:AESCTE AID:4801 /FLA

[m5G; v1.246; Prn:24/10/2018; 13:32] P.13 (1-13)

G. Chunlin et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6

National Natural Science Foundation of China (No. 11672225, 11872293), and the Program of Introducing Talents of Discipline to Universities (known as the “111” Program (B18040)). The authors will also acknowledge the Joint Natural Science Foundation of China with Guangdong Province for TianHe-II supercomputer resources support.

[28] Mahmud Ashraf, J.C.S. Lai, J. Young, Numerical analysis of flapping wing aerodynamics, in: 16th Australasian Fluid Mechanics Conference (AFMC) School of Engineering, The University of Queensland, 2007, pp. 1283–1290. [29] Kosuke Suzuki, M. Yoshino, Aerodynamic comparison of a butterfly-like flapping wing-body model and a revolving-wing model, Fluid Dyn. Res. 49 (3) (2017) 035512.

7 8

References

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 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

RO

16

DP

15

TE

14

EC

13

CO RR

12

UN

11

[1] Brian P. Danowsky, et al., Incorporation of feedback control into a high-fidelity aeroservoelastic fighter aircraft model, J. Aircr. 47 (4) (2012) 1274–1282. [2] Paolo Domenici, R.W. Blake, The kinematics and performance of the escape response in the angelfish (Pterophyllum eimekei), Can. J. Zool. 71 (11) (1991) 2319–2326. [3] M.F. Platzer, et al., Flapping wing aerodynamics: progress and challenges, AIAA J. 46 (9) (2008) 2136–2149. [4] W. Shyy, et al., Aerodynamics, sensing and control of insect-scale flapping-wing flight, Proc. Math. Phys. Eng. Sci. 472 (2186) (2016) 20150712. [5] Farooq Umar, M. Sun, Aerodynamic force and power for flapping wing at low Reynolds number in ground effect, in: International Bhurban Conference on Applied Sciences and Technology, 2018, pp. 553–558. [6] George V. Lauder, et al., Design and Performance of a Fish Fin-Like Propulsor for AUVs, 2005. [7] Earl H. Dowell, Some recent advances in nonlinear aeroelasticity, in: A Modern Course in Aeroelasticity, 2015, pp. 609–649. [8] Wolfgang Geissler, B.G.V.D. Wall, Dynamic stall control on flapping wing airfoils, Aerosp. Sci. Technol. 62 (2017) 1–10. [9] J. Gray, Studies in animal locomotion V: resistance reflexes in the Eel, J. Exp. Biol. 13 (2) (1936) 181–191. [10] M.J. Lighthill, Large-amplitude elongated-body theory of fish locomotion, Proc. R. Soc. Lond. 179 (1055) (1971) 125–138. [11] Che Shu Lin, C. Hwu, W.B. Young, The thrust and lift of an ornithopter’s membrane wings with simple flapping motion, Aerosp. Sci. Technol. 10 (2) (2006) 111–119. [12] Shuanghou Deng, B.V. Oudheusden, Wake structure visualization of a flappingwing micro-air-vehicle in forward flight, Aerosp. Sci. Technol. 50 (2016) 204–211. [13] Mittal Dong, Najjar, Wake topology and hydrodynamic performance of lowaspect-ratio flapping foils, J. Fluid Mech. 566 (566) (2006) 309–343. [14] Yun Fei Zhang, Y.E. Zheng-Yin, Three-dimensional Navier–Stokes simulation of flapping wing with pitching elasticity, J. Aerosp. Power 26 (4) (2011) 880–889. [15] Chengyu Li, H. Dong, Three-dimensional wake topology and propulsive performance of low-aspect-ratio pitching-rolling plates, Phys. Fluids 28 (7) (2016) 2413. [16] M. Moriche, O. Flores, M. García-Villalba, Three-dimensional instabilities in the wake of a flapping wing at low Reynolds number, Int. J. Heat Fluid Flow (2016). [17] Kosuke Suzuki, M. Yoshino, Aerodynamic comparison of a butterfly-like flapping wing-body model and a revolving-wing model, Fluid Dyn. Res. 49 (3) (2017) 035512. [18] S.M. Dash, et al., Enhanced thrust performance of a two dimensional elliptic airfoil at high flapping frequency in a forward flight, J. Fluids Struct. 76 (2018) 37–59. [19] Dengtao Yu, et al., Numerical study of the effect of motion parameters on propulsive efficiency for an oscillating airfoil, J. Fluids Struct. 68 (2017) 245–263. [20] J.M. Anderson, et al., Oscillating foils of high propulsive efficiency, J. Fluid Mech. 360 (360) (1998) 41–72. [21] Michael S. Triantafyllou, A.H. Techet, F.S. Hover, Review of experimental work in biomimetic foils, IEEE J. Ocean. Eng. 29 (3) (2004) 585–594. [22] R. Godoy-Diana, J.L. Aider, J.E. Wesfreid, Transitions in the wake of a flapping foil, Phys. Rev. E, Stat. Nonlinear Soft Matter Phys. 77 (2) (2008) 016308. [23] Kim Boon Lua, H. Lu, T.T. Lim, Rotating elliptic cylinders in a uniform cross flow, J. Fluids Struct. 78 (2018) 36–51. [24] K.B. Lua, et al., On the thrust performance of a flapping two-dimensional elliptic airfoil in a forward flight, J. Fluids Struct. 66 (2016) 91–109. [25] Sunetra Sarkar, K. Venkatraman, Numerical simulation of thrust generating flow past a pitching airfoil, Comput. Fluids 35 (1) (2006) 16–42. [26] Jie Zhang, N. Liu, L.U. Xiyun, Locomotion of a passively flapping flat plate, J. Fluid Mech. 659 (659) (2010) 43–68. [27] R. Ramamurti, W. Sandberg, A computational investigation of the threedimensional unsteady aerodynamics of drosophila hovering and maneuvering, J. Exp. Biol. 210 (5) (2007) 881–896.

[30] D.A. Read, F.S. Hover, M.S. Triantafyllou, Forces on oscillating foils for propulsion and maneuvering, J. Fluids Struct. 17 (1) (2003) 163–183. [31] F.S. Hover, Ø. Haugsdal, M.S. Triantafyllou, Effect of angle of attack profiles in flapping foil propulsion, J. Fluids Struct. 19 (1) (2004) 37–47. [32] L. Schouveiler, F.S. Hover, M.S. Triantafyllou, Performance of flapping foil propulsion, J. Fluids Struct. 20 (7) (2005) 949–959. [33] M. Rosén, G.R. Spedding, A. Hedenström, Wake structure and wingbeat kinematics of a house-martin Delichon urbica, J. R. Soc. Interface 4 (15) (2007) 659–668. [34] Hikaru Aono, F. Liang, H. Liu, Near- and far-field aerodynamics in insect hovering flight: an integrated computational study, J. Exp. Biol. 211 (2) (2008) 239–257. [35] Chengyu Li, H. Dong, Three-dimensional wake topology and propulsive performance of low-aspect-ratio pitching-rolling plates, Phys. Fluids 28 (7) (2016) 2413. [36] Charles S. Peskin, Flow patterns around heart valves, in: Lecture Notes in Physics, vol. 10, 1973, pp. 214–221. [37] F.B. Tian, et al., Interaction between a flexible filament and a downstream rigid body, Phys. Rev. E, Stat. Nonlinear Soft Matter Phys. 82 (2) (2010) 026301. [38] Kourosh, Leading edge strengthening and the propulsion performance of flexible ray fins, J. Fluid Mech. 693 (2) (2012) 402–432. [39] Cyrus K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (1) (2010) 439–472. [40] Jourabian Mahmoud, et al., Lattice Boltzmann simulation of melting phenomenon with natural convection from an eccentric annulus, Therm. Sci. 17 (3) (2013) 877–890. [41] Antonio F. Miguel, Non-Darcy porous media flow in no-slip and slip regimes, Therm. Sci. 16 (1) (2012) 167–176. [42] Zhiwei Tian, et al., Reactive transport LBM model for CO2 injection in fractured reservoirs, Comput. Geosci. 86 (C) (2016) 15–22. [43] Zhi Gang Feng, E.E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems, J. Comput. Phys. 195 (2) (2004) 602–628. [44] Zhi Gang Feng, E.E. Michaelides, Proteus: a direct forcing method in the simulations of particulate flows, J. Comput. Phys. 202 (1) (2005) 20–51. [45] Jie Wu, et al., Fluid dynamics of flapping insect wing in ground effect, J. Bionics Eng. 11 (1) (2014) 52–60. [46] Jie Zhang, N. Liu, L.U. Xiyun, Locomotion of a passively flapping flat plate, J. Fluid Mech. 659 (659) (2010) 43–68. [47] Takaji Inamuro, Y. Kimura, K. Suzuki, Flight simulations of a two-dimensional flapping wing by the IB-LBM, Int. J. Mod. Phys. C 25 (01) (2014) 4191. [48] S. Succi, The lattice Boltzmann equation – for fluid dynamics and beyond, Phys. Today 55 (12) (2002) 58–60. [49] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479–517. [50] D. Goldstein, R. Handler, L. Sirovich, Modeling a no-slip flow boundary with an external force field, J. Comput. Phys. 105 (2) (1993) 354–366. [51] Wei Xi Huang, S.J. Shin, H.J. Sung, Simulation of flexible filaments in a uniform flow by the immersed boundary method, J. Comput. Phys. 226 (2) (2007) 2206–2228. [52] J. Mohd-Yusof, Combined Immersed Boundary/B-Spline Methods for Simulation of Flow in Complex Geometries, Annual Research Briefs Center for Turbulence Research, 1997. [53] Kosuke Suzuki, K. Minami, T. Inamuro, Lift and thrust generation by a butterflylike flapping wing–body model: immersed boundary-lattice Boltzmann simulations, J. Fluid Mech. 767 (2015) 659–695. [54] C.L. Gong, et al., Numerical investigation of unsteady flows past flapping wings with immersed boundary-lattice Boltzmann method, J. Mech. (2017) 1–15. [55] Zong Jing Yuan, J.I. Xing, G. Chen, Numerical simulation of unsteady flow past biomimetic wing with IB-LBM, Phys. Gases 2 (1) (2017) 39–47. [56] James Buchholz, M. Green, A. Smits, Pressure distribution, thrust performance, and wake structure of a low-aspect ratio pitching panel, in: Meeting of the APS Division of Fluid Dynamics, American Physical Society, 2006. [57] K. Parker, K.D.V. Ellenrieder, J. Soria, Using stereo multigrid DPIV (SMDPIV) measurements to investigate the vortical skeleton behind a finite-span flapping wing, Exp. Fluids 39 (2) (2005) 281–298.

OF

9 10

13

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126

61

127

62

128

63

129

64

130

65

131

66

132