- Email: [email protected]
Ground-induced lift enhancement in a tandem of symmetric flapping wings: Lattice Boltzmann-immersed boundary simulations
Computers and Structures 153 (2015) 230–238
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
Ground-induced lift enhancement in a tandem of symmetric flapping wings: Lattice Boltzmann-immersed boundary simulations Alessandro De Rosis Laboratoire de Mécanique des Fluides et d’Acoustique (LMFA), École Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Écully Cedex, France
a r t i c l e
i n f o
Article history: Received 20 August 2014 Accepted 15 February 2015
Keywords: Fluid–structure interaction Flapping wings Lattice Boltzmann method Immersed boundary
a b s t r a c t The behavior of a tandem of symmetric flapping wings immersed in a quiescent viscous fluid is numerically dissected. The attention focuses on the effect on the flight performance of a solid surface which idealizes the presence of the ground. A wide numerical campaign is carried out. The author demonstrates that the presence of a solid surface can drastically modify the lift force, thus giving a remarkable advantage for the vertical take-off. Therefore, a proper governing parameter is proposed, which accounts for the ratio between the initial gap from the solid surface and the length of the wing. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Flapping wing dynamics represents one of the most intriguing topic affecting contemporary numerical modeling within the computation fluid dynamics (CFD) framework. In fact, such problem shows a lot of practical applications, such as the design of microair vehicles [1–6] or the improvement of the performance of underwater energy harvesters [7–9]. In the last decade, a lot of works developed aiming at dissecting the dynamics of flapping wings immersed in an unbounded quiescent viscous fluid by considering rigid [10–13], flexible [14–20], and even composite [21–23] structures. An interesting aspect worth to be investigated is represented by the force enhancement induced by the presence of a solid surface, as the ground [24]. For example, the ground effect involving wings composed of single and double elements has been investigated in [25,26], respectively, while the aerodynamics of Gurney flaps has been dissected in [27], showing that the presence of the ground can considerably modify the dynamics of a flapping wing. Motion control, structural design and materials involving the technology of the wing undergoing the ground effect has been reviewed in [28], by devoting with special attention for the mathematical modeling. Moreover, in [29] a wind tunnel investigation of a low aspect ratio wing has been carried out, showing the lift augmentation induced by the ground effect. In addition, it has been demonstrated that the fluid forces acting upon a rigid sharp-edged lamina undergoing harmonic oscillations in a quiescent viscous fluid strongly depends E-mail address: [email protected] http://dx.doi.org/10.1016/j.compstruc.2015.02.016 0045-7949/Ó 2015 Elsevier Ltd. All rights reserved.
on the gap-over-length ratio, that is the ratio between the distance from the solid surface and the length of the body [30–32]. Aiming at predicting the flow physics induced by the motion of flapping wings in a fluid, numerical methods can be successfully adopted. Among the possible approaches able to compute the fluid flow, the lattice Boltzmann (LB) method [33–35] represents a relatively new, computationally efficient, and accurate alternative to the standard continuum-based CFD solvers. In particular, the author remarks that the LB method recovers the Navier–Stokes equations for incompressible flow with second-order of accuracy [36]. Moreover, the computation of the forces acting upon an immersed solid body represents a quite easy task, as in the numerical simulations carried out in this paper. In order to account for the presence of the wings which are immersed in the fluid domain, the Immersed Boundary (IB) method [37–40] is employed. Such method has been preferred to the well consolidated interpolated bounce-back scheme [41–43], due to its superior properties in terms of stability and involved computational effort, as recently demonstrated by the author [44]. Once the forces that the fluid exerts upon the wings are computed, wings solid dynamics is predicted through the time discontinuous Galerkin (TDG) scheme, which exhibits superior properties with respect to standard Newmark or a schemes in terms of stability, accuracy and convergence, as discussed in [45]. The LB, IB and TDG methods have been combined by the author within a proper coupling strategy [46], whose effectiveness has already been widely tested against problems involving flapping wings dynamics [13,20], shallow waters [47], multiphase flows [48], hull slamming [49] and even non-Newtonian fluids [50], among the others.
231
A. De Rosis / Computers and Structures 153 (2015) 230–238
In this work, findings from numerical simulations concerning a tandem of rigid symmetric flapping wings undergoing harmonic oscillations are discussed. Specifically, wings can move only in vertical direction and the flight performance is elucidated by highlighting the role of the gap-over-length ratio. Moreover, scenarios characterized by different fluid viscosity are investigated, together with the influence of the reciprocal distance between the wings. Results in terms of time history of the position of the centers of mass are discussed, together with considerations about the velocity field. This paper presents new insights with respect to previous efforts carried out by different authors [13,20,51–54]. Specifically, the behavior of a tandem of butterfly-like wings is investigated, whereas the literature focuses mainly on isolated bodies. In this way, the mutual interaction between the wings is highlighted, together with the role of the encompassing hydrodynamics. Moreover, the effect of a solid ground-like surface is widely elucidated, showing its huge impact on the flight performance. The paper is organized as follows. In Section 2, the problem is stated and the adopted numerical methods are briefly recalled. In Section 3, results of a numerical campaign are discussed. Some conclusions are drawn in Section 4. Finally, in Appendix A the convergence and accuracy properties of the proposed algorithm are discussed.
y
20 L
A
B
g
z
2. Problem statement A fluid of viscosity m and density q surrounds a tandem of two symmetric flapping wings. The following parameters are adopted, which correspond to a Pieris melete butterfly [55,56]: wing mass 3:5 106 kg, body mass 5:0 105 kg, hinge-wing distance 5:0 103 m, wing length L ¼ 3:0 102 m. Making reference to Fig. 1, each wing undergoes a harmonic motion, that is
hðtÞ ¼ Dh cos
2pt ; T
ð1Þ
where Dh ¼ 46:8 is the amplitude, T ¼ 0:1 s is the period of the harmonic oscillation and t is the time. Notice that the total mass M ¼ 5:7 105 kg is considered to be concentrated at the hinge. This means that the update of the position of the wings reduces to the one of the corresponding center of mass. At the bottom-most section, a rigid solid surface is enforced, whereas outflow boundary conditions are prescribed at the remaining edges. The fluid and the wings are initially at rest. Wings are initially horizontally aligned. Moreover, the initial distance between the wings and the solid surface is denoted by g. Thus, it is possible to define a governing parameter, the gap-over-length ratio, as G ¼ g=L. The dependence of the solution of the problem on G is discussed in the following, together with the influence of the Reynolds number Re and the incidence of the reciprocal distance between the wings. The problem set-up is sketched in Fig. 2.
L
L
M Fig. 1. Symmetric flapping wings of length L and mass M undergoing a harmonic motion of maximum amplitude Dh. The total mass M is concentrated at the hinge (red circle). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
10 L Fig. 2. Sketch of the problem set-up. Two symmetric flapping wings (A and B, red solid lines) of length L and mass M are immersed in a viscous fluid. These are placed symmetrically with respect to the vertical axis y. A solid surface is prescribed at the bottom-most section (blue solid line). Outflow conditions are enforced at the other boundaries. The centers of mass are initially horizontally aligned and these are placed at a distance g from the solid surface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
2.1. Governing equations The flow physics is predicted by solving the Navier–Stokes equations for an incompressible flow, which can be written as
r u ¼ 0;
ð2Þ
@u 1 þ ðu rÞðuÞ ¼ rp þ mr2 u; @t q
ð3Þ
where u is the flow velocity and p is the fluid pressure. From the solid dynamics point of view, the wings (i.e the center of mass) is forced to move only in the vertical direction according to equation of the solid motion, that is
€ðtÞ ¼ FðtÞ; My
ð4Þ
where y is the vertical component of the displacement of the center of mass of the wings, FðtÞ is the time-dependent vertical component of the external forces which are exerted by the encompassing fluid on the solid body. Superimposed dots indicate the time derivatives. At the fluid-wing interface and at the bottom section, i.e. at the solid surface, the no-slip condition is prescribed, i.e. uðy ¼ 0Þ ¼ 0. The following initial conditions complete the definition of the problem: _ ¼ 0Þ ¼ 0 and uðt ¼ 0Þ ¼ 0. yðt ¼ 0Þ ¼ 0; yðt 2.2. Numerical methods The lattice D2Q9 model [35] is adopted to compute the space– time evolution of the particle distribution functions f i , which are
A. De Rosis / Computers and Structures 153 (2015) 230–238
20
20
15
15
y/L
y/L
232
10 5 0
10 5
0
5
10
0
15
0
5
t/T
20
20
15
15
10 5
10 5
0
5
10
0
15
0
5
t/T
20
15
15
y/L
y/L
(d) B. C1.
20
10 5
10 5
0
5
10
0
15
0
5
t/T
(f) B. C2.
20
20
15
15
y/L
y/L
15
10
t/T
(e) A. C2.
10 5 0
15
10
t/T
(c) A. C1.
0
15
(b) B. C1.
y/L
y/L
(a) A. C1.
0
10
t/T
10 5
0
5
10
15
0
5
0
t/T
10
15
t/T
(g) A. C2.
(h) B. C2.
Fig. 3. Re ¼ 50. Time history of the position of the centers of mass of A and B in the configurations C1 and C2 for different values of G, i.e.1.25 (red solid line), 2.5 (green solid line), 3.75 (blue solid line), 5 (magenta solid line), 6.25 (red dotted line), 7.5 (green dotted line), 8.75 (blue dotted line), 10 (magenta dotted line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
forced to move along nine prescribed directions i with velocities ci . This equation reads as follows:
f i ðx þ Dt ci ; t þ DtÞ ¼ f i ðx; tÞ þ
1
s
Mach number [34], the equilibrium particle distribution functions eq f i can be computed as
"
eq ½f i ðx; tÞ
f i ðx; tÞ þ 3wi c i gðx; tÞ; ð5Þ
where x is the position, Dt is the time step and s is the so-called relaxation parameter. It is worth to notice that the relaxation parameter s is strictly related to the fluid viscosity as P m ¼ ðs 1=2Þc2s , being c2s ¼ wi c2i ¼ 1=3 with wi a set of 9 weights defined as w0 ¼ 4=9; w1 ¼ w2 ¼ w3 ¼ w4 ¼ 1=9 and w5 ¼ w6 ¼ w7 ¼ w8 ¼ 1=36. By performing a second-order expansion in the local
eq f i ðx; tÞ
¼ wi q 1 þ
v ci c2s
# ðv c i Þ2 v v þ : 2c4s 2c2s
ð6Þ
The macroscopic fluid density q and the flow velocities v are defined as
qðx; tÞ ¼
X f i ðx; tÞ; i
v ðx; tÞ ¼
P
i f i ðx; tÞci
qðx; tÞ
;
ð7Þ
respectively. In order to account for the presence of the wings in the fluid lattice background, the immersed boundary method [38] is
233
20
20
15
15
y/L
y/L
A. De Rosis / Computers and Structures 153 (2015) 230–238
10 5
5
0
5
10
0
15
(b) B. C1.
15
15
10 5
0
5
10
0
15
0
5
10
t/T
t/T
(c) A. C1.
(d) B. C1. 20
15
15
y/L
y/L
5
10 5
15
10 5
0
5
10
0
15
0
5
10
t/T
t/T
(e) A. C2.
(f) B. C2.
20
20
15
15
y/L
y/L
15
10
20
10 5 0
10
(a) A. C1. 20
0
5
t/T
20
0
0
t/T
y/L
y/L
0
10
15
10 5
0
5
10
15
0
0
5
10
t/T
t/T
(g) A. C2.
(h) B. C2.
15
Fig. 4. Re ¼ 100. Time history of the position of the centers of mass of A and B in the configurations C1 and C2 for different values of G, i.e.1.25 (red solid line), 2.5 (green solid line), 3.75 (blue solid line), 5 (magenta solid line), 6.25 (red dotted line), 7.5 (green dotted line), 8.75 (blue dotted line), 10 (magenta dotted line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
adopted. Specifically, each wing is idealized by a set of Lagrangian solid points, namely X j , whose reciprocal spacing is set to a half of the fluid mesh size. The IB is adopted in an implicit velocitycorrection strategy, leading to a correction term gðx; tÞ which is used in the right-hand side of Eq. (5) [44,46]. It is worth to notice that fluid forces, which are used in Eq. (4), can be easily and immediately computed as
Fj ¼
X gðxÞWðx X j Þ;
dynamics. The adopted numerical methods are combined within a staggered explicit coupling algorithm. Its high properties in terms of stability and accuracy have already been widely tested by the author by computing the error in terms of the energy that is artificially introduced in the system at the fluid–structure interface [46,47]. Finally, notice that the problem is solved in LB units. This means that particular attention should be paid to the conversion from the physical system to the LB world.
ð8Þ
x
3. Results and discussion W being the IB kernel [38]. Once fluid forces acting upon the wings are computed, the solid position is updated by adopting the TDG method. The interested reader can refer to [45,57,58] for further details about such time integration scheme for predicting solid
In this section, results from a numerical campaign are discussed. Each wing is idealized by Llb ¼ 100 lattice nodes and it is represented by 200 Lagrangian immersed boundary points.
A. De Rosis / Computers and Structures 153 (2015) 230–238
20
20
15
15
y/L
y/L
234
10 5 0
10 5
5
0
10
0
15
5
0
10
t/T
(b) B. C1.
20
20
15
15
y/L
y/L
(a) A. C1.
10 5 0
10 5
0
5
10
0
15
0
5
10
t/T
20
15
15
y/L
y/L
(d) B. C1.
20
10 5
10 5
0
5
10
0
15
0
5
10
t/T
(f) B. C2.
20
20
15
15
y/L
y/L
15
t/T
(e) A. C2.
10 5 0
15
t/T
(c) A. C1.
0
15
t/T
10 5
0
5
10
0
15
0
5
10
t/T
15
t/T
(g) A. C2.
(h) B. C2.
10
8
8
6
6
yt /L
yt /L
Fig. 5. Re ¼ 200. Time history of the position of the center of mass of A and B in the configurations C1 and C2 for different values of G, i.e.1.25 (red solid line), 2.5 (green solid line), 3.75 (blue solid line), 5 (magenta solid line), 6.25 (red dotted line), 7.5 (green dotted line), 8.75 (blue dotted line), 10 (magenta dotted line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4
2
2 0
4
2.5
7.5
5
10
0
2.5
5
7.5
g/L
g/L
(a) A.
(b) B.
10
Fig. 6. Position achieved at t=T ¼ 6, namely yt =L, by A and B vs the gap-over-length ratio, g=L, for different values of the Reynolds number, i.e. Re ¼ 50 (red lines), 100 (green) and 200 (blue). Configurations C1 and C2 are represented by solid and dotted lines, respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
235
A. De Rosis / Computers and Structures 153 (2015) 230–238 Table 1 A. Position achieved at t=T ¼ 6 for different values of G and Re in the configurations C1 and C2. The peak value is highlighted by bold. G
1.25 2.5 3.75 5 6.25 7.5 8.75 10
Re ¼ 50
Re ¼ 100
Re ¼ 200
C1
C2
C1
C2
C1
C2
8.6353 5.1355 3.4593 2.5957 2.2740 2.2103 2.2401 2.3745
3.0715 1.9865 1.6076 1.3911 1.2902 1.2785 1.2955 1.3176
8.9457 3.0395 1.9508 1.3675 1.0388 0.9179 0.9777 1.1439
1.9823 1.3662 1.2304 1.0391 0.8819 0.8615 0.9162 1.0105
9.8773 4.8161 2.9350 1.8429 1.3233 1.0071 1.0805 1.1974
1.0885 0.6463 0.3788 0.3296 0.2858 0.1156 0.1639 0.2214
Table 2 B. Position achieved at t=T ¼ 6 for different values of G and Re in the configurations C1 and C2. The peak value is highlighted by bold. G
1.25 2.5 3.75 5 6.25 7.5 8.75 10
Re ¼ 50
Re ¼ 100
Re ¼ 200
C1
C2
C1
C2
C1
C2
7.7949 4.4192 3.0824 2.3642 2.0768 2.0117 2.0430 2.1779
3.0334 1.9457 1.5816 1.3712 1.2682 1.2528 1.2725 1.3014
5.7789 2.6741 1.8413 1.3442 1.0850 1.0060 1.0359 1.1925
1.9365 1.3516 1.2287 1.0473 0.9001 0.8820 0.9310 1.0173
5.3232 2.9726 1.8029 1.4218 1.1144 0.8423 0.9508 1.1802
0.7726 0.3873 0.2639 0.2581 0.2029 0.1117 0.1116 0.1628
(a) t/T = 0.76.
(b) t/T = 2.29.
(c) t/T = 6.89.
(d) t/T = 8.42.
Moreover, the grid consists of 1000 and 2000 lattice nodes in the horizontal and vertical directions, respectively. All the LB simulations are carried out at low Mach number Ma, i.e. Ma ¼ v tip;lb =cs ¼ 0:0173, where v tip;lb ¼ 4Llb Dh=T lb ¼ 0:01. Moreover, the initial density is set to q ¼ 1. In order to help the reader in the unit conversion, the author defines the lattice grid spacing as
Dx ¼ L=Llb ¼ 3 104 m and the time step Dt ¼ 3:061 106 s. Therefore, the period of oscillation in the LB simulations is set to T lb ¼ T=Dt ¼ 32672:56 lattice time steps. Simulations stop when A or B overcomes the fluid domain limits or if t=T > 15. Scenarios characterized by different values of the relaxation parameter s are investigated, thus corresponding to distinct values of the Reynolds number, i.e. Re ¼ v tip;lb Llb =m ¼ 50; 100; 200. Moreover, two configurations are dissected. In the former (C1), each center of mass is far 100 lattice nodes from the vertical axis of symmetry, whereas in the latter (C2) such distance is set to 200 lattice nodes. At the bottom-most section the no-slip condition is enforced, i.e. v ¼ 0. At the remaining boundaries, outflow conditions are prescribed, i.e. n rf i ¼ 0; n being the outer normal to the boundary. In Fig. 3, the time history of the position of the centers of mass of A and B are depicted in the configurations C1 and C2 for different values of G, i.e. 1.25, 2.5, 3.75, 5, 6.25, 7.5, 8.75, 10. As it is possible to observe, two different patterns arise. The former collects the scenarios characterized by G 6 5 (solid lines), where the lower G, the faster the take-off is. At a first glance, such behavior can be confusing. In fact, A and B achieve higher values of the position if these start from a lower initial abscissa. On the other hand, if G P 5 (dotted lines), such trend drastically changes, demonstrating that the flight performance are remarkably dominated by the initial position. Indeed, in the latter pattern, the higher the initial position, the higher the final position is. This dramatic variation has to be addressed to the role of the solid surface which is located at the bottom-most section. In fact, as the wings are close to the boundary, these exhibit a progressively enhancement in terms of generated lift force, thus corresponding to higher values of the achieved position. As demonstrated in [32], if a solid body is immersed in a
15 (e) Color legend. Fig. 7. Most advantageous scenarios: velocity field at different time instants for Re ¼ 200 and G ¼ 1:25 in C1.
viscous fluid and it oscillates close to a solid surface, a growth of the fluid force arises, consisting of one order of magnitude with respect to the case in which the solid surface is neglected. Therefore, it is possible to assess that the flapping wings considerably benefit from the presence of this surface. Such behavior is known as the so-called ground effect. Moreover, the behavior is even dominated by the reciprocal distance between A and B. In fact, the wings appear to benefit from being relatively close. This
236
A. De Rosis / Computers and Structures 153 (2015) 230–238
(a) t/T = 0.76.
(b) t/T = 2.29.
reported in Fig. 5. Again, the deleterious effects experienced at G > 5 are still present. On the other hand, configurations C1 and C2 show different trends. In the former, curves appears to be shifted back in time, whereas the latter presents a time shift in advance with respect to the scenarios characterized by Re ¼ 100. This means that the solution of the problem is highly non-linear, since it is strongly dependent on G; Re and the horizontal reciprocal distance. Different patterns arose due to the distinct tip-induced vorticity which involves the flow field. As well known, it is related to the fluid viscosity and to the distance between solid surfaces. In order to summarize the above mentioned effects, in Fig. 6 the position achieved at t=T ¼ 6, namely yt =L, by A and B is plotted against the gap-over-length ratio, G ¼ g=L, for different values of the Reynolds number in the configurations C1 and C2. As it is possible to observe, scenarios corresponding to the configuration C1 are considerably higher with respect to the ones related to C2. As G decreases, the benefit is higher, whereas curves plateau for progressively larger values of the gap-over-length ratio. In addition, notice that a more marked advantage is experienced by A, since B shows lower values of yt =L. In order to help the reader, data depicted in Fig. 6 are reported in Tables 1 and 2 for A and B, respectively. Finally, making reference to the flight performance of A, in Figs. 7 and 8 the velocity field is depicted at different time instants in the most advantageous configuration, i.e Re ¼ 200 and G ¼ 1:25 in C1, and the most arduous one, Re ¼ 200 and G ¼ 7:5 in C2. The above mentioned considerations about the role of the reciprocal distance are confirmed. In particular, notice that if the wings are closer, a more marked vorticity arises in the flow field. Vortexes are mainly located below the wing, thus considerably contributing to the lift force enhancement. On the other hand, the scenario manifesting an arduous take-off is characterized by low values of the velocity. 4. Conclusions
(c) t/T = 6.89.
(d) t/T = 8.42.
16 (e) Color legend. Fig. 8. Most arduous scenarios: velocity field at different time instants for Re ¼ 200 and G ¼ 7:5 in C2.
trend has to be addressed to the tip-induced vorticity, which highly supports wings traveling. Configurations corresponding to Re ¼ 100 are shown in Fig. 4. Notice that curves appear to be shifted in advance in time with respect to the previous case. The threshold represented by G 5 is already experienced, since two different trend in the flight performance are found. In addition, the configuration C2 exhibits a remarkably arduous take-off, except for the peak value of G. Finally, scenarios at Re ¼ 200 are
In this paper, the dynamics involving a tandem of symmetric flapping wings immersed in a viscous fluid near a solid surface has been investigated. Specifically, the dependence of the flight performance in terms of time history of the position of the centers of mass on the fluid viscosity, i.e. the Reynolds number Re, has been elucidated. Moreover, the so-called ground effect has been dissected by performing simulations characterized by progressively higher values of the gap-over-length ratio, G. The author demonstrated that below a certain threshold, i.e. G 6 5, wings dynamics is positively affected by the presence of the solid surface. In fact, a remarkable enhancement of the flight performance has been found. If G > 5, the advantageous ground effect vanishes. Finally, the incidence of wings reciprocal distance has been elucidated, thus underlining the mutual interaction between A and B, especially in terms of tip-induced vorticity. It has been demonstrated that wings benefit from a closer arrangement, thus justifying the natural attitude to move in flocks. Finally, it is possible to assess that the above adopted numerical algorithms can be successfully employed to solve problems involving the dynamics of flapping wings which are immersed in a viscous fluid. Acknowledgement The author would like to thank Prof. Stefano Ubertini and Prof. Francesco Ubertini for the useful and motivating discussions. Appendix A The convergence and accuracy properties of the adopted numerical methods are evaluated in the following convergence
A. De Rosis / Computers and Structures 153 (2015) 230–238
Error
-2
-3
-4
5
5.5
6
6.5
Lattice nodes Fig. 9. Convergence analysis: relative error vs number of lattice nodes in a log–log scale.
analysis. Specifically, the value of yt =L at Re ¼ 200 and G ¼ 1:25 in C1 for A is computed by progressively refining the lattice fluid grid. For a given grid resolution, such value is denoted ecurr . It is compared to the same quantity computed on a very fine grid composed of 4000 2000 lattice nodes and it is denoted eref . Then, it is possible to define the relative error as
Error ¼
eref ecurr : eref
ð9Þ
In Fig. 9, the relative error is depicted against the number of lattice nodes. Both the quantities are given in a logarithmic scale. As it is possible to observe, a good convergence rate equal to 1.78 is experienced, thus allowing us to assess that the proposed strategy shows high properties in terms of accuracy and convergence. References [1] Ellington CP. The novel aerodynamics of insect flight: applications to micro-air vehicles. J Exp Biol 1999;202:3439–48. [2] Shyy W, Berg M, Ljungqvist D. Flapping and flexible wings for biological and micro air vehicles. Prog Aerosp Sci 1999;35:455–505. [3] Ettinger SM, Nechyba MC, Ifju PG, Waszak M. Vision-guided flight stability and control for micro air vehicles. Adv Robot 2003;17:617–40. [4] Ansari S, Abikowski R, Knowles K. Aerodynamic modelling of insect-like flapping flight for micro air vehicles. Prog Aerosp Sci 2006;42:129–72. [5] Shyy W, Aono H, Chimakurthi SK, Trizila P, Kang C-K, Cesnik CE, et al. Recent progress in flapping wing aerodynamics and aeroelasticity. Prog Aerosp Sci 2010;46:284–327. [6] Orlowski CT, Girard AR. Dynamics, stability, and control analyses of flapping wing micro-air vehicles. Prog Aerosp Sci 2012;51:18–30. [7] Aureli M, Prince C, Porfiri M, Peterson SD. Energy harvesting from base excitation of ionic polymer metal composites in fluid environments. Smart Mater Struct 2010;19:015003. [8] Giacomello A, Porfiri M. Underwater energy harvesting from a heavy flag hosting ionic polymer metal composites. J Appl Phys 2011;109:084903. [9] Erturk A, Delporte G. Underwater thrust and power generation using flexible piezoelectric composites: an experimental investigation toward self-powered swimmer-sensor platforms. Smart Mater Struct 2011;20:125013. [10] Rozhdestvensky KV, Ryzhov VA. Aerohydrodynamics of flapping-wing propulsors. Prog Aerosp Sci 2003;39:585–633. [11] Trizila P, Kang C-K, Aono H, Shyy W, Visbal M. Low-reynolds-number aerodynamics of a flapping rigid flat plate. AIAA J 2011;49:806–23. [12] Fuchiwaki M, Kuroki T, Tanaka K, Tababa T. Dynamic behavior of the vortex ring formed on a butterfly wing. Exp Fluids 2013;54:1–12. [13] De Rosis A. On the dynamics of a tandem of asynchronous flapping wings: lattice Boltzmann-immersed boundary simulations. Physica A 2014;410:276–86. [14] Heathcote S, Gursul I. Flexible flapping airfoil propulsion at low reynolds numbers. AIAA J 2007;45:1066–79. [15] Heathcote S, Wang Z, Gursul I. Effect of spanwise flexibility on flapping wing propulsion. J Fluids Struct 2008;24:183–99. [16] Majumder L, Rao S. Interval-based optimization of aircraft wings under landing loads. Comput Struct 2009;87:225–35. [17] Masoud H, Alexeev A. Resonance of flexible flapping wings at low reynolds number. Phys Rev E 2010;81:056304. [18] Mountcastle AM, Combes SA. Wing flexibility enhances load-lifting capacity in bumblebees. Proc R Soc B 2013;280. [19] Kang C-k, Shyy W. Scaling law and enhancement of lift generation of an insectsize hovering flexible wing. J R Soc Interface 2013;10.
237
[20] De Rosis A, Falcucci G, Ubertini S, Ubertini F. Aeroelastic study of flexible flapping wings by a coupled lattice Boltzmann finite element approach with immersed boundary method. J Fluids Struct 2014;49:516–33. [21] Kim T-U, Hwang IH. Optimal design of composite wing subjected to gust loads. Comput Struct 2005;83:1546–54. [22] Kameyama M, Fukunaga H. Optimum design of composite plate wings for aeroelastic characteristics using lamination parameters. Comput Struct 2007;85:213–24. [23] Kim T-U, Shin JW, Hwang IH. Stacking sequence design of a composite wing under a random gust using a genetic algorithm. Comput Struct 2007;85:579–85. Stochastic Structural Analysis, Optimization and Re-Analysi. [24] Jones K, Platzer M. Numerical computation of flapping-wing propulsion and power extraction. AIAA Paper 1997;97:0826. [25] Zerihan J, Zhang X. Aerodynamics of a single element wing in ground effect. J Aircraft 2000;37:1058–64. [26] Zhang X, Zerihan J. Aerodynamics of a double-element wing in ground effect. AIAA J 2003;41:1007–16. [27] Zerihan J, Zhang X. Aerodynamics of gurney flaps on a wing in ground effect. AIAA J 2001;39:772–80. [28] Rozhdestvensky KV. Wing-in-ground effect vehicles. Prog Aerosp Sci 2006;42:211–83. [29] Ahmed N, Goonaratne J. Lift augmentation of a low-aspect-ratio thick wing in ground effect. J Aircraft 2002;39:381–4. [30] Green CP, Sader JE. Frequency response of cantilever beams immersed in viscous fluids near a solid surface with applications to the atomic force microscope. J Appl Phys 2005;98:114913. [31] Grimaldi E, Porfiri M, Soria L. Finite amplitude vibrations of a sharp-edged beam immersed in a viscous fluid near a solid surface. J Appl Phys 2012;112:104907. [32] De Rosis A. Harmonic oscillations of a lamina in a viscous fluid near a solid surface: a lattice Boltzmann-immersed boundary approach. Physica A 2014;415:386–97. [33] Higuera FJ, Succi S, Benzi R. Lattice gas dynamics with enhanced collisions. Europhys Lett 1989;9:345–9. [34] Benzi R, Succi S, Vergassola M. The lattice Boltzmann equation: theory and applications. Phys Rep 1992;222:145–97. [35] Succi S. The lattice Boltzmann equation for fluid dynamics and beyond. Clarendon; 2001. [36] Chen H, Chen S, Matthaeus WH. Recovery of the Navier–Stokes equations using a lattice-gas Boltzmann method. Phys Rev Lett 1992;45:R5339–42. [37] Fadlun E, Verzicco R, Orlandi P, Mohd-Yusof J. Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. J Comput Phys 2000;161:35–60. [38] Peskin C. The immersed boundary method. Acta Numer 2002;11:479–517. [39] Wu J, Shu C. Implicit velocity correction-based immersed boundarylattice Boltzmann method and its applications. J Comput Phys 2009;228: 1963–79. [40] Sotiropoulos F, Yang X. Immersed boundary methods for simulating fluid– structure interaction. Prog Aerosp Sci 2014;65:1–21. [41] Filippova O, Hänel D. Lattice Boltzmann simulation of gas-particle flow in filters. Comput Fluids 1997;26:697–712. [42] Mei R, Luo L, Shyy W. An accurate curved boundary treatment in the lattice Boltzmann method. J Comput Phys 1999;155:307–30. [43] Mei R, Yu D, Shyy W, Luo L. Force evaluation in the lattice Boltzmann method involving curved geometry. Phys Rev Lett E 2002;65:041203. [44] De Rosis A, Ubertini S, Ubertini F. A comparison between the interpolated bounce-back scheme and the immersed boundary method to treat solid boundary conditions for laminar flows in the lattice Boltzmann framework. J Sci Comput 2014;61:477–89. [45] Mancuso M, Ubertini F. An efficient integration procedure for linear dynamics based on a time discontinuous Galerkin formulation. Comput Mech 2003;32:154–68. [46] De Rosis A, Ubertini S, Ubertini F. A partitioned approach for two-dimensional fluid–structure interaction problems by a coupled lattice Boltzmannfinite element method with immersed boundary. J Fluids Struct 2014;45: 202–15. [47] De Rosis A. A lattice Boltzmann-finite element model for two-dimensional fluid–structure interaction problems involving shallow waters. Adv Water Resour 2014;65:18–24. [48] De Rosis A. A lattice Boltzmann model for multiphase flows interacting with deformable bodies. Adv Water Resour 2014;73:55–64. [49] De Rosis A, Falcucci G, Porfiri M, Ubertini F, Ubertini S. Hydroelastic analysis of hull slamming coupling lattice Boltzmann and finite element methods. Comput Struct 2014;138:24–35. [50] De Rosis A. Harmonic oscillations of laminae in non-Newtonian fluids: a lattice Boltzmann-immersed boundary approach. Adv Water Resour 2014;73:97–107. [51] Srinivas Rao K, Subrahmanya M, Kulkarni D, Rajani B. Numerical simulation of turbulent flow past micro air vehicle wing sections. Bangalore: National Aerospace Laboratories; 2009. [52] Ota K, Suzuki K, Inamuro T. Lift generation by a two-dimensional symmetric flapping wing: immersed boundary-lattice Boltzmann simulations. Fluid Dyn Res 2012;44:045504. [53] Inamuro T, Kimura Y, Suzuki K. Flight simulations of a two-dimensional flapping wing by the ib-lbm. Bull Am Phys Soc 2012;57.
238
A. De Rosis / Computers and Structures 153 (2015) 230–238
[54] Kimura Y, Suzuki K, Inamuro T. Flight simulations of a two-dimensional flapping wing by the ib-lbm. Int J Mod Phys C 2014;25. [55] Sunada S, Kawachi K, Watanabe I, Azuma A. Performance of a butterfly in takeoff flight. J Exp Biol 1993;183:249–77. [56] Iima T, Yanagita M. Asymmetric motion of a two-dimensional symmetric flapping model. Fluid Dyn Res 2005;36:407–25.
[57] Mancuso M, Ubertini F. An efficient time discontinuous Galerkin procedure for non-linear structural dynamics. Comput Methods Appl Mech Eng 2006;195: 6391–406. [58] De Rosis A, Falcucci G, Ubertini S, Ubertini F, Succi S. Lattice Boltzmann analysis of fluid–structure interaction with moving boundaries. Commun Comput Phys 2012;13:823–34.
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
Ground-induced lift enhancement in a tandem of symmetric flapping wings: Lattice Boltzmann-immersed boundary simulations Alessandro De Rosis Laboratoire de Mécanique des Fluides et d’Acoustique (LMFA), École Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Écully Cedex, France
a r t i c l e
i n f o
Article history: Received 20 August 2014 Accepted 15 February 2015
Keywords: Fluid–structure interaction Flapping wings Lattice Boltzmann method Immersed boundary
a b s t r a c t The behavior of a tandem of symmetric flapping wings immersed in a quiescent viscous fluid is numerically dissected. The attention focuses on the effect on the flight performance of a solid surface which idealizes the presence of the ground. A wide numerical campaign is carried out. The author demonstrates that the presence of a solid surface can drastically modify the lift force, thus giving a remarkable advantage for the vertical take-off. Therefore, a proper governing parameter is proposed, which accounts for the ratio between the initial gap from the solid surface and the length of the wing. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Flapping wing dynamics represents one of the most intriguing topic affecting contemporary numerical modeling within the computation fluid dynamics (CFD) framework. In fact, such problem shows a lot of practical applications, such as the design of microair vehicles [1–6] or the improvement of the performance of underwater energy harvesters [7–9]. In the last decade, a lot of works developed aiming at dissecting the dynamics of flapping wings immersed in an unbounded quiescent viscous fluid by considering rigid [10–13], flexible [14–20], and even composite [21–23] structures. An interesting aspect worth to be investigated is represented by the force enhancement induced by the presence of a solid surface, as the ground [24]. For example, the ground effect involving wings composed of single and double elements has been investigated in [25,26], respectively, while the aerodynamics of Gurney flaps has been dissected in [27], showing that the presence of the ground can considerably modify the dynamics of a flapping wing. Motion control, structural design and materials involving the technology of the wing undergoing the ground effect has been reviewed in [28], by devoting with special attention for the mathematical modeling. Moreover, in [29] a wind tunnel investigation of a low aspect ratio wing has been carried out, showing the lift augmentation induced by the ground effect. In addition, it has been demonstrated that the fluid forces acting upon a rigid sharp-edged lamina undergoing harmonic oscillations in a quiescent viscous fluid strongly depends E-mail address: [email protected] http://dx.doi.org/10.1016/j.compstruc.2015.02.016 0045-7949/Ó 2015 Elsevier Ltd. All rights reserved.
on the gap-over-length ratio, that is the ratio between the distance from the solid surface and the length of the body [30–32]. Aiming at predicting the flow physics induced by the motion of flapping wings in a fluid, numerical methods can be successfully adopted. Among the possible approaches able to compute the fluid flow, the lattice Boltzmann (LB) method [33–35] represents a relatively new, computationally efficient, and accurate alternative to the standard continuum-based CFD solvers. In particular, the author remarks that the LB method recovers the Navier–Stokes equations for incompressible flow with second-order of accuracy [36]. Moreover, the computation of the forces acting upon an immersed solid body represents a quite easy task, as in the numerical simulations carried out in this paper. In order to account for the presence of the wings which are immersed in the fluid domain, the Immersed Boundary (IB) method [37–40] is employed. Such method has been preferred to the well consolidated interpolated bounce-back scheme [41–43], due to its superior properties in terms of stability and involved computational effort, as recently demonstrated by the author [44]. Once the forces that the fluid exerts upon the wings are computed, wings solid dynamics is predicted through the time discontinuous Galerkin (TDG) scheme, which exhibits superior properties with respect to standard Newmark or a schemes in terms of stability, accuracy and convergence, as discussed in [45]. The LB, IB and TDG methods have been combined by the author within a proper coupling strategy [46], whose effectiveness has already been widely tested against problems involving flapping wings dynamics [13,20], shallow waters [47], multiphase flows [48], hull slamming [49] and even non-Newtonian fluids [50], among the others.
231
A. De Rosis / Computers and Structures 153 (2015) 230–238
In this work, findings from numerical simulations concerning a tandem of rigid symmetric flapping wings undergoing harmonic oscillations are discussed. Specifically, wings can move only in vertical direction and the flight performance is elucidated by highlighting the role of the gap-over-length ratio. Moreover, scenarios characterized by different fluid viscosity are investigated, together with the influence of the reciprocal distance between the wings. Results in terms of time history of the position of the centers of mass are discussed, together with considerations about the velocity field. This paper presents new insights with respect to previous efforts carried out by different authors [13,20,51–54]. Specifically, the behavior of a tandem of butterfly-like wings is investigated, whereas the literature focuses mainly on isolated bodies. In this way, the mutual interaction between the wings is highlighted, together with the role of the encompassing hydrodynamics. Moreover, the effect of a solid ground-like surface is widely elucidated, showing its huge impact on the flight performance. The paper is organized as follows. In Section 2, the problem is stated and the adopted numerical methods are briefly recalled. In Section 3, results of a numerical campaign are discussed. Some conclusions are drawn in Section 4. Finally, in Appendix A the convergence and accuracy properties of the proposed algorithm are discussed.
y
20 L
A
B
g
z
2. Problem statement A fluid of viscosity m and density q surrounds a tandem of two symmetric flapping wings. The following parameters are adopted, which correspond to a Pieris melete butterfly [55,56]: wing mass 3:5 106 kg, body mass 5:0 105 kg, hinge-wing distance 5:0 103 m, wing length L ¼ 3:0 102 m. Making reference to Fig. 1, each wing undergoes a harmonic motion, that is
hðtÞ ¼ Dh cos
2pt ; T
ð1Þ
where Dh ¼ 46:8 is the amplitude, T ¼ 0:1 s is the period of the harmonic oscillation and t is the time. Notice that the total mass M ¼ 5:7 105 kg is considered to be concentrated at the hinge. This means that the update of the position of the wings reduces to the one of the corresponding center of mass. At the bottom-most section, a rigid solid surface is enforced, whereas outflow boundary conditions are prescribed at the remaining edges. The fluid and the wings are initially at rest. Wings are initially horizontally aligned. Moreover, the initial distance between the wings and the solid surface is denoted by g. Thus, it is possible to define a governing parameter, the gap-over-length ratio, as G ¼ g=L. The dependence of the solution of the problem on G is discussed in the following, together with the influence of the Reynolds number Re and the incidence of the reciprocal distance between the wings. The problem set-up is sketched in Fig. 2.
L
L
M Fig. 1. Symmetric flapping wings of length L and mass M undergoing a harmonic motion of maximum amplitude Dh. The total mass M is concentrated at the hinge (red circle). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
10 L Fig. 2. Sketch of the problem set-up. Two symmetric flapping wings (A and B, red solid lines) of length L and mass M are immersed in a viscous fluid. These are placed symmetrically with respect to the vertical axis y. A solid surface is prescribed at the bottom-most section (blue solid line). Outflow conditions are enforced at the other boundaries. The centers of mass are initially horizontally aligned and these are placed at a distance g from the solid surface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
2.1. Governing equations The flow physics is predicted by solving the Navier–Stokes equations for an incompressible flow, which can be written as
r u ¼ 0;
ð2Þ
@u 1 þ ðu rÞðuÞ ¼ rp þ mr2 u; @t q
ð3Þ
where u is the flow velocity and p is the fluid pressure. From the solid dynamics point of view, the wings (i.e the center of mass) is forced to move only in the vertical direction according to equation of the solid motion, that is
€ðtÞ ¼ FðtÞ; My
ð4Þ
where y is the vertical component of the displacement of the center of mass of the wings, FðtÞ is the time-dependent vertical component of the external forces which are exerted by the encompassing fluid on the solid body. Superimposed dots indicate the time derivatives. At the fluid-wing interface and at the bottom section, i.e. at the solid surface, the no-slip condition is prescribed, i.e. uðy ¼ 0Þ ¼ 0. The following initial conditions complete the definition of the problem: _ ¼ 0Þ ¼ 0 and uðt ¼ 0Þ ¼ 0. yðt ¼ 0Þ ¼ 0; yðt 2.2. Numerical methods The lattice D2Q9 model [35] is adopted to compute the space– time evolution of the particle distribution functions f i , which are
A. De Rosis / Computers and Structures 153 (2015) 230–238
20
20
15
15
y/L
y/L
232
10 5 0
10 5
0
5
10
0
15
0
5
t/T
20
20
15
15
10 5
10 5
0
5
10
0
15
0
5
t/T
20
15
15
y/L
y/L
(d) B. C1.
20
10 5
10 5
0
5
10
0
15
0
5
t/T
(f) B. C2.
20
20
15
15
y/L
y/L
15
10
t/T
(e) A. C2.
10 5 0
15
10
t/T
(c) A. C1.
0
15
(b) B. C1.
y/L
y/L
(a) A. C1.
0
10
t/T
10 5
0
5
10
15
0
5
0
t/T
10
15
t/T
(g) A. C2.
(h) B. C2.
Fig. 3. Re ¼ 50. Time history of the position of the centers of mass of A and B in the configurations C1 and C2 for different values of G, i.e.1.25 (red solid line), 2.5 (green solid line), 3.75 (blue solid line), 5 (magenta solid line), 6.25 (red dotted line), 7.5 (green dotted line), 8.75 (blue dotted line), 10 (magenta dotted line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
forced to move along nine prescribed directions i with velocities ci . This equation reads as follows:
f i ðx þ Dt ci ; t þ DtÞ ¼ f i ðx; tÞ þ
1
s
Mach number [34], the equilibrium particle distribution functions eq f i can be computed as
"
eq ½f i ðx; tÞ
f i ðx; tÞ þ 3wi c i gðx; tÞ; ð5Þ
where x is the position, Dt is the time step and s is the so-called relaxation parameter. It is worth to notice that the relaxation parameter s is strictly related to the fluid viscosity as P m ¼ ðs 1=2Þc2s , being c2s ¼ wi c2i ¼ 1=3 with wi a set of 9 weights defined as w0 ¼ 4=9; w1 ¼ w2 ¼ w3 ¼ w4 ¼ 1=9 and w5 ¼ w6 ¼ w7 ¼ w8 ¼ 1=36. By performing a second-order expansion in the local
eq f i ðx; tÞ
¼ wi q 1 þ
v ci c2s
# ðv c i Þ2 v v þ : 2c4s 2c2s
ð6Þ
The macroscopic fluid density q and the flow velocities v are defined as
qðx; tÞ ¼
X f i ðx; tÞ; i
v ðx; tÞ ¼
P
i f i ðx; tÞci
qðx; tÞ
;
ð7Þ
respectively. In order to account for the presence of the wings in the fluid lattice background, the immersed boundary method [38] is
233
20
20
15
15
y/L
y/L
A. De Rosis / Computers and Structures 153 (2015) 230–238
10 5
5
0
5
10
0
15
(b) B. C1.
15
15
10 5
0
5
10
0
15
0
5
10
t/T
t/T
(c) A. C1.
(d) B. C1. 20
15
15
y/L
y/L
5
10 5
15
10 5
0
5
10
0
15
0
5
10
t/T
t/T
(e) A. C2.
(f) B. C2.
20
20
15
15
y/L
y/L
15
10
20
10 5 0
10
(a) A. C1. 20
0
5
t/T
20
0
0
t/T
y/L
y/L
0
10
15
10 5
0
5
10
15
0
0
5
10
t/T
t/T
(g) A. C2.
(h) B. C2.
15
Fig. 4. Re ¼ 100. Time history of the position of the centers of mass of A and B in the configurations C1 and C2 for different values of G, i.e.1.25 (red solid line), 2.5 (green solid line), 3.75 (blue solid line), 5 (magenta solid line), 6.25 (red dotted line), 7.5 (green dotted line), 8.75 (blue dotted line), 10 (magenta dotted line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
adopted. Specifically, each wing is idealized by a set of Lagrangian solid points, namely X j , whose reciprocal spacing is set to a half of the fluid mesh size. The IB is adopted in an implicit velocitycorrection strategy, leading to a correction term gðx; tÞ which is used in the right-hand side of Eq. (5) [44,46]. It is worth to notice that fluid forces, which are used in Eq. (4), can be easily and immediately computed as
Fj ¼
X gðxÞWðx X j Þ;
dynamics. The adopted numerical methods are combined within a staggered explicit coupling algorithm. Its high properties in terms of stability and accuracy have already been widely tested by the author by computing the error in terms of the energy that is artificially introduced in the system at the fluid–structure interface [46,47]. Finally, notice that the problem is solved in LB units. This means that particular attention should be paid to the conversion from the physical system to the LB world.
ð8Þ
x
3. Results and discussion W being the IB kernel [38]. Once fluid forces acting upon the wings are computed, the solid position is updated by adopting the TDG method. The interested reader can refer to [45,57,58] for further details about such time integration scheme for predicting solid
In this section, results from a numerical campaign are discussed. Each wing is idealized by Llb ¼ 100 lattice nodes and it is represented by 200 Lagrangian immersed boundary points.
A. De Rosis / Computers and Structures 153 (2015) 230–238
20
20
15
15
y/L
y/L
234
10 5 0
10 5
5
0
10
0
15
5
0
10
t/T
(b) B. C1.
20
20
15
15
y/L
y/L
(a) A. C1.
10 5 0
10 5
0
5
10
0
15
0
5
10
t/T
20
15
15
y/L
y/L
(d) B. C1.
20
10 5
10 5
0
5
10
0
15
0
5
10
t/T
(f) B. C2.
20
20
15
15
y/L
y/L
15
t/T
(e) A. C2.
10 5 0
15
t/T
(c) A. C1.
0
15
t/T
10 5
0
5
10
0
15
0
5
10
t/T
15
t/T
(g) A. C2.
(h) B. C2.
10
8
8
6
6
yt /L
yt /L
Fig. 5. Re ¼ 200. Time history of the position of the center of mass of A and B in the configurations C1 and C2 for different values of G, i.e.1.25 (red solid line), 2.5 (green solid line), 3.75 (blue solid line), 5 (magenta solid line), 6.25 (red dotted line), 7.5 (green dotted line), 8.75 (blue dotted line), 10 (magenta dotted line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4
2
2 0
4
2.5
7.5
5
10
0
2.5
5
7.5
g/L
g/L
(a) A.
(b) B.
10
Fig. 6. Position achieved at t=T ¼ 6, namely yt =L, by A and B vs the gap-over-length ratio, g=L, for different values of the Reynolds number, i.e. Re ¼ 50 (red lines), 100 (green) and 200 (blue). Configurations C1 and C2 are represented by solid and dotted lines, respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
235
A. De Rosis / Computers and Structures 153 (2015) 230–238 Table 1 A. Position achieved at t=T ¼ 6 for different values of G and Re in the configurations C1 and C2. The peak value is highlighted by bold. G
1.25 2.5 3.75 5 6.25 7.5 8.75 10
Re ¼ 50
Re ¼ 100
Re ¼ 200
C1
C2
C1
C2
C1
C2
8.6353 5.1355 3.4593 2.5957 2.2740 2.2103 2.2401 2.3745
3.0715 1.9865 1.6076 1.3911 1.2902 1.2785 1.2955 1.3176
8.9457 3.0395 1.9508 1.3675 1.0388 0.9179 0.9777 1.1439
1.9823 1.3662 1.2304 1.0391 0.8819 0.8615 0.9162 1.0105
9.8773 4.8161 2.9350 1.8429 1.3233 1.0071 1.0805 1.1974
1.0885 0.6463 0.3788 0.3296 0.2858 0.1156 0.1639 0.2214
Table 2 B. Position achieved at t=T ¼ 6 for different values of G and Re in the configurations C1 and C2. The peak value is highlighted by bold. G
1.25 2.5 3.75 5 6.25 7.5 8.75 10
Re ¼ 50
Re ¼ 100
Re ¼ 200
C1
C2
C1
C2
C1
C2
7.7949 4.4192 3.0824 2.3642 2.0768 2.0117 2.0430 2.1779
3.0334 1.9457 1.5816 1.3712 1.2682 1.2528 1.2725 1.3014
5.7789 2.6741 1.8413 1.3442 1.0850 1.0060 1.0359 1.1925
1.9365 1.3516 1.2287 1.0473 0.9001 0.8820 0.9310 1.0173
5.3232 2.9726 1.8029 1.4218 1.1144 0.8423 0.9508 1.1802
0.7726 0.3873 0.2639 0.2581 0.2029 0.1117 0.1116 0.1628
(a) t/T = 0.76.
(b) t/T = 2.29.
(c) t/T = 6.89.
(d) t/T = 8.42.
Moreover, the grid consists of 1000 and 2000 lattice nodes in the horizontal and vertical directions, respectively. All the LB simulations are carried out at low Mach number Ma, i.e. Ma ¼ v tip;lb =cs ¼ 0:0173, where v tip;lb ¼ 4Llb Dh=T lb ¼ 0:01. Moreover, the initial density is set to q ¼ 1. In order to help the reader in the unit conversion, the author defines the lattice grid spacing as
Dx ¼ L=Llb ¼ 3 104 m and the time step Dt ¼ 3:061 106 s. Therefore, the period of oscillation in the LB simulations is set to T lb ¼ T=Dt ¼ 32672:56 lattice time steps. Simulations stop when A or B overcomes the fluid domain limits or if t=T > 15. Scenarios characterized by different values of the relaxation parameter s are investigated, thus corresponding to distinct values of the Reynolds number, i.e. Re ¼ v tip;lb Llb =m ¼ 50; 100; 200. Moreover, two configurations are dissected. In the former (C1), each center of mass is far 100 lattice nodes from the vertical axis of symmetry, whereas in the latter (C2) such distance is set to 200 lattice nodes. At the bottom-most section the no-slip condition is enforced, i.e. v ¼ 0. At the remaining boundaries, outflow conditions are prescribed, i.e. n rf i ¼ 0; n being the outer normal to the boundary. In Fig. 3, the time history of the position of the centers of mass of A and B are depicted in the configurations C1 and C2 for different values of G, i.e. 1.25, 2.5, 3.75, 5, 6.25, 7.5, 8.75, 10. As it is possible to observe, two different patterns arise. The former collects the scenarios characterized by G 6 5 (solid lines), where the lower G, the faster the take-off is. At a first glance, such behavior can be confusing. In fact, A and B achieve higher values of the position if these start from a lower initial abscissa. On the other hand, if G P 5 (dotted lines), such trend drastically changes, demonstrating that the flight performance are remarkably dominated by the initial position. Indeed, in the latter pattern, the higher the initial position, the higher the final position is. This dramatic variation has to be addressed to the role of the solid surface which is located at the bottom-most section. In fact, as the wings are close to the boundary, these exhibit a progressively enhancement in terms of generated lift force, thus corresponding to higher values of the achieved position. As demonstrated in [32], if a solid body is immersed in a
15 (e) Color legend. Fig. 7. Most advantageous scenarios: velocity field at different time instants for Re ¼ 200 and G ¼ 1:25 in C1.
viscous fluid and it oscillates close to a solid surface, a growth of the fluid force arises, consisting of one order of magnitude with respect to the case in which the solid surface is neglected. Therefore, it is possible to assess that the flapping wings considerably benefit from the presence of this surface. Such behavior is known as the so-called ground effect. Moreover, the behavior is even dominated by the reciprocal distance between A and B. In fact, the wings appear to benefit from being relatively close. This
236
A. De Rosis / Computers and Structures 153 (2015) 230–238
(a) t/T = 0.76.
(b) t/T = 2.29.
reported in Fig. 5. Again, the deleterious effects experienced at G > 5 are still present. On the other hand, configurations C1 and C2 show different trends. In the former, curves appears to be shifted back in time, whereas the latter presents a time shift in advance with respect to the scenarios characterized by Re ¼ 100. This means that the solution of the problem is highly non-linear, since it is strongly dependent on G; Re and the horizontal reciprocal distance. Different patterns arose due to the distinct tip-induced vorticity which involves the flow field. As well known, it is related to the fluid viscosity and to the distance between solid surfaces. In order to summarize the above mentioned effects, in Fig. 6 the position achieved at t=T ¼ 6, namely yt =L, by A and B is plotted against the gap-over-length ratio, G ¼ g=L, for different values of the Reynolds number in the configurations C1 and C2. As it is possible to observe, scenarios corresponding to the configuration C1 are considerably higher with respect to the ones related to C2. As G decreases, the benefit is higher, whereas curves plateau for progressively larger values of the gap-over-length ratio. In addition, notice that a more marked advantage is experienced by A, since B shows lower values of yt =L. In order to help the reader, data depicted in Fig. 6 are reported in Tables 1 and 2 for A and B, respectively. Finally, making reference to the flight performance of A, in Figs. 7 and 8 the velocity field is depicted at different time instants in the most advantageous configuration, i.e Re ¼ 200 and G ¼ 1:25 in C1, and the most arduous one, Re ¼ 200 and G ¼ 7:5 in C2. The above mentioned considerations about the role of the reciprocal distance are confirmed. In particular, notice that if the wings are closer, a more marked vorticity arises in the flow field. Vortexes are mainly located below the wing, thus considerably contributing to the lift force enhancement. On the other hand, the scenario manifesting an arduous take-off is characterized by low values of the velocity. 4. Conclusions
(c) t/T = 6.89.
(d) t/T = 8.42.
16 (e) Color legend. Fig. 8. Most arduous scenarios: velocity field at different time instants for Re ¼ 200 and G ¼ 7:5 in C2.
trend has to be addressed to the tip-induced vorticity, which highly supports wings traveling. Configurations corresponding to Re ¼ 100 are shown in Fig. 4. Notice that curves appear to be shifted in advance in time with respect to the previous case. The threshold represented by G 5 is already experienced, since two different trend in the flight performance are found. In addition, the configuration C2 exhibits a remarkably arduous take-off, except for the peak value of G. Finally, scenarios at Re ¼ 200 are
In this paper, the dynamics involving a tandem of symmetric flapping wings immersed in a viscous fluid near a solid surface has been investigated. Specifically, the dependence of the flight performance in terms of time history of the position of the centers of mass on the fluid viscosity, i.e. the Reynolds number Re, has been elucidated. Moreover, the so-called ground effect has been dissected by performing simulations characterized by progressively higher values of the gap-over-length ratio, G. The author demonstrated that below a certain threshold, i.e. G 6 5, wings dynamics is positively affected by the presence of the solid surface. In fact, a remarkable enhancement of the flight performance has been found. If G > 5, the advantageous ground effect vanishes. Finally, the incidence of wings reciprocal distance has been elucidated, thus underlining the mutual interaction between A and B, especially in terms of tip-induced vorticity. It has been demonstrated that wings benefit from a closer arrangement, thus justifying the natural attitude to move in flocks. Finally, it is possible to assess that the above adopted numerical algorithms can be successfully employed to solve problems involving the dynamics of flapping wings which are immersed in a viscous fluid. Acknowledgement The author would like to thank Prof. Stefano Ubertini and Prof. Francesco Ubertini for the useful and motivating discussions. Appendix A The convergence and accuracy properties of the adopted numerical methods are evaluated in the following convergence
A. De Rosis / Computers and Structures 153 (2015) 230–238
Error
-2
-3
-4
5
5.5
6
6.5
Lattice nodes Fig. 9. Convergence analysis: relative error vs number of lattice nodes in a log–log scale.
analysis. Specifically, the value of yt =L at Re ¼ 200 and G ¼ 1:25 in C1 for A is computed by progressively refining the lattice fluid grid. For a given grid resolution, such value is denoted ecurr . It is compared to the same quantity computed on a very fine grid composed of 4000 2000 lattice nodes and it is denoted eref . Then, it is possible to define the relative error as
Error ¼
eref ecurr : eref
ð9Þ
In Fig. 9, the relative error is depicted against the number of lattice nodes. Both the quantities are given in a logarithmic scale. As it is possible to observe, a good convergence rate equal to 1.78 is experienced, thus allowing us to assess that the proposed strategy shows high properties in terms of accuracy and convergence. References [1] Ellington CP. The novel aerodynamics of insect flight: applications to micro-air vehicles. J Exp Biol 1999;202:3439–48. [2] Shyy W, Berg M, Ljungqvist D. Flapping and flexible wings for biological and micro air vehicles. Prog Aerosp Sci 1999;35:455–505. [3] Ettinger SM, Nechyba MC, Ifju PG, Waszak M. Vision-guided flight stability and control for micro air vehicles. Adv Robot 2003;17:617–40. [4] Ansari S, Abikowski R, Knowles K. Aerodynamic modelling of insect-like flapping flight for micro air vehicles. Prog Aerosp Sci 2006;42:129–72. [5] Shyy W, Aono H, Chimakurthi SK, Trizila P, Kang C-K, Cesnik CE, et al. Recent progress in flapping wing aerodynamics and aeroelasticity. Prog Aerosp Sci 2010;46:284–327. [6] Orlowski CT, Girard AR. Dynamics, stability, and control analyses of flapping wing micro-air vehicles. Prog Aerosp Sci 2012;51:18–30. [7] Aureli M, Prince C, Porfiri M, Peterson SD. Energy harvesting from base excitation of ionic polymer metal composites in fluid environments. Smart Mater Struct 2010;19:015003. [8] Giacomello A, Porfiri M. Underwater energy harvesting from a heavy flag hosting ionic polymer metal composites. J Appl Phys 2011;109:084903. [9] Erturk A, Delporte G. Underwater thrust and power generation using flexible piezoelectric composites: an experimental investigation toward self-powered swimmer-sensor platforms. Smart Mater Struct 2011;20:125013. [10] Rozhdestvensky KV, Ryzhov VA. Aerohydrodynamics of flapping-wing propulsors. Prog Aerosp Sci 2003;39:585–633. [11] Trizila P, Kang C-K, Aono H, Shyy W, Visbal M. Low-reynolds-number aerodynamics of a flapping rigid flat plate. AIAA J 2011;49:806–23. [12] Fuchiwaki M, Kuroki T, Tanaka K, Tababa T. Dynamic behavior of the vortex ring formed on a butterfly wing. Exp Fluids 2013;54:1–12. [13] De Rosis A. On the dynamics of a tandem of asynchronous flapping wings: lattice Boltzmann-immersed boundary simulations. Physica A 2014;410:276–86. [14] Heathcote S, Gursul I. Flexible flapping airfoil propulsion at low reynolds numbers. AIAA J 2007;45:1066–79. [15] Heathcote S, Wang Z, Gursul I. Effect of spanwise flexibility on flapping wing propulsion. J Fluids Struct 2008;24:183–99. [16] Majumder L, Rao S. Interval-based optimization of aircraft wings under landing loads. Comput Struct 2009;87:225–35. [17] Masoud H, Alexeev A. Resonance of flexible flapping wings at low reynolds number. Phys Rev E 2010;81:056304. [18] Mountcastle AM, Combes SA. Wing flexibility enhances load-lifting capacity in bumblebees. Proc R Soc B 2013;280. [19] Kang C-k, Shyy W. Scaling law and enhancement of lift generation of an insectsize hovering flexible wing. J R Soc Interface 2013;10.
237
[20] De Rosis A, Falcucci G, Ubertini S, Ubertini F. Aeroelastic study of flexible flapping wings by a coupled lattice Boltzmann finite element approach with immersed boundary method. J Fluids Struct 2014;49:516–33. [21] Kim T-U, Hwang IH. Optimal design of composite wing subjected to gust loads. Comput Struct 2005;83:1546–54. [22] Kameyama M, Fukunaga H. Optimum design of composite plate wings for aeroelastic characteristics using lamination parameters. Comput Struct 2007;85:213–24. [23] Kim T-U, Shin JW, Hwang IH. Stacking sequence design of a composite wing under a random gust using a genetic algorithm. Comput Struct 2007;85:579–85. Stochastic Structural Analysis, Optimization and Re-Analysi. [24] Jones K, Platzer M. Numerical computation of flapping-wing propulsion and power extraction. AIAA Paper 1997;97:0826. [25] Zerihan J, Zhang X. Aerodynamics of a single element wing in ground effect. J Aircraft 2000;37:1058–64. [26] Zhang X, Zerihan J. Aerodynamics of a double-element wing in ground effect. AIAA J 2003;41:1007–16. [27] Zerihan J, Zhang X. Aerodynamics of gurney flaps on a wing in ground effect. AIAA J 2001;39:772–80. [28] Rozhdestvensky KV. Wing-in-ground effect vehicles. Prog Aerosp Sci 2006;42:211–83. [29] Ahmed N, Goonaratne J. Lift augmentation of a low-aspect-ratio thick wing in ground effect. J Aircraft 2002;39:381–4. [30] Green CP, Sader JE. Frequency response of cantilever beams immersed in viscous fluids near a solid surface with applications to the atomic force microscope. J Appl Phys 2005;98:114913. [31] Grimaldi E, Porfiri M, Soria L. Finite amplitude vibrations of a sharp-edged beam immersed in a viscous fluid near a solid surface. J Appl Phys 2012;112:104907. [32] De Rosis A. Harmonic oscillations of a lamina in a viscous fluid near a solid surface: a lattice Boltzmann-immersed boundary approach. Physica A 2014;415:386–97. [33] Higuera FJ, Succi S, Benzi R. Lattice gas dynamics with enhanced collisions. Europhys Lett 1989;9:345–9. [34] Benzi R, Succi S, Vergassola M. The lattice Boltzmann equation: theory and applications. Phys Rep 1992;222:145–97. [35] Succi S. The lattice Boltzmann equation for fluid dynamics and beyond. Clarendon; 2001. [36] Chen H, Chen S, Matthaeus WH. Recovery of the Navier–Stokes equations using a lattice-gas Boltzmann method. Phys Rev Lett 1992;45:R5339–42. [37] Fadlun E, Verzicco R, Orlandi P, Mohd-Yusof J. Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. J Comput Phys 2000;161:35–60. [38] Peskin C. The immersed boundary method. Acta Numer 2002;11:479–517. [39] Wu J, Shu C. Implicit velocity correction-based immersed boundarylattice Boltzmann method and its applications. J Comput Phys 2009;228: 1963–79. [40] Sotiropoulos F, Yang X. Immersed boundary methods for simulating fluid– structure interaction. Prog Aerosp Sci 2014;65:1–21. [41] Filippova O, Hänel D. Lattice Boltzmann simulation of gas-particle flow in filters. Comput Fluids 1997;26:697–712. [42] Mei R, Luo L, Shyy W. An accurate curved boundary treatment in the lattice Boltzmann method. J Comput Phys 1999;155:307–30. [43] Mei R, Yu D, Shyy W, Luo L. Force evaluation in the lattice Boltzmann method involving curved geometry. Phys Rev Lett E 2002;65:041203. [44] De Rosis A, Ubertini S, Ubertini F. A comparison between the interpolated bounce-back scheme and the immersed boundary method to treat solid boundary conditions for laminar flows in the lattice Boltzmann framework. J Sci Comput 2014;61:477–89. [45] Mancuso M, Ubertini F. An efficient integration procedure for linear dynamics based on a time discontinuous Galerkin formulation. Comput Mech 2003;32:154–68. [46] De Rosis A, Ubertini S, Ubertini F. A partitioned approach for two-dimensional fluid–structure interaction problems by a coupled lattice Boltzmannfinite element method with immersed boundary. J Fluids Struct 2014;45: 202–15. [47] De Rosis A. A lattice Boltzmann-finite element model for two-dimensional fluid–structure interaction problems involving shallow waters. Adv Water Resour 2014;65:18–24. [48] De Rosis A. A lattice Boltzmann model for multiphase flows interacting with deformable bodies. Adv Water Resour 2014;73:55–64. [49] De Rosis A, Falcucci G, Porfiri M, Ubertini F, Ubertini S. Hydroelastic analysis of hull slamming coupling lattice Boltzmann and finite element methods. Comput Struct 2014;138:24–35. [50] De Rosis A. Harmonic oscillations of laminae in non-Newtonian fluids: a lattice Boltzmann-immersed boundary approach. Adv Water Resour 2014;73:97–107. [51] Srinivas Rao K, Subrahmanya M, Kulkarni D, Rajani B. Numerical simulation of turbulent flow past micro air vehicle wing sections. Bangalore: National Aerospace Laboratories; 2009. [52] Ota K, Suzuki K, Inamuro T. Lift generation by a two-dimensional symmetric flapping wing: immersed boundary-lattice Boltzmann simulations. Fluid Dyn Res 2012;44:045504. [53] Inamuro T, Kimura Y, Suzuki K. Flight simulations of a two-dimensional flapping wing by the ib-lbm. Bull Am Phys Soc 2012;57.
238
A. De Rosis / Computers and Structures 153 (2015) 230–238
[54] Kimura Y, Suzuki K, Inamuro T. Flight simulations of a two-dimensional flapping wing by the ib-lbm. Int J Mod Phys C 2014;25. [55] Sunada S, Kawachi K, Watanabe I, Azuma A. Performance of a butterfly in takeoff flight. J Exp Biol 1993;183:249–77. [56] Iima T, Yanagita M. Asymmetric motion of a two-dimensional symmetric flapping model. Fluid Dyn Res 2005;36:407–25.
[57] Mancuso M, Ubertini F. An efficient time discontinuous Galerkin procedure for non-linear structural dynamics. Comput Methods Appl Mech Eng 2006;195: 6391–406. [58] De Rosis A, Falcucci G, Ubertini S, Ubertini F, Succi S. Lattice Boltzmann analysis of fluid–structure interaction with moving boundaries. Commun Comput Phys 2012;13:823–34.