Propeller Slipstream Effect on Aerodynamic Characteristics of Micro Air Vehicle at Low Reynolds Number

Abstract

A numerical investigation on propeller-induced flow effects in tractor configurations on a Zimmerman wing-fuselage using the cambered thin airfoil is presented in this paper. The Reynolds number based on the mean aerodynamic chord was 1.3 × 10⁵. Significant aerodynamic performance benefits could be found for a propeller in the tractor configuration. The numerical results showed that the propeller slipstream effect on the wings was highly dependent on the size of the propeller, and the major slipstream effect was working at 60% inboard wingspan, whereas less effects were observed towards the wingtip. The propeller slipstream increased the local angle of attack on the up-going blade side. This effect simultaneously augmented the section lift. The unsteady Reynolds-averaged Navier–Stokes (URANS) simulations helped to improve understanding of the interaction of the propeller wake and the wing-fuselage, which is an important aspect to guide the design of future efficient and controllable micro air vehicles. The results indicated that, in MAV designs, the slipstream from the propeller had a significant effect on the wing aerodynamics, regarding both performance and stability of the vehicle.

1. Introduction

In recent years, with the rapid development of micro-electro-mechanical system (MEMS) technology, the concept of miniaturization of fixed wing [1], rotary wing [2], and flapping wing [3,4] designs have been rapidly carried out. With the implementation and promotion of the civil-military integration strategy worldwide, the application of micro air vehicle (MAV) technology in the civil field has made great progress. Substantial reduction in manufacturing cost was achieved due to the gradual maturity of MAV research and development of technologies, and such small flying vehicles have been widely used in various fields including agricultural plant protection, power line patrol, police and fire protection and search and rescue, damage assessment, surveillance, and reconnaissance [5,6]. However, the miniaturization of MAV designs still faces multiple challenges, including low Reynolds number aerodynamics, miniaturized structure design, low-speed propeller propulsion system, flight stability and control, micro remote control, and sensor devices, etc. The key issue for a MAV toward miniaturization in general is: airfoil performance deteriorates rapidly as the chord Reynolds number decreases below 100,000 [7]. This is because of the much stronger viscous effects including separation, transition, and reattachments, influencing the lifting surface performance. The formation of the laminar separation bubbles (LSBs), one of the pronounced phenomena which occurs in low Reynolds number aerodynamics, can significantly alter the boundary layer behavior and stall characteristics [8]. Early research focused on the aerodynamic characteristics of isolated wings, with more representative work such as Torres and Mueller [9], Mizoguchi et al. [10], Chen et al. [8,11], and Traub et al. [12]. However, the effects of propeller slipstream on aerodynamics and flow physics have not yet been fully investigated and understood. Conventionally, it is common for the propulsion system to be considered separately from the wing and the airframe.

For MAVs, the interaction between the propeller and the aerodynamic surfaces can be much more significant due to the relative size of the propeller rotor disc to the MAV wing lifting surface, affecting its stability, performance, power consumption, noise, and overall endurance. Most of the MAVs are equipped with an electric motored propulsion system that contributes to the simplicity of operation. An obvious feature is that the aerodynamic efficiency of the wing control surface is significantly improved due to the flow separation that is prevented based on an extra momentum added to the flow from the propeller [13]. Most recently, propeller effects have been investigated on the vertical take-off and landing of MAVs [14], in which the authors found that the control surfaces became more effective as the slipstream strengthened. Therefore, the aerodynamic performance for the isolated wings is no longer applicable after a propeller is installed. The propeller swirling flow significantly modifies the surface pressure distribution and a considerable shifting of the center of pressure occurs, resulting in a change in overall pitching moment [15]. The swirling flow produced by the propeller also generates an additional yawing moment and a side slip can occur [16]. This is due to the propeller swirling flow effectively modifying the angle of attack of the downstream wing, thereby changing the wing’s circulation [17]. On the other hand, the local angle of attack of the propeller blades is also modified by the upwash of the wing behind the propeller. The downward rotating blade increases the local angle of attack and increases the local lift and blade loading, which augments the thrust and torque on the blade.

The pitching moment coefficient (C_m), as an important characteristic, is often overlooked when considering airfoil performance. The flying wing configuration is normally desired for MAV designs and such a flying wing configuration must make a negative slope of the pitching moment curve (C_m,α), and Cm at α = 0° (C_m,0) must be positive (+) (Nickel and Wohlfahrt, 1994) [18]. The requirements cannot be achieved by a traditional cambered airfoil, and hence a reflex camber must be designed which has a concave shape and locates near the trailing edge of the airfoil, shifting the pressure distribution aft and placing more upward force on the lower surface of the airfoil. However, the addition of reflex camber reduces the overall lift generation capabilities of the airfoil but is necessary for stable MAV flight. The limited amount of literature on cambered plate wings with a reflexed camber designed at trailing edge for low Reynolds numbers suggests that there is a need to expand the understanding of propeller-induced flow effects on reflex camber aerodynamic characteristics. To accomplish this goal, a computational study was performed to investigate the mutual aerodynamic influences between an MAV configuration and its tractored propeller. This is based on our previous studies on laminar separation bubbles [8] and planform effects [11] in a similar low Reynolds number range. The propeller–wing mutual interaction was studied and the comparison of wing aerodynamic characteristics of propeller on and off was conducted.

2. Computational Framework

2.1. Governing Equations and Solution Details

All simulations for this study solved Reynolds-averaged Navier–Stokes solutions. The flow equations were solved using the finite volume method, and the sliding mesh technique was applied for the interface between rotational and stationary domains. This led the computational domain to update the deformed mesh but change in time. Hence, the moving mesh integrated over a control volume was coupled with continuity and momentum equations, which can be written as follows:

(1)$\frac{\partial }{\partial \mathrm{t}}{{\displaystyle \int }}_{{\mathrm{V}}_{\mathrm{c}}}^{}{\mathsf{\rho }}^{\mathrm{f}}{\mathrm{dV}}_{\mathrm{c}}+{{\displaystyle \oint }}_{\mathrm{S}}^{}{\mathsf{\rho }}^{\mathrm{f}}\left(\upsilon -{\upsilon }_g\right)·n \mathrm{dS}=0$

(2)$\frac{\partial }{\partial \mathrm{t}}{{\displaystyle \int }}_{{\mathrm{V}}_{\mathrm{c}}}^{}{\mathsf{\rho }}^{\mathrm{f}}\upsilon {\mathrm{dV}}_{\mathrm{c}}+{{\displaystyle \oint }}_{\mathrm{S}}^{}{\mathsf{\rho }}^{\mathrm{f}}\upsilon \left(\upsilon -{\upsilon }_g\right)·n\mathrm{dS}=-{{\displaystyle \oint }}_{\mathrm{S}}^{}\tau ·n\mathrm{dS}-{{\displaystyle \oint }}_{\mathrm{S}}^{}\mathrm{p}·n\mathrm{dS}$

where, ${\mathsf{\rho }}^{\mathrm{f}}$ is the fluid density, $\upsilon$ is the fluid velocity vector which is based on the Cartesian coordinates, ${\upsilon }_g$ is the moving grid velocity, $\tau$ is the molecular momentum transport tensor, $\mathrm{p}$ indicates the pressure gradients term, ${\mathrm{V}}_{\mathrm{c}}$ is the mesh cell control volume, n is the normal vector to the control volume surface, and S is the control volume surface area.

The $\mathsf{\gamma }-{\mathrm{Re}}_{\mathsf{\theta }}$ transition model of Menter et al., 2006 [19], coupled the SST model with transport equations for the intermittency and ${\mathrm{Re}}_{\mathsf{\theta }}$ (i.e., momentum-thickness Reynolds number). The major improvement of this transition model was that it did not rely on nonlocal parameters, hence it was more suitable for modern CFD codes and complicated transitional-flow simulations. Further, a special modification to the intermittency was included to allow for separation-induced transition prediction. The SST $\mathsf{\gamma }-{\mathrm{Re}}_{\mathsf{\theta }}$ has been used widely at both low [20] and high Reynolds numbers in aerodynamic applications. This paper adopted the $\mathsf{\gamma }-{\mathrm{Re}}_{\mathsf{\theta }}$ transition model, which is proportional to the maximum strain-rate Reynolds number, i.e., presents the advantage of being a local property. The vorticity Reynolds number (${\mathrm{Re}}_{\mathrm{v}}$) is defined as:

(3)${\mathrm{Re}}_{\mathrm{v}}=\frac{{\mathsf{\rho }\mathrm{y}}^2}{\mathsf{\mu }}\frac{\partial \mathrm{u}}{\partial \mathrm{y}}=\frac{{\mathsf{\rho }\mathrm{y}}^2}{\mathsf{\mu }}\mathsf{\Omega }$

The $\mathsf{\Omega }$ and $\mathrm{y}$ in Equation (3), are the vorticity and the wall normal distance, respectively, and the maximum value of ${\mathrm{Re}}_{\mathrm{v}}$ is dependent on the ${\mathrm{Re}}_{\mathsf{\theta }}$. The momentum-thickness Reynolds number transport equation was used to capture the nonlocal effect of freestream-turbulence intensity and pressure gradient at the boundary-layer edge, which indicates where transition onset occurs, and is defined as:

(4)$\frac{\partial \left(\mathsf{\rho }\tilde{\mathrm{R}}{e}_{\mathsf{\theta }\mathrm{t}}\right)}{\partial \mathrm{t}}+\frac{\partial \left({\mathsf{\rho }\mathrm{U}}_{\mathrm{j}}\tilde{\mathrm{R}}{e}_{\mathsf{\theta }\mathrm{t}}\right)}{\partial {\mathrm{x}}_{\mathrm{j}}}={\mathrm{P}}_{\mathsf{\theta }\mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left[{\mathsf{\sigma }}_{\mathsf{\theta }\mathrm{t}}\left(\mathsf{\mu }+{\mathsf{\mu }}_{\mathrm{t}}\right)\frac{\partial \tilde{\mathrm{R}}{e}_{\mathsf{\theta }\mathrm{t}}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right]$

The transport equation for intermittency is used to trigger the transition process ($\mathsf{\gamma }>0$) and is defined as:

(5)$\frac{\partial \left(\mathsf{\rho }\mathsf{\gamma }\right)}{\partial \mathrm{t}}+\frac{\partial \left({\mathsf{\rho }\mathrm{U}}_{\mathrm{j}}\mathsf{\gamma }\right)}{\partial {\mathrm{x}}_{\mathrm{j}}}={\mathrm{P}}_{\mathsf{\gamma }}-{\mathrm{E}}_{\mathsf{\gamma }}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left[\left(\mathsf{\mu }+\frac{{\mathsf{\mu }}_{\mathrm{t}}}{{\mathsf{\sigma }}_{\mathrm{f}}}\right)\frac{\partial \mathsf{\gamma }}{\partial {\mathrm{x}}_{\mathrm{j}}}\right]$

However, when the boundary separates, the modification of intermittency to unity is defined as:

(6)${\mathsf{\gamma }}_{\mathrm{sep}}=\min \left\{8\max \left[\left(\frac{{\mathrm{Re}}_{\mathrm{v}}}{2.193{\mathrm{Re}}_{\mathsf{\theta }\mathrm{c}}}\right)-1,0\right]{\mathrm{e}}^{-{\left(\frac{{\mathrm{R}}_{\mathrm{T}}}{15}\right)}^4},5\right\}{\mathrm{F}}_{\mathsf{\theta }\mathrm{t}}$

In the above equations, (5), ${\mathrm{F}}_{\mathrm{length}}$, and ${\mathrm{Re}}_{\mathsf{\theta }\mathrm{c}}$ are two key functions, in which the former controls the transition extent and the later determines the onset of transition. The source term ${\mathrm{P}}_{\mathsf{\gamma }}$ is activated when the local strain-rate Reynolds number exceeds the local transition-onset criterion. The destruction source term ${\mathrm{E}}_{\mathsf{\gamma }}$ enables the relaminarization prediction when the transition-onset criterion is no longer satisfied, and vanishes in the fully turbulent regime. A complete description of the model is available in the article by Menter et al., 2006 [19]. This model has been used by numbers of researchers for low Re transitional flows. For example, a detailed study on two specific parameters (${\mathrm{F}}_{\mathrm{length}}{\mathrm{and}\mathrm{F}}_{\mathrm{onset}}$ are used in the intermittency equation for controlling the length of transition region and onset location of transition, respectively) was shown by Suluksna [21]. Benyahia [22] conducted a validation study for the model for low Re number flows. According to the comparison between the numerical and experimental data, it was shown that the $\mathsf{\gamma }-{\mathrm{Re}}_{\mathsf{\theta }}$ model accurately predicted the location and extent of the two-dimensional laminar separation bubble. Counsil [23] also studied the two-dimensional airfoils using the transition model, showing that the transition model was accurate in the preturbulent regions. A comparative study of four airfoils at low Reynolds numbers using the different transition modeling methods, namely, $\mathsf{\gamma }-{\mathrm{Re}}_{\mathsf{\theta }}$ and ${\mathrm{e}}^{\mathrm{N}}$, was performed by Seyfert [24].

A blended first- and second-order accurate scheme was implemented for the spatial discretization, which switched from the latter to the former in regions of steep spatial gradients, based on the boundedness principle of Barth and Jespersen [25]. SIMPLEC algorithm [26] with a staggered grid and the second order-accurate central difference method was enforced, and such pressure–velocity coupling used coordinated under-relaxations for the momentum and pressure corrections to improve the convergence which was inherently slow in the original SIMPLE method. An implicit second-order quadratic backward approximation with an iterative procedure was realized for temporal discretization. The nonlinear coefficients were updated within each inner loop while the outer loop advanced the solution in time.

The entire computational domain consisted of a rotational zone containing the propeller and a stationary zone containing the MAV wing and fuselage, as shown in Figure 1. The multiple zones were connected with each other through non-conformal mesh [27]. The rotating domain contained a cylinder-type boundary, and the central axis which was coincident with the rotational axis of the propeller. The height and radius of the rotational domain were determined based on the propeller diameter (D_ia), indicating 2.55D_ia and 7.6D_ia, respectively. The stationary domain, on the other hand, was set as a cubic block with a distance being roughly about 12.75D_ia, 20.5D_ia, and 12.75D_ia for upstream, downstream, and height, respectively; see Figure 1e,f. Structured mesh for both rotational and stational domain was generated, shown in Figure 1a,b. Mesh with high quality was considered and generated due to the importance of interpolation relationship between the interface surfaces, and a mapped mesh topology was proceeded between the interface boundaries. For the rotating domain, two O-topologies were created to cover the propeller and the center spinner segments. The O-grid included 10 cells normal to the propeller wall surface with a first cell distance of 2 × 10⁻⁵ m. There were 30 grid points in propeller radial direction and 56 grid points in the circumferential direction. A cylindrical wake block was used to ensure a good resolution of the blade wakes and the tip vortices, which have a significant influence on the MAV domain. For the outer domain, a similar mesh topology as we showed in the validation case section was used. The total size of the mesh was about 8 × 10⁶ nodes.

The standard characteristic boundary conditions were applied on the farfield boundaries. For the case of low Reynolds number flow, total pressure, incoming freestream velocity, freestream turbulent intensity ($T_i$), and the turbulence length scale were imposed at the inlet, whereas pressure was prescribed at the outlet boundary. Furthermore, at the solid wall, the non-slip boundary condition was applied. The turbulent kinetic energy was set to zero, and the pressure on the wall had zero normal gradients. High performance computing (HPC) with a parallel-processing implementation over 48 partitions and efficient message-passing interface between the partitions was adopted.

2.2. Specifications of Validation Cases

A three-dimensional Zimmerman wing planform was selected for the validation purposes, and the geometry was investigated by Torres and Mueller experimentally [9]. It has a zero camber and the aspect ratio (AR) is two, and the aerodynamic mean chord is 0.1725 m which gives a corresponding Reynolds number of 100,000. The model has a thickness-to-chord ratio of 1.96% and 5-to-1 elliptical leading- and trailing-edges. The aerodynamic center was assumed to be at the 25% point of the mean aerodynamic chord. To minimize the boundary condition effects, two different types of boundary settings were tested, and results were compared with the experiment data. Model 1 had an H-type mesh topology and the boundary conditions were set as velocity-inlet, pressure-outlet, and a non-slip wall boundary was applied on the wing. The farfield boundary condition was applied for the outer boundary. On the other hand, Model 2 was set as wing tunnel model boundaries. It has a velocity-inlet, pressure-outlet, and a wall boundary was applied for the wing and the outer domain faces, respectively (listed in Table 1). Model 2 was set as the wing tunnel model which had a contraction ratio of 20.6 to 1 and a rectangular inlet contraction cone. The freestream turbulence intensity in the test section was measured to be less than 0.05%. The test section was 182 cm long with a square cross-sectional area of 61 cm × 61 cm [9].

An O-grid mesh topology was applied around the wing planform to capture the detailed flow field near the boundary layer for Model 1 and Model 2 configurations. Therefore, the grid points were concentrated near the wing planform edge, as shown in Figure 2d. For the wind tunnel setting case, a symmetric boundary condition was applied at the wing root plane, as shown in Figure 2c. The no-slip wall boundary condition was enforced on the wing surface. To minimize the farfield boundary condition effects, Model 1 had the domain which is set at 25$\overline{c}$ upstream, 35$\overline{c}$ downstream, and the upper and lower boundaries were placed at 25$\overline{c}$ away from the airfoil leading-edge. However, Model 2 was set as wing tunnel model boundaries which had exactly same dimensions as the real wind tunnel.

The mesh sensitivity was studied for the three-dimensional cases for α = 4° and the results are shown in Table 1. The baseline mesh (G1) had about 1.5 million mesh elements and model with a fine mesh (G2) had a doubled size of 3 million. Close to the wall, there were about 60 grid points within the boundary layer, and in the turbulent region, the y⁺ value of the first cell distance was ensured to be in the order of O(1). The stretching ratio for the mesh was less than or equal to 1.2. Validation results are summarized in Table 1. The aerodynamic results comparisons, in Figure 2e,f, showed that Model 2 (wind tunnel settings) gave a better lift coefficient as compared with the experimental data. The results from Model 1 were close to the experimental data at low incidences but under-predicted at high incidences. The drag values showed a similar conclusion, showing potentially the stronger wall interference at higher incidences (Mueller [9]). From the mesh sensitivity study, the larger mesh size showed a reasonably better comparison than the coarse mesh. Therefore, for the propeller-wing-fuselage model presented later in this paper, a similar mesh topology, as shown for the validation case, was chosen and applied on the wing-fuselage part. The mesh for the propeller was integrated with the wing-fuselage mesh. The general topology is that an O-grid mesh was used around the propeller, and an H-type mesh was on the outer zone inside the rest of the rotating domain (details are shown in Figure 1).

2.3. Present Investigation Case

The model used in the current study was based on the flying wing MAV developed at the University of Sheffield [8,28]. This model is composed of a Zimmerman planform wing, a fuselage, a vertical stabilizer, and a propeller in a tractor configuration. The model is shown in Figure 3c with a coordinate system with x in the chordwise direction, y normal, and z spanwise. The wing has a mean aerodynamic chord of 0.221 m and an aspect ratio of 2.12. The freestream speed is 8.4 m/s, and the corresponding Reynolds number is 1.3 × 10⁵.

The airfoil used in this investigation was a simple cambered thin airfoil, as shown in Figure 3b. It included both positive and negative cambers and the relevant parameters are listed in Table 2. The positive camber was designed to have better aerodynamic performance whereas the reflex camber was designed to maintain stable level flight, resulting in decreased flight times. Not only the static stability but also the dynamic stabilities (for control handling) were important for the design. Similar to Torres and Mueller [9], the lift center location was calculated from the normal force and pitching moment taken from the $\overline{c}$/4 location of the wing. Although the location of lift center was close to $\overline{c}$/4 point at low angles of attack, it shifted towards the trailing edge as the incidence increased. The reason behind this is from both the increasing trailing edge separation and the strengthening wing tip vortices. In this investigation, the mean aerodynamic chord,$\overline{c}$, was used (Figure 3a).

The wing had a constant thickness of 2 mm and the shape of the cross section was rectangular with no sharp leading- and trailing-edge, and the airfoil is shown in Figure 3b. The camber was defined as maximum convex camber (h₁/$\overline{c}$), maximum concave and reflex camber (h₂/$\overline{c}$), maximum concave camber location (d₁/$\overline{c}$), and the maximum reflex camber location (d₁/$\overline{c}$); further details are shown in Ref. [11].

Considering the MAV configuration, the fuselage was a substantial component to accommodate the payload and propulsion devices, such as battery, motor, and servos, etc. Figure 3d shows the fuselage dimension for our model. It had a front height, h_f, of 0.059 m, rear height, h_t, of 0.021 m, and a total length h_L = 0.216 m (Table 2). This design of the fuselage was dictated by the size and placement of the components. The fuselage layout affects the center of gravity margin and hence the static stability. For this purpose, the battery was designed to be movable to adjust the center of gravity. Another interesting point is how the fuselage affects the overall aerodynamics. Brion [29] simulated the fuselage and wing separately and the relevant aerodynamic forces are shown and discussed. However, the authors did not mention anything about the interaction between the wing and the fuselage. Ramamurti’s [30] numerical results showed MAV with fuselage reduced the lift-to-drag ratio dramatically and the drag for all configurations considered was nearly the same. The effects of the fuselage were also investigated in the present study.

The slow-flyer propeller was chosen for the current MAV study. The propeller was installed at a distance d_t = 0.068 m ahead of the wing planform. It had a diameter of 8 in, the pitch was 4 in, and the hub diameter was 0.014 m. Figure 3e shows the propeller geometry and the blade azimuth angle. The propeller rotated in an anti-clockwise direction viewing from the front, and Figure 3f shows the propeller geometric characteristics. The aerodynamic balance determined a horizontal component X and a vertical component Y of the total force acted on the MAV model (Figure 3d). To obtain the overall lift and drag force on the model, horizontal and vertical components (i.e., X and Y components) were transferred into L and D components (i.e., based on the incoming freestream coordinates), as shown in Equation (7).

(7)$$\begin{gathered} L=L_T+L_W,andD=\left(-T_T\right)+D_W \\ {L}_T=Y_{T}cos\alpha +X_{T}sin\alpha ,and{L}_{WF}=Y_{WF}cos\alpha -X_{WF}sin\alpha \\ {T}_T=-Y_{T}sin\alpha +X_{T}cos\alpha ,and{D}_{WF}=Y_{WF}sin\alpha +X_{WF}cos\alpha \end{gathered}$$

A time step of 1.5724 × 10⁻⁵ s with 30 sub-iterations was applied for this study (is equivalent to 0.5 degrees per time step). To have reasonable numerical results, the y⁺ value of the first grid point in the order of O(1) was required. Figure 4 shows the numerical aerodynamic force coefficients versus the blade azimuth angle. It shows that periodic pulses were produced. This type of signature was found to be relatively independent of the advance ratio and appeared to be mainly associated with the local loading on the propeller itself.

A steady state flight condition is defined as one for which all motion variables remain constant with time relative to the body-fixed axis system XYZ. Mathematically speaking, steady state flight conditions are ${\overrightarrow{\dot{V}}}_p=0$, and $\overrightarrow{\dot{\omega }}=0$, and implies that MAV does not have any acceleration in any direction (i.e.,$\frac{du}{dt}=\frac{dv}{dt}=\frac{dw}{dt}=0$), and that the roll, yaw, and pitch rates are zero. Therefore, a propeller with rotational speed of 555 rad/s satisfied the steady state flight condition, as shown in Figure 5. Namely, the lift equaled to the weight and drag equaled to the propeller thrust.

3. Results and Discussion

3.1. Propeller Slipstream Effects on Aerodynamic Performance

From previous investigations on both pusher and tractor propeller configurations, the latter was generally found to have higher aerodynamic efficiency and produce better control surface’s effectiveness at landing and take-off. The MAV model considered here had a tractor configuration, and this was mainly to ensure the safety of the hand-throwing technique used on MAV take-off. Figure 6 illustrates the tractor propeller configuration effect on the left and right side of the Zimmerman wing. The propeller wake flow generated areas of upwash and downwash on the wing. Figure 6a shows the wing-airfoil aerodynamic loading in the upwash case where the propeller tangential velocity was greater than the downwash. The resultant force was tilted forward, and a localized wing thrust was created at this section. Figure 6b shows the aerodynamics of wing section under the propeller downwash. The wing downwash and the propeller tangential velocity acted together, resulting in a further reduction in the local angle of attack. The resultant force, F, shifted further backwards and resulted in an increased sectional drag.

Figure 7a–f shows the aerodynamic coefficient of the propeller at a rotational speed of 555 rad/s. The forces and moments for different components clearly showed in the plots individually (propeller, wing-fuselage, and the vertical stabilizer). The lift coefficient in Figure 7a showed a linear variation for α < 6° and became nonlinear as the angle of attack was further increased. The maximum lift coefficient, C_L,max, was about 0.81 at the incidence of 12°. The propeller had negligible lift contribution at low angles of attack. The vertical stabilizer also showed a near zero contribution on the lift coefficient. At high angles of attack, the propeller had a significant contribution on the lift, which improved the maximum lift coefficient.

The drag coefficient, in Figure 7b, increased continuously with the angle of attack when the rotational speed was fixed at 555 rad/s. A nearly constant thrust was produced by the propeller with various angles of attack (C_T is roughly around −0.09). The main drag force was formed from the wing-fuselage. The vertical stabilizer, however, showed a negligible drag contribution for the entire range of incidences.

Figure 7c shows the side force generated by the MAV model. In general, the MAV side force increased continuously as the angle of attack was increased. At 0°, the main side force coefficient was generated by the wing-fuselage, approximately 0.0037 (or 0.0143N). However, the propeller generated a negative side force coefficient of −0.0016 (or −0.0062N), and the side force from the vertical stabilizer was nearly zero at a 0° angle of attack. In general, the side force from the MAV had a small effect on the MAV as compared with the total weight.

The rolling, yawing, and pitching moments are shown in Figure 7d–f, respectively. A continuous increase of the rolling moment coefficient is shown in Figure 7d. At 0°, the overall rolling moment coefficient was 0.0014 (or 0.0024 Nm). A peak value was found at 18° mainly from the wing-fuselage part. The yawing moment coefficient, C_My, indicated that the MAV with a propeller spun clockwise as viewed from the rear and the moments caused a left yaw. However, it was negligible at low incidences around α = 0° with a value of approximately 9e⁻⁵. A linear increase in the propeller yawing moment can be found from the plot.

The pitching moment for analyzing the static longitudinal stability had a natural statically longitudinal stable region at α < 0° and positive contribution on the longitudinal stability side as the incidence increased. The pitching moment coefficient slope was closer to linear than that for the isolated wing-fuselage model in Ref. [11]. In contrast, the pitching moment slope, C_Mα, was about −0.0033 (for the linear section at incidence between 6° and 36°). The propeller showed a negative contribution on the statically longitudinal stability and a positive pitching moment slope can be identified in Figure 7f.

For the cruise condition, the lift force at zero degree of angle of attack was equal to the MAV weight and the thrust force could overcome the drag force, shown in Figure 7b. Figure 8a shows the locations of C_p distribution on the Zimmerman wing at r/R = 1 (or z/b = ± 0.25). The C_p distribution locations are indicated by the dash line for the up-going and down-going blade sides. Both C_p with and without slipstream at 0° flight condition are presented in Figure 8b–e: the instantaneous C_p distribution is plotted with various azimuth angles (in solid line). The C_p plots show that the positive lift force was generated at the positive camber area around x/c = 0.25 and the negative lift force was obtained from the negative camber, which was located near the trailing edge. The negative lift was also generated at near leading-edge area on the down-going blade side. The averaged C_p at 2z/b = 0.5, in Figure 8f, shows that a similar amount of lift was produced by the isolated MAV model and the MAV propeller model. However, the negative lift generated at the reflex camber region by the MAV propeller model was quite significant, and more negative lift was produced due to the propeller slipstream effect. In comparison, the amount of negative lift almost doubled as compared with the isolated MAVs (dashed line). Therefore, less reflex camber could be used for the MAV with the propeller installed. The turning point is the point where the negative lift force started to occur, which shifted toward the leading edge slightly due to the propeller slipstream effect.

Figure 9 shows the two-dimensional lift and drag coefficient distribution at various spanwise locations. The difference of the lift coefficients distribution for both wing sides is clearly shown. The drag coefficient, on the other hand, showed a higher value on the down-going blade and a slightly lowered value on the up-going blade. The possible explanation could be that in the up-going blade region, the propeller swirl counteracted the effects of the wing downwash such that the local angles of attack were increased. This effect simultaneously augmented the section lift and rotated the force vector forwards, which reduced the drag component at the section, as detailed in Figure 6. The sudden increase in lift or decrease in drag coefficient indicated the aerodynamic interaction on the wing and the fuselage from the propeller wake. In general, the major slipstream effect was working at 60% inboard wingspan,$2z/b\le 0.6$, and less effects were observed toward the wingtip.

In general, the propeller slipstream had a significant effect on both positive and negative lift. Identical negative lift was mainly generated at reflex camber area for both downward- and upward-going blade sides and asymmetrical forces were produced at the positive camber region due to stronger propeller slipstream, as shown clearly from the pressure distribution in Figure 8.

3.2. Propeller Slipstream Effect on the Flow Structure

Flow structures around and downstream of the propeller are presented and discussed in this section. Figure 10 shows the asymmetric flow structure of the MAV propeller model for the propeller at two different azimuth angles, $\phi =0^{°} \mathrm{and} {90}^{°}$, respectively, with $\mathsf{\alpha }=0^{°}$. With the propeller slipstream, the wingtip vortex was reduced as compared with the MAV without the propeller slipstream, shown in the previous paper [11]. A clear asymmetric flow appeared, which followed the airfoil local inflow angle of attack evolution. Larger local angle of attack was formed on the upward-going blade side than that on the downward-going blade side, as shown in Figure 10e,g, respectively. Toward the wingtip, the local angles of attack on both sides were about the same. Details are shown in Figure 10f,h.

At $\phi =0^{°}$, the leading edge separation bubble formed on the upper wing surface and it was located near the center plane region, as shown in Figure 10i. The low pressure region, as marked with dark blue color, can be clearly seen on the upper surface, indicating the asymmetric flow due to the propeller slipstream effects. On the lower wing surface, two separated and asymmetric bubbles formed underneath the wing near the leading edge. The smaller bubble formed on the upward-going blade side and the larger one on the downward-going blade side. The shape of the bubble at different spanwise locations can be found in Figure 10e–h. The areas with relative larger pressure force were on the lower surface at the up-going blade side (Figure 10d). This was due to the propeller tangential flow impinging on the lower wing surface. However, on the down-going blade side, the tangential flow impinged on the upper surface instead.

4. Conclusions

The effect of a tractor propeller on the flow over a Zimmerman wing-fuselage model was investigated. The MAV propeller configuration was investigated with various advance ratios, and the equivalent rotational speed of 555 rad/s for the steady state condition was designed and the results were analyzed. Strong unsteady fluctuations were found in the slipstream region, showing development of periodic aerodynamic forces and moments. The sinusoidal fluctuation resulted in a large impact on the wing-fuselage aerodynamics. A severe mutual influence on the flow structure between the propeller and the MAV was observed as the angle of attack varied.

With the presence of the propeller, periodic aerodynamic force and moment on wing-fuselage-propeller configuration was observed. The contribution on the lift from the propeller increased as the angle of attack increased. While the propeller thrust and the drag force from the vertical stabilizer remained nearly constant for the incidence range, the drag from the wing-fuselage increased dramatically as the angle of attack increased.

The propeller slipstream effects on two-dimensional lift and drag coefficient distribution at various spanwise locations were clearly shown. The drag coefficient showed a higher value on the down-going blade and a slightly lowered value on the up-going blade. The possible reason could be that in the up-going blade region, the propeller swirl counteracted the effects of the wing downwash such that the local angles of attack were increased. This effect simultaneously augmented the section lift and rotated the force vector forwards, which reduced the drag component at the section. The sudden increase in lift or decrease in drag coefficient indicated the aerodynamic interaction on the wing and the fuselage from the propeller wake. In general, the major slipstream effect was working at 60% inboard wingspan, and less propeller slipstream effects were observed toward the wingtip.

The propeller slipstream showed a significant effect on both positive (for positive camber) and negative (for reflex camber) lift. The pressure distribution comparison between the isolated wing planform and propeller-installed model indicated that the down force caused by reflex camber due to the propeller slipstream effect was significant. In other words, the reflex camber of the wing needs to be reduced to achieve better aerodynamic performances for wing-propeller configuration.

Author Contributions

Conceptualization, methodology, and Writing—Original draft preparation, Z.C.; validation, F.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was sponsored by the National Natural Science Foundation of China under Grant No. 11672132, 11672133, 12002161, and 12032011.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors wish to thank Gabriel Torres and Thomas J Mueller of University of Notre Dame for providing the wind tunnel data.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Non-conformal mesh interface for propeller MAV geometry, (a) propeller-MAV configuration, (b) detailed structured mesh for propeller-MAV, (c) Mesh slice cut though domains, (d) propeller structured mesh, (e) computational domain in front view, and (f) computational domain in side view.

Figure 2. Mesh topology and aerodynamic coefficients for wing planform, (a) C-H mesh topology, (b) mesh at leading-edge, (c) O-grid mesh topology, (d) spanwise mesh spacing, (e) lift coefficient, (f) drag coefficient.

Figure 3. MAV and propeller geometrical descriptions, (a) Mav in top view, (b) airfoil profile, (c) MAV in 3D view, (d) aerodynamic force decomposition, (e) propeller azimuth angle, and (f) propeller geometrical information(chord and pitch angle distribution).

Figure 4. Propeller time history aerodynamic forces coefficient development, (a) drag coefficient, (b) lift coefficient, (c) side force coefficient, (d) rolling moment coefficient, (e) yawing moment coefficient, and (f) pitching moment coefficient.

Figure 5. Drag and tractor coefficient with a function of advance ratio.

Figure 6. Propeller wing in tractor configuration, (a,b) wing local effective angle of attack due to propeller slipstream effect.

Figure 7. Averaged aerodynamic coefficients versus α at J = 0.468, (a) lift coefficient, (b) drag coefficient, (c) side force coefficient, (d) rolling moment coefficient, (e) pitching moment coefficient, and (f) yawing moment coefficient.

Figure 8. Instantaneous and averaged Cp distributions at various spanwise location with α = 0° and ω = 555 rad/s, (a) MAV wing spanwise locations, (b,c) wing spanwise location of 2z/b = 0.25 and propeller zimuth angle of 0° and 90°, (d,e) wing spanwise location of 2z/b = 0.5 and propeller zimuth angle of 0° and 90°, (f) pressure coefficient comparison between isolated wing and wing-MAV at spanwise location of [Forumla omitted. See PDF.].

Figure 9. Instantaneous Cl, Cd distribution at various propeller azimuth angles at ω = 555 rad/s, (a–d) lift and drag coefficient distribution along the wing span with changing of propeller azimuth angle from [Forumla omitted. See PDF.].

Figure 10. Instantaneous flow structure around the MAV at [Forumla omitted. See PDF.], (a–d) flow field around MAV, (e,f) airfoil sectional flow structure with up-going blade at spanwise location of 2z/b = 0.25 and 0.5, (g,h) airfoil sectional flow structure with down-going blade at spanwise location of 2z/b = 0.25 and 0.5, (i) flow structure at propeller’s azimuth angle [Forumla omitted. See PDF.].

Table 1

Grid-sensitivity analysis for α = 4°.

	Case	Grid Size	CL	CD
Model 1	G1	190 × 135 × 60	0.1744	0.01913
	G2	220 × 165 × 90	0.1775	0.02068
	G3	250 × 195 × 120	0.1781	0.02052
Model 2	G1	190 × 135 × 60	0.1899	0.01890
	G2	220 × 165 × 90	0.1900	0.01877
	G3	250 × 195 × 120	0.1901	0.01879
	Experiment [1]	NA	0.1906 ± 0.02	0.0220 ± 0.003

Table 2

Propeller location and fuselage specifications.

MAV model	${\mathbf{d}}_{\mathbf{t}}\mathbf{(}{\mathbf{d}}_{\mathbf{t}}/\overline{\mathit{c}}$)	${\mathbf{h}}_{\mathbf{f}}\mathbf{(}{\mathbf{h}}_{\mathbf{f}}/\overline{\mathit{c}}$)	${\mathbf{h}}_{\mathbf{L}}\mathbf{(}{\mathbf{h}}_{\mathbf{L}}/\overline{\mathit{c}}$)	${\mathbf{h}}_{\mathbf{t}}\mathbf{(}{\mathbf{h}}_{\mathbf{t}}/\overline{\mathit{c}}$)	$\mathbf{t}(\mathbf{t}/\overline{\mathit{c}}$)
	0.068(31.7%)	0.059(26.6%)	0.216(97.6%)	0.021(9.49%)	2e−3(0.93%)

References

1. Shi, L.; Uy, C.; Huang, S.K.; Yang, Z.; Huang, G.P.; Wen, C.-Y. Experimental investigation on flow characteristics of a four-wing micro air vehicle. Int. J. Micro Air Veh.; 2016; 8, pp. 181-193. [DOI: https://dx.doi.org/10.1177/1756829316667770]

2. Désert, T.; Moschetta, J.; Bézard, H. Numerical and experimental investigation of an airfoil design for a Martian micro rotorcraft. Int. J. Micro Air Veh.; 2018; 10, pp. 262-272. [DOI: https://dx.doi.org/10.1177/1756829318794171]

3. De Wagter, C.; Karásek, M.; de Croon, G. Quad-thopter: Tailless flapping wing robot with four pairs of wings. Int. J. Micro Air Veh.; 2018; 10, pp. 244-253. [DOI: https://dx.doi.org/10.1177/1756829318794972]

4. Lankford, J.; Mayo, D.; Chopra, I. Computational investigation of insect-based flapping wings for micro air vehicle applications. Int. J. Micro Air Veh.; 2016; 8, pp. 64-78. [DOI: https://dx.doi.org/10.1177/1756829316646640]

5. Hassanalian, M.; Abdelkefi, A. Classifications, applications, and design challenges of drones: A review. Prog. Aerosp. Sci.; 2017; 91, pp. 99-131. [DOI: https://dx.doi.org/10.1016/j.paerosci.2017.04.003]

6. Hassanalian, M.; Abdelkefi, A. Design, manufacturing, and flight testing of a fixed wing micro air vehicle with Zimmerman planform. Meccanica; 2017; 52, pp. 1265-1282. [DOI: https://dx.doi.org/10.1007/s11012-016-0475-2]

7. Mueller, T.J.; DeLaurier, J.D. Aerodynamics of Small Vehicles. Annu. Rev. Fluid Mech.; 2003; 35, pp. 89-111. [DOI: https://dx.doi.org/10.1146/annurev.fluid.35.101101.161102]

8. Chen, Z.J.; Qin, N.; Nowakowski, A.F. Three-Dimensional Laminar-Separation Bubble on a Cambered Thin Wing at Low Reynolds Numbers. J. Aircr.; 2013; 50, pp. 152-163. [DOI: https://dx.doi.org/10.2514/1.C031829]

9. Torres, G.E.; Mueller, T.J. Low-aspect-ratio aerodynamics at low reynolds numbers. AIAA J.; 2004; 42, pp. 865-873. [DOI: https://dx.doi.org/10.2514/1.439]

10. Mizoguchi, M.; Itoh, H. Effect of Aspect Ratio on Aerodynamic Characteristics at Low Reynolds Numbers. AIAA J.; 2013; 51, pp. 1631-1639. [DOI: https://dx.doi.org/10.2514/1.J051915]

11. Chen, Z.J.; Qin, N. Planform Effects for Low-Reynolds-Number Thin Wings with Positive and Reflex Cambers. J. Aircr.; 2013; 50, pp. 952-964. [DOI: https://dx.doi.org/10.2514/1.C032102]

12. Traub, L.W.; Botero, E.; Waghela, R.; Callahan, R.; Watson, A. Effect of Taper Ratio at Low Reynolds Number. J. Aircr.; 2015; 52, pp. 734-747. [DOI: https://dx.doi.org/10.2514/1.C032559]

13. Ikami, T.; Kanou, K.; Takahashi, K.; Fujita, K.; Nagai, H. Flow Field on Wing Surface with Control Surface in Propeller Slipstream at Low Reynolds Number. Proceedings of the AIAA Scitech 2019 Forum; San Diego, CA, USA, 7–11 January 2019; American Institute of Aeronautics and Astronautics: Washington, DC, USA, 2019.

14. Lyu, X.; Gu, H.; Zhou, J.; Li, Z.; Shen, S.; Zhang, F. Simulation and flight experiments of a quadrotor tail-sitter vertical take-off and landing unmanned aerial vehicle with wide flight envelope. Int. J. Micro Air Veh.; 2018; 10, pp. 303-317. [DOI: https://dx.doi.org/10.1177/1756829318813633]

15. Marretta, R.M.A. Different Wings Flow Fields Interaction on the Wing-Propeller Coupling. J. Aircr.; 1997; 34, pp. 740-747. [DOI: https://dx.doi.org/10.2514/2.2252]

16. Phillips, W.F.; Niewoehner, R.J. Effect of Propeller Torque on Minimum-Control Airspeed. J. Aircr.; 2006; 43, pp. 1393-1398. [DOI: https://dx.doi.org/10.2514/1.20303]

17. Witkowski, D.P.; Lee, A.K.H.; Sullivan, J.P. Aerodynamic interaction between propellers and wings. J. Aircr.; 1989; 26, pp. 829-836. [DOI: https://dx.doi.org/10.2514/3.45848]

18. Nickel, K.; Wohlfahrt, M. Tailless Aircraft in Theory and Practice; American Institute of Aeronautics and Astronautics: Washington, DC, USA, 1994.

19. Menter, F.R.; Langtry, R.B.; Likki, S.R.; Suzen, Y.B.; Huang, P.G.; Volker, S. A Correction-Based Transition Model Using Local Variables- Part I: Model Formulation. J. Turbomach.; 2006; 128, pp. 413-422. [DOI: https://dx.doi.org/10.1115/1.2184352]

20. Colonia, S.; Leble, V.; Steijl, R.; Barakos, G. Assessment and Calibration of the γ-Equation Transition Model at Low Mach. AIAA J.; 2017; 55, pp. 1126-1139. [DOI: https://dx.doi.org/10.2514/1.J055403]

21. Suluksna, K.; Dechaumphai, P.; Juntasaro, E. Correlations for modeling transitional boundary layers under influences of freestream turbulence and pressure gradient. Int. J. Heat Fluid Flow; 2009; 30, pp. 66-75. [DOI: https://dx.doi.org/10.1016/j.ijheatfluidflow.2008.09.004]

22. Benyahia, A.; Houdeville, R. Transition prediction in transonic turbine configurations using a correlation-based transport equation model. Int. J. Eng. Syst. Model. Simul.; 2011; 3, 36. [DOI: https://dx.doi.org/10.1504/IJESMS.2011.038748]

23. Counsil, J.N.N.; Boulama, K.G. Low-Reynolds-Number Aerodynamic Performances of the NACA 0012 and Selig–Donovan 7003 Airfoils. J. Aircr.; 2013; 50, pp. 204-216. [DOI: https://dx.doi.org/10.2514/1.C031856]

24. Seyfert, C.; Krumbein, A. Evaluation of a Correlation-Based Transition Model and Comparison with the eN Method. J. Aircr.; 2012; 49, pp. 1765-1773. [DOI: https://dx.doi.org/10.2514/1.C031448]

25. Barth, T.; Jespersen, D. The design and application of upwind schemes on unstructured meshes. Proceedings of the 27th Aerospace Sciences Meeting; Reno, NV, USA, 9–12 January 1989; American Institute of Aeronautics and Astronautics: Washington, DC, USA, 1989.

26. Van Doormaal, J.P.; Raithby, G.D. Enhancedments of the SIMPLE Method for Predicting Incompressible Fluid Flows. Numer. Heat Transf.; 1984; 7, pp. 147-163. [DOI: https://dx.doi.org/10.1080/01495728408961817]

27. Domino, S.P. Design-order, non-conformal low-Mach fluid algorithms using a hybrid CVFEM/DG approach. J. Comput. Phys.; 2018; 359, pp. 331-351. [DOI: https://dx.doi.org/10.1016/j.jcp.2018.01.007]

28. Boller, C.; Kuo, C.M.; Qin, N. Biologically Inspired Shape Changing Aerodynamic Profiles and their Effect on Flight Performance of Future Aircraft. Adv. Sci. Technol.; 2008; 56, pp. 534-544.

29. Brion, V.; Aki, M.; Shkarayev, S. Numerical simulation of low Reynolds number flows around micro air vehicles and comparison against wind tunnel data. Proceedings of the 24th Applied Aerodynamics Conference; San Francisco, CA, USA, 5–8 June 2006.

30. Ramamurti, R.; Sandberg, W.; Lohner, R. Simulation of the dynamic of micro air vehicles. Proceedings of the 38th Aerodynamic Sciences Meeting & Exhibit; Reno, NV, USA, 10–13 January 2000.