4. Grid dependence study and vortical structures visualization

When conducting numerical simulations, it is well known that the grid resolution can affect the results from a quantitative and qualitative point of view. In this section, we present the outcome of the grid dependence study used to determine the best overlapping grid system G in terms of computing time, stability, and accuracy of the solution. To conduct this study, we used the grid convergence index method or GCI, as described by Roache in [25, 26].

During this study, fixed and moving wings were considered, but for simplicity, we will only present the results related to pure heaving motion (where the uncertainties are higher). Several simulations were run at a Strouhal number St = 0.3 and at a reduced frequency k = 1.570795, and unsteadiness was observed to disappear typically after 4 to 5 cycles of wing heaving motion, and further calculations show negligible nonperiodicity. Each simulation was checked for acceptable iterative convergence.

The different grid sizes, layouts, and mesh stretching ratios considered during this study are shown in Table 1. In Table 1, the grid spacing ratio (GSR) is the refinement ratio from the finer grid to the coarser grid and is expressed in reference to the position of the first node normal to the wing surface (1NW). Therefore, the grid spacing refinement ratio r is equal to 2. In Figure 2, we illustrate a typical overlapping grid system G used in this study. In Figure 2, the background grid is shown in images a, b, and c (rectangular grid). The wing grid is made up of


BG stands for background grid, WG-CS stands for wing grid center-section, WG-TS stands for wing grid tip-section (only one tip-section), GSR stands for grid refinement spacing ratio, and 1NW stands for position of the first node normal to the wall.

In all the cases studied, a rectangular wing with an aspect ratio AR equal to 2 was used. The cross-section of the wing is an ellipse, with a corresponding major axis a = 0.25 and a minor axis

The initial conditions used for all the heaving wings simulations are those of a fully converged solution of the corresponding fixed wing case. In Figure 1 (left), the left boundary of the BG corresponds to an inflow boundary condition (u = (1.0, 0.0, 0.0), ∂n^ p ¼ 0), and the top, bottom, and right boundaries of the BG are outflow boundaries (velocity extrapolated from the interior points). In Figure 1 (right), all the boundaries of the BG correspond to outflow conditions. On the wing surface (which is a moving body), we impose a no-slip boundary condition for

we used a nonconservative Lagrange interpolation scheme. The Reynolds number (defined as

In all the simulations conducted, we assumed that the wing is undergoing pure heaving motion, wherein the wing cross-section heaves in the vertical direction (or y axis in Figure 1)

where y(t) is the heaving motion (and is defined positive upwards), h is the heaving amplitude, f is the heaving oscillating frequency, φ is the phase angle of the heaving motion (0 in this case),

When conducting numerical simulations, it is well known that the grid resolution can affect the results from a quantitative and qualitative point of view. In this section, we present the outcome of the grid dependence study used to determine the best overlapping grid system G in terms of computing time, stability, and accuracy of the solution. To conduct this study, we

During this study, fixed and moving wings were considered, but for simplicity, we will only present the results related to pure heaving motion (where the uncertainties are higher). Several simulations were run at a Strouhal number St = 0.3 and at a reduced frequency k = 1.570795, and unsteadiness was observed to disappear typically after 4 to 5 cycles of wing heaving motion, and further calculations show negligible nonperiodicity. Each simulation was checked

The different grid sizes, layouts, and mesh stretching ratios considered during this study are shown in Table 1. In Table 1, the grid spacing ratio (GSR) is the refinement ratio from the finer grid to the coarser grid and is expressed in reference to the position of the first node normal to the wing surface (1NW). Therefore, the grid spacing refinement ratio r is equal to 2. In Figure 2, we illustrate a typical overlapping grid system G used in this study. In Figure 2, the background grid is shown in images a, b, and c (rectangular grid). The wing grid is made up of

4. Grid dependence study and vortical structures visualization

used the grid convergence index method or GCI, as described by Roache in [25, 26].

<sup>z</sup>). The rest of the boundaries is interpolation boundaries, where

y tðÞ¼ <sup>h</sup> � sin 2 � <sup>π</sup> � <sup>f</sup> � <sup>t</sup> <sup>þ</sup> <sup>φ</sup> (15)

b = 0.025. Therefore, the wing's chord c is equal to 2 � a = 0.5.

moving walls (u ¼ G

and t is the physical time.

for acceptable iterative convergence.

� x, G � y, G �

124 Flight Physics - Models, Techniques and Technologies

and according to the following function,

Re = U � c/ν) is equal to 250 for all the simulations.

Figure 2. Typical overlapping grid system G used in this study. (A) Side view. (B) Front view. (C) Bottom view. (D) Perspective view.

three component grids (image d), namely, center section, left wing grid tip-section, and right wing grid tip-section.

Since in the study of heaving wings propulsion, the main task is thrust production, and it is more convenient to think in terms of thrust force T instead of drag force D. The thrust force is equal in magnitude but opposite in direction to the drag force (T = � D). To quantify the unsteady aerodynamics performance, we computed the lift coefficient cl, the drag coefficient cd, and thrust coefficient ct as follows,

$$c\_l = \frac{L}{\frac{1}{2}\rho lI^2A}, \quad c\_d = \frac{D}{\frac{1}{2}\rho lI^2A}, \quad c\_l = \frac{T}{\frac{1}{2}\rho lI^2A} \tag{16}$$

In Eq. (16), the lift force L and the thrust force T (where T = � D) are computed by integrating the viscous and pressure forces over the wing surface.


Table 2. Observed values of cd and crms <sup>l</sup> for the grid dependence study.

In Table 2, we present the average drag coefficient cd and the root mean square of the lift coefficient crms <sup>l</sup> for the grid dependence study. These values were used to compute the observed order of convergence pobs, the extrapolated value of the observed quantity at zero grid spacing χspacing = 0, the fine grid convergence index GCI12 and GCI23, and the constancy of GCI<sup>23</sup> = r pobs GCI12.

Based on the outcome of the GCI study (refer to Table 3), cd at zero grid spacing is estimated to be 0.057708, with an error band of 0.397203% for grid G<sup>1</sup> and an error band of 1.641617% for grid G2, whereas crms <sup>l</sup> at zero grid spacing is estimated to be 0.908738 with an error band of 0.402503% for grid G<sup>1</sup> and an error band of 1.567203% for grid G2. The constancy of GCI<sup>23</sup> = r pobs GCI<sup>12</sup> indicates that the results for grids <sup>G</sup><sup>1</sup> and <sup>G</sup><sup>2</sup> are within the asymptotic range of convergence. Finally, the observed order of convergence pobs represents a direct indication of the accuracy of the numerical method. In this study, we obtained a value of pobs close to 2 for both quantities of interest, which indicate that the method is second-order accurate in space and time.

To supplement the quantitative GCI study, we conducted an additional grid dependence study but from the qualitative point of view (vortical structures resolution on the wake of the wing). The qualitative results for grids G<sup>1</sup> and G<sup>2</sup> (refer to Table 1) are illustrated in Figure 3. In Figure 3, we use the iso-surfaces of vorticity magnitude (∣ω∣-criterion), and the iso-surfaces of Q-criterion [27, 28], to capture the vortices and their corresponding cores. As depicted in this figure, there are no discernible differences between the solutions, indicating that both grids are adequate for vortical structures resolution. It can also be seen that the ∣ω∣-criterion, although capable of capturing the general vortical structures, has the disadvantage of also showing the shear layers near the wing surface and between the vortices. On the other hand, the Q-criterion shows the details of the vortical structures more clearly as it does not show the shear layers close to the wing (as illustrated in Figure 4). Apart from this, both methods provide nearly


Table 3. GCI study results.

In Table 2, we present the average drag coefficient cd and the root mean square of the lift

<sup>l</sup> for the grid dependence study.

Grid Gg<sup>g</sup> GSR cd crms

G<sup>1</sup> 1 0.057892 0.905822 G<sup>2</sup> 2 0.058476 0.897486 G<sup>3</sup> 4 0.060914 0.865326

observed order of convergence pobs, the extrapolated value of the observed quantity at zero grid spacing χspacing = 0, the fine grid convergence index GCI12 and GCI23, and the constancy

Based on the outcome of the GCI study (refer to Table 3), cd at zero grid spacing is estimated to be 0.057708, with an error band of 0.397203% for grid G<sup>1</sup> and an error band of 1.641617% for

0.402503% for grid G<sup>1</sup> and an error band of 1.567203% for grid G2. The constancy of

range of convergence. Finally, the observed order of convergence pobs represents a direct indication of the accuracy of the numerical method. In this study, we obtained a value of pobs close to 2 for both quantities of interest, which indicate that the method is second-order

To supplement the quantitative GCI study, we conducted an additional grid dependence study but from the qualitative point of view (vortical structures resolution on the wake of the wing). The qualitative results for grids G<sup>1</sup> and G<sup>2</sup> (refer to Table 1) are illustrated in Figure 3. In Figure 3, we use the iso-surfaces of vorticity magnitude (∣ω∣-criterion), and the iso-surfaces of Q-criterion [27, 28], to capture the vortices and their corresponding cores. As depicted in this figure, there are no discernible differences between the solutions, indicating that both grids are adequate for vortical structures resolution. It can also be seen that the ∣ω∣-criterion, although capable of capturing the general vortical structures, has the disadvantage of also showing the shear layers near the wing surface and between the vortices. On the other hand, the Q-criterion shows the details of the vortical structures more clearly as it does not show the shear layers close to the wing (as illustrated in Figure 4). Apart from this, both methods provide nearly

Outcome cd crms

pobs 2.061650 1.947838 χspacing = 0 0.057708 0.908738 GCI12 (%) 0.397203 0.402503 GCI23 (%) 1.641617 1.567203 (GCI23/GCI12) (1/rpobs) 0.990012 1.009249

pobs GCI<sup>12</sup> indicates that the results for grids <sup>G</sup><sup>1</sup> and <sup>G</sup><sup>2</sup> are within the asymptotic

<sup>l</sup> for the grid dependence study. These values were used to compute the

l

l

<sup>l</sup> at zero grid spacing is estimated to be 0.908738 with an error band of

coefficient crms

of GCI<sup>23</sup> = r

GCI<sup>23</sup> = r

grid G2, whereas crms

accurate in space and time.

Table 3. GCI study results.

pobs GCI12.

Table 2. Observed values of cd and crms

126 Flight Physics - Models, Techniques and Technologies

Figure 3. Vortex topology at the beginning of the upstroke (t = 5.0). Iso-surfaces of ∣ω∣-criterion are shown in light gray, and iso-surfaces of Q-criterion are shown in dark gray. Simulation parameters: Re = 250, St = 0.3, and k = 1.570795. A corresponds to overlapping grid system G<sup>1</sup> and B corresponds to overlapping grid system G2.

identical structures in the far wake. Based on these qualitative results, we chose the Q-criterion as the main criterion for wake topology characterization.

Summarizing the quantitative and qualitative results previously presented, we can conclude that the solutions obtained by using the overlapping grid systems G<sup>1</sup> and G<sup>2</sup> are grid independent. Taking into account the computational resources available, CPU time restrictions, and

Figure 4. Shear layers close to the wing at the beginning of the upstroke (t = 5.0). Left column corresponds to iso-surfaces of Q-criterion, and the right column corresponds to iso-surfaces of ∣ω∣-criterion. Simulation parameters: Re = 250, St = 0.3, and k = 1.570795.

solution accuracy, G<sup>2</sup> with a value of 1 NW equal to 0.001 2c is used as the base grid to perform all further computations. In the case of a smaller or bigger computational domain, the grid dimensions are scaled in order to keep the same grid spacing as for this domain.

#### 5. Simulation results

Hereafter, we carry out a comprehensive parametric study to assess the wake signature and aerodynamic performance of heaving rigid wings. In Table 4, we present the kinematics

Wake Topology and Aerodynamic Performance of Heaving Wings http://dx.doi.org/10.5772/intechopen.71517 129


Table 4. Simulation results for the pure heaving parametric study (positive ct indicates thrust production whereas negative ct indicates drag production).

parameters governing the heaving motion (described by Eq. (15)). In Table 4, we also present the quantitative results obtained, where ct is the average thrust coefficient and <sup>b</sup>cl is the maximum lift coefficient (which was measured during the downstroke). By inspecting these results, we observe that as we increase St and <sup>k</sup>, the values of ct and <sup>b</sup>cl also increase. This result is expected because as we increase St and k (therefore, the oscillating frequency and heaving amplitude), the vertical velocity of the wing is higher. Hence, the forces excerpted on the wing's surface are larger. We also observe two different behaviors of the aerodynamic forces for high and low reduced frequencies k values. Hence, it seems that for heaving wings, the oscillating frequency plays an important role in the vortex generation and shedding and, henceforth, on the aerodynamic forces. These frequency dependence observations are similar to those of Wang [29], Young and Lai [30], and Guerrero [31], but here, we extend them to three-dimensional cases.

In Table 4, we can read that for values of St < 0.30, the heaving wing produces drag; for values of 0.30 < St < 0.35, the heaving wing produces little or no drag (or thrust), whereas for values of St > 0.35, the heaving wing produces thrust. These results suggest that there is a range of St values, where the heaving wing generates thrust. The results also point at the presence of a combination of St and k values, where propulsive efficiency peaks. In general, the results are inline with the hypothesis that flying and swimming animals cruise at Strouhal numbers corresponding to a regime of vortex growth and shedding in which propulsion efficiency is high.

solution accuracy, G<sup>2</sup> with a value of 1 NW equal to 0.001 2c is used as the base grid to perform all further computations. In the case of a smaller or bigger computational domain, the

Figure 4. Shear layers close to the wing at the beginning of the upstroke (t = 5.0). Left column corresponds to iso-surfaces of Q-criterion, and the right column corresponds to iso-surfaces of ∣ω∣-criterion. Simulation parameters: Re = 250, St = 0.3,

Hereafter, we carry out a comprehensive parametric study to assess the wake signature and aerodynamic performance of heaving rigid wings. In Table 4, we present the kinematics

grid dimensions are scaled in order to keep the same grid spacing as for this domain.

5. Simulation results

128 Flight Physics - Models, Techniques and Technologies

and k = 1.570795.

To understand how the heaving wing generates thrust, we need to take a closer look at the vortex shedding process. In Figures 5 and 6, we illustrate the vortex topology for a heaving case in the thrust producing regime (case 3DH-11 in Table 4). These figures show that the downstream wake (or far field wake) of this case consists of two sets of doughnut-shaped vortex rings (VR) which convect at oblique angles about the centerline of the heaving motion. Thus, the flow induced by each vortex ring along its axis is expected to have a net streamwise component linked to thrust production. This net streamwise momentum excess in the wake is connected with the thrust production of heaving wings. The process by which the vortex rings are formed can be explained by examining the vortex formation and shedding close to the

Figure 5. Vortex wake topology at the beginning of the upstroke for case 3DH-11 (t = 5.0). Heaving parameters: Re = 250, St = 0.4, and k = 1.570795 (thrust producing wake).

To understand how the heaving wing generates thrust, we need to take a closer look at the vortex shedding process. In Figures 5 and 6, we illustrate the vortex topology for a heaving case in the thrust producing regime (case 3DH-11 in Table 4). These figures show that the downstream wake (or far field wake) of this case consists of two sets of doughnut-shaped vortex rings (VR) which convect at oblique angles about the centerline of the heaving motion. Thus, the flow induced by each vortex ring along its axis is expected to have a net streamwise component linked to thrust production. This net streamwise momentum excess in the wake is connected with the thrust production of heaving wings. The process by which the vortex rings are formed can be explained by examining the vortex formation and shedding close to the

Figure 5. Vortex wake topology at the beginning of the upstroke for case 3DH-11 (t = 5.0). Heaving parameters: Re = 250,

St = 0.4, and k = 1.570795 (thrust producing wake).

130 Flight Physics - Models, Techniques and Technologies

Figure 6. Vortex wake topology at the beginning of the upstroke for case 3DH-11 (t = 5.0). Heaving parameters: Re = 250, St = 0.4, and k = 1.570795 (thrust producing wake).

wing's surface. During the heaving motion and close to the wing, four vortices are formed; namely, one leading edge vortex (LEV), one trailing edge vortex (TEV), and two wing-tip vortices or WTV (one WTV on the left wing-tip or LH-WTV and another WTV on the right wing-tip or RH-WTV). These four vortices are all connected and form a closed vortex loop (VL). During the heaving motion and as this vortex loop is convected downstream, it disconnects from the wing, creating in this way the doughnut-shaped vortex rings. It is also of interest in the fact that each vortex loop has two sets of thin contrails (TC) that connect the VL to the VR generated in the previous stroke; these structures are segments of the wing-tip vortices, and, as the vortex rings are convected downstream, they become weaker and ultimately disappear (as in vortex ring VR2 in Figure 5). During a complete heaving cycle, two VRs are formed, one at the end of the upstroke and the other one at the end of the downstroke.

In Figure 7, we present the wake topology for a drag-producing case (which corresponds to case 3DH-3 in Table 4). It is clear from this figure that the wake topology is very different from the one of the thrust producing case. In this case, as the vortex loops are convected downstream, they do not morph into vortex rings. Instead, they keep their original shape, and as they are convected, they diffuse. We can also observe that the wake height is very compact, in comparison to that of the thrust production case (as depicted in Figure 8). Finally, notice how the flow induced by each vortex loop is inclined in the same direction of the wing's travel direction, resulting this in a momentum surfeit linked to drag production. The momentum

Figure 7. Vortex wake topology at the beginning of the upstroke for case 3DH-3 (t = 5.0). Heaving parameters: Re = 250, St = 0.2, and k = 1.570795 (drag-producing wake).

Figure 8. Wake comparison of a drag-producing case (A) and thrust producing case (B). The wake height was measured approximately at 3 c behind the trailing edge.

deficit and momentum excess scenarios can be better appreciated in Figure 8, where A corresponds to drag production (momentum deficit) and B corresponds to trust production (momentum excess).

Animations of several simulation cases are available at the author's web site: http://www.dicat. unige.it/guerrero/flapsim1.html (last accessed: Sep 2017).
