Experimental Investigation of Propeller Induced Flow on Flying Wing Micro Aerial Vehicle for Improved 6DOF Modeling

In this research effect of propeller induced flow on aerodynamic characteristics of low aspect ratio flying wing micro aerial vehicle has been investigated experimentally in subsonic wind tunnel. Left turning tendencies of right-handed propellers have been discussed in literature, but not much work has been done to quantify them. In this research, we have quantified these tendencies as a change in aerodynamic coefficient with a change in advance ratio at a longitudinal trim angle of attack using subsonic wind tunnel. For experimental testing, three fixed pitch propeller diameters (5 inch, 6 inch and 7 inch), three propeller rotational speeds (7800, 10800 and 12300 RPMs) and three wind tunnel speeds (10, 15 and 20 m/s) have been considered to form up 27 advance ratios. Additionally, wind tunnel tests of 9 wind mill cases were conducted and considered as baseline. Experimental uncertainty assessment for measurement of forces and moments was carried out before conduct of wind tunnel tests. Large variation in lift, drag, yawing moment and rolling moment was captured at low advance ratios, which indicated their significance at high propeller rotational speeds and large propeller diameters. Side force and pitching moment did not reflect any significant change. <inline-formula> <tex-math notation="LaTeX">$\frac {L}{D}$ </tex-math></inline-formula> at trim point was found a nonlinear function of propeller diameter to wingspan ratio <inline-formula> <tex-math notation="LaTeX">$\frac {D}{b}$ </tex-math></inline-formula>, and propeller rotational speed. Rate and control derivatives were obtained using unsteady vortex lattice method with propeller induced flow effect modeled by Helical Vortex Modeling approach. In this research, we have proposed improved 6-DOF equations of motion, with a contribution of advance ratio <inline-formula> <tex-math notation="LaTeX">$J$ </tex-math></inline-formula>. It is concluded, that propeller induced flow effects have a significant contribution in flight dynamic modeling for vehicles with large propeller diameter to wingspan ratio, <inline-formula> <tex-math notation="LaTeX">$\frac {D}{b}$ </tex-math></inline-formula> of 22% or more.


I. INTRODUCTION
Powered flight because of Wright brothers' efforts has changed the world since early 20th century. The secret of Wright brothers' success lies in their appreciation, that what can generate lift for wings, can also generate thrust for propellers. The propeller blade is a screw that generates thrust by moving a large mass of air with less velocity [1]. Propeller efficiency is the measure of its effectiveness of converting engine power into propulsive thrust. The number of blades, propeller rotational speed, radius, configuration (trac-The associate editor coordinating the review of this manuscript and approving it for publication was Wen-Sheng Zhao . tor, pusher or both), pitch, solidity and twist, all affect efficiency of the propeller. By design, the propeller has a thicker airfoil at high angle of incidence near axis of rotation, while the incidence angle keeps on decreasing when approaching tip, refer figure 1. Additionally the tip has thinner airfoil as compared to airfoil near the hub [2]. This is to ensure even distribution of lift all along the length of a propeller blade. Early researchers appreciated, that the basics of propeller lie in the physics of a simple wing in a two-dimensional plane, which is essentially an airfoil. Like lift, drag and pitching moment are the main characteristics of airfoil; thrust, torque and power are considered to be the main characteristics for propellers [3]. Propellers' thrust propel the vehicle forward but, rotating propellers also generate unwanted moments. These unwanted moments need correction by the pilot to keep the airplane flying in desired direction [4]. In discussing forces produced by propellers the concept of actuator disc is more common which can be thought of as a converging tunnel, where an abrupt pressure jump exists and flow accelerates. This concept was proposed in classical momentum theory separately by Rankine [6] and Froude [7]. Momentum theory concluded, that fluid's axial velocity passing through the propeller disc is higher than the vehicle's speed. The accelerating velocity through this actuator disc is called an induced velocity denoted by ν. Propeller thus induces additional velocity in the flow by infusing more energy due to its rotation. Induce velocity far downstream of actuator disc is double the induce velocity just downstream of disc. The theory proposed propeller induce velocity as ν which can be calculated using equation 1. Momentum theory, however, ignored various important effects of viscosity and compressibility. Further, it lacked geometrical details of the blade.
Blade Element Theory, also known as strip theory, proposed by Drzewiecki addressed some discrepancies of momentum theory by calculating thrust, power required and torque on propeller blade cross section. The theory integrated forces on individual blade elements which were considered as airfoils [8]. It, however, ignored compressibility effects, tip loss effect, interference effect and downwash effect. In figure 2, thrust is parallel to whereas torque is perpendicular to the forward velocity of airplane. Lift is perpendicular to while drag is parallel to the vectorial sum of V ∞ , ν and 2rω.
Analogous to wing elliptical lift distribution, Betz Vortex theory considered propeller blades as lifting surfaces with circulation associated with bound vortex and vortex sheet that were continuously shed from trailing edge [9]. The shed vortices formed a rigid helicoidal vortex sheet that moves backward behind the propeller. At the same time, Prandtl proposed the tip correction factor for an approximate solution to the flow around helicoidal vortex sheets [10]. Tip correction factor presented acceptable results for low advance ratios and large number of blades. Goldstein solved the ideal distribution of circulation for helicoidal vortex system for small advance ratios by using Goldstein functions [11]. Larrabee [12] used previous works of Theodorsen [13] regarding the treatment of propellers with load distribution and presented design of the propeller blade for MIL (minimum induce loss). The theory was capable to calculate the induced velocities due to helical vortex sheets which constitute the slipstream of a propeller. Strong flow gradients both in the stream-wise and radial direction are produced due to self-induced velocities of the vortex system in the propeller slipstream. Important flow quantities to characterize slipstream are swirl velocity, axial velocity, pressure distribution, vorticity, helicity and contraction, refer figure 3. These characteristics of helical system induce unwanted forces and moments on a vehicle which needs to be quantified specially for vehicles with large propeller diameter to wingspan ratio, D b .

II. EFFECTS OF PROPELLER HELICAL SYSTEM ON AIRCRAFT AERODYNAMIC MOMENTS
The helicoidal vortex system because of peculiar characteristics of vorticity and helicity generates undesirable forces and moments especially at high propeller rotations (lower advance ratios ). Advance ratio is defined as the ratio of forward velocity to propeller rotational speed and propeller diameter, as depicted in equation 2.
The effect of propeller rotation is an increased drag (scrubbing drag), which tends to decrease propeller efficiency and induces reactionary moments, for instance unwanted yaw and unwanted roll [14]. These moments are generally known as left-turning tendencies on right-handed propellers. Such left turning effects are more pronounced in MAVs as majority of the wing span in covered by helical nature of airflow coming from the propeller [15]- [19]. These turning effects are slipstream effect, torque reaction, gyroscopic precession and P factor [5]. The propeller is rotated in a clockwise direction once viewed from the rear, therefore equal and opposite reaction gives an aircraft a left rolling tendency. This turning effect is known as torque reaction, refer figure 4.
The second propeller turning effect is a gyroscopic precession effect. It is known from rotating mass as precession effect. If a force is applied to the circumference of rotating mass, the reaction will act at 90 • in the direction of rotation. The propeller is a rotating mass so it too has a property of precession. If an aircraft pitches down, the force will act at the top of the propeller and due to precession, reaction will act 90 • to the point of force application. This will generate left yawing moment and aircraft nose will turn left. However, when an aircraft pitches up, force will act at the bottom of propeller and due to precession, reaction force will act at 90 • in the direction of rotation, which will generate right yawing moment. Likewise, whenever a propeller aircraft is yawed to left, force is applied to the right of propeller and due to precession, its reaction will act at the bottom of propeller in the direction of rotation and aircraft pitches up. Therefore yaw to left is coupled with pitching up tendency and yaw FIGURE 2. Velocity vector diagram on a cross section of propeller blade as proposed by blade element theory [8]. Aerodynamic forces on the propeller blade are also shown. What lift is to a fixed wing, is thrust to the propeller and what drag is to a fixed wing, is torque to the propeller. ω is the propeller rotational speed in cycles per second. V ∞ is the forward velocity of the airplane and V R is the resultant velocity vector of forward velocity and tangential velocity, V t . Propeller induced flow velocity component in the thrust direction is ν i cosφ. β is the geometric pitch of the propeller blade and α is the angle of attack which propeller blade encounters. φ is called the helix angle while θ is the induced angle. Figure adapted from [3].
to right is coupled with pitching down tendency in right handed propeller aircraft. In totality, gyroscopic precession effect of pitching the aircraft up, is a right yawing moment and gyroscopic effect of yawing the aircraft to right, is a nose-down pitching moment. Rolling the aircraft will not give any gyroscopic effect because the propeller disc is not being inclined. Schematic of gyroscopic precession effect is shown in figure 5. The third propeller turning effect is 'P' factor, refer figure 6. At low indicated speeds, the propeller disc is angled upwards with the direction of flight to make flow strike propeller disc at some positive angle of attack. With upwards inclination of the propeller disc specially during take-off, the down going blade moves slightly forward. Hence, the down going blade experiences an increase angle   of attack and therefore generates more thrust. The up going blade generates less thrust. If aircraft is viewed from above, this asymmetric blade effect can be appreciated as shown in figure 6. More thrust from down going blade and less thrust from up going blade will be generated at positive angles of attack. This asymmetric thrust loading on the propeller disc generates left yawing moment which is more pronounced at higher propeller rotational speed. FIGURE 6. Schematic of P factor effect. This figure shows that for right-handed propellers, more thrust is created by down moving blade as compared to up going blade when propeller disc is inclined at positive angle of attack. This asymmetry of force generates left yawing moment due to differential in forces on propeller disc. Figure adapted from [5].
Apart from these turning tendencies, there are more unwanted moments that are generated in right-handed propellers. The propeller thrust line can cause pitching moment once it is not aligned with Fuselage Reference Line (FRL). A side slip may induce rolling moment as more of the propwash passes over one wing than another. In addition to this when slipstream strikes vertical fin, force is generated at the aerodynamic center of fin which is at the vertical and longitudinal distance from aircraft's center of gravity. Slipstream therefore generates positive rolling moment in addition to negative yawing moment. A list of undesired moments with their primary effect along remedial action is mentioned in table 1.

III. EFFECTS OF PROPELLER HELICAL SYSTEM ON AIRCRAFT AERODYNAMIC FORCES
Propellers can be attached with an airplane in two basic configurations, tractor or pusher. Tractor configuration pulls an airplane, whereas pusher configuration pushes them through fluidic medium. Wing to propeller interactions exist in pusher configuration whereas propeller to wing interactions exist in tractor configuration, refer figure 7. It cannot be said with certainty, that which propeller configuration suits best for flying wing bodies. Both configurations have advantages and disadvantages associated with them.
The slipstream of pusher configuration propeller does not interact with aircraft wings which eliminate their interference effects, but still, there are some inherent problems. The main concern is the ingestion of highly turbulent airflow coming from fore body and wings into the propellers which deteriorates propeller efficiency and performance. There is an increase in noise level when pusher propeller cuts through the wings wake. Backward CG shift gives problem in longitudi- nal stability whereas ground clearance is yet another problem [21]- [23].
Tractor propeller does not give problems in ground clearance, CG shift or noise levels but wing aerodynamic efficiency is directly compromised since swirl velocity, vorticity and helicity destroy laminar nature of air flow over the wings [24]. Previous researchers never selected one particular configuration and therefore relied upon extensive computational or experimental studies for propeller / wing interactions. The task of configuration selection becomes even more challenging when low aspect ratio wings are integrated with propellers. Periodic nature of viscous propeller wake covers the majority of wing area, thereby affecting wing aerodynamic performance [25].
Two schools of thought have been reported in the literature. Some researchers reported an increase in lift along with lift curve slope increment, delay in stall due to flow separation and a decrease in drag in case of tractor configuration. However, some other researchers reported the same benefits in case of pusher configuration. Table 3 mentions major conclusions from available literature from which no conclusive selection can be adopted.
The literature of propeller interactions with low aspect ratio wings like micro aerial vehicles is found with mixed views from researchers' side. Dave Witkowski [39] and RT Johnson [40] found, that the drag is reduced and the lift is increased in tractor configuration. With the advance ratio, wing drag coefficient increased whereas thrust and power coefficient decreased. Increased axial velocity along with swirl velocity can enhance wings aerodynamic efficiency. However, Sergey [29] found, that maximum figure of merit could be obtained in pusher configuration which Gamble negated [32]. FOM is usually defined for helicopter rotor during hover phase. FOM is the ratio of ideal power required for hover to actual power required. LLM Veldhuis in his PhD thesis analyzed tractor configuration of propeller with low aspect ratio semispan wing at Delft University Wind Tunnel [41]. He concluded, that a strong reduction in wing induced drag was noticed VOLUME 8, 2020 TABLE 1. Undesired moments on right-handed propeller aircraft with remedial actions. Secondary effects to unwanted moments are also stated. It is highlighted that roll causes yaw in opposite direction whereas yaw causes roll in same direction [20]. Physical phenomenons are extracted for airplanes which have large propeller diameter to wingspan ratio, D b of around 12% or more Recommended configuration has negative C mα , highest α stall , least C D0 , highest L D , highest α L/D max , highest C Lmax for a fixed wing micro aerial vehicle. Tractor configuration was identified as most suited option and therefore, adopted for analysis purposes in this research once propeller was installed at a negative tilt angle with respect to the wing. Catalano concluded, that the pusher configuration affects wing aerodynamic characteristics by delaying boundary layer separation [42]. W Null did research on propulsive-induced flow from the tractor propeller on the aerodynamics of micro air vehicles at higher angles of attack [26]. He found, that induced flow from the propeller resulted in delayed stalling behavior with a decrease in the lift to drag ratios at low angles of attack. From open source literature on propeller to wing interaction, it was found, that only a few researchers conducted experimental work on propeller effects for low aspect ratio wings concerning both configurations. Thipyopas found, that the stall angle was delayed and maximum lift improved in both configurations [27] and [28]. He determined, that the increase in maximum lift coefficient associated with a high stalling angle was more pronounced for biplane airplane, as compared to single winged micro aerial vehicle. Sergey concluded, that 20% more thrust could be obtained for the pusher propeller for the same power input as compared to the tractor propeller [29]. Jon Ahn concluded, that compared to the pusher propeller, the tractor propeller increased lift to drag ratio along with an increase in pitching moment [17]. Gavin Kumar Ananda also compared both configurations and commented that the lift curve slope varied in tractor configuration, which was not observed in pusher configuration [30]. Chinwicharnam found, that delay in stall angle was linked with energizing air flow by the propeller propwash. The impact of propwash on aerodynamic efficiency was larger as compared to wing wash for low aspect ratio wings. He concluded, that tractor configuration was more preferred as far as aerodynamic efficiency was concerned [31]. Other researchers did experimental or computational work on one of the configuration. Gamble conducted experimental work on tractor configurations and found, that rolling moment of about 38% of motor torque was induced by the propeller [32]. Sungjin Choi found, that aerodynamic efficiency metered by lift to drag ratio decreased by 2% in the case of pusher propeller on MAV [16]. Arivoli Durai found significant effect of propeller slipstream during experimental studies of MAV on lift, drag and stall angle once the propeller was installed in tractor configuration [33]. Deng found, that at lower AoA, micro aerial vehicle was submerged with propeller vortex ring whereas, at higher AoA, wing tip vortex structure had a major contribution to flow physics [34]. K A Kasim found that pusher propeller configuration enhanced vortex structure better than tractor configuration [43]. Pampala employed CFD techniques to study propeller thrust, power and slipstream characteristics whereas Premkumar did experimental work on the efficiency of MAV propeller [44] and [45]. Researchers found, that a thrust calculation became significant for larger diameter propellers. LW Traub studied pusher propeller slipstream effect in the wind tunnel for delta wing of 65 • leading edge sweep. He concluded, that pusher configuration could cause a delay in vortex breakdown consequently increasing aerodynamic efficiency. Propeller increased suction level on leeward side which increased lift at all angles of attack. Lift was favored at low advance ratios [46]. Shuvrangshu Jana modeled propeller-induced flow effects as a function of motor rotation speed and mathematical analysis was performed to quantify their effects [19]. He concluded, that propeller flow contributed to the rolling moment and to the pitching moment, while it had negligible effects on the yawing moment. Ahmed Aboelezz conducted Wind tunnel testing with propeller effects and a non-linear flight simulation model was proposed for true prediction of MAV dynamics [47]. Chen improved lift to drag ratio ratio of the wing by conducting optimization study on propeller slipstream using Genetic Algorithm and Kriging surrogate modeling technique for distributed electric propulsion [48]. Favorable propeller interference was used to promote aerodynamic performance of double wing by Hongbo [49]. Numerical simulations were carried out to analyze complex flow structures around the wing. He concluded, that suction upstream of the propeller and acceleration downstream of the propeller can be used to improve wing performance [49]. Much of the work has been done in the past relating micro aerial vehicles flight dynamics and researchers have been employing either computational or experimental techniques to model correct forces and moments. However, only propeller rotations had been found significant for calculations of forces and moments [19]

IV. NOVELTY OF RESEARCH
In this research we discuss the effects of tractor propeller induced flow effects on overall flight dynamics of a low aspect ratio micro aerial vehicle. In the case of micro aerial vehicle, ratio of propeller diameter to wingspan D b , is significantly large which specifies large effect of propwash on the vehicle aerodynamics. For this research, large D b is taken as 0.2 or higher. A major portion of micro aerial vehicle wings are submerged in a highly rotating vortex flow emerging from propeller, which has high levels of swirl velocity, vorticity and helicity. This helical flow structure causes aerodynamic instabilities and fluctuations in aerodynamic coefficients, which is not fully understood until today. Though propeller effects as a function of propeller rotational speed have been reported earlier by Jana [19], but our research focuses on the propeller effects on overall aerodynamics of vehicle as a function of advance ratio, J . Contribution of each factor of advance ratio on aerodynamic coefficinets is discussed. Consequently, 6-DOF equations of motion, valid for longitudinal trim point corresponding to angle of attack of 2 • are proposed for flying wing micro aerial vehicle.
The paper is managed in the following manner. In section V, FWMAV geometry and wind tunnel testing description is discussed. Section VI discussed extensive discussions on the results of wind tunnel testing. Section VII discusses proposed 6-DOF model with inclusion of advance ratio terms. Section VIII summarizes all results before presenting conclusion in section IX, recommendations in section X and future work in section XII.

V. FLYING WING MICRO AERIAL VEHICLE GEOMETRY AND WIND TUNNEL TESTING
FWMAV is flying wing micro aerial vehicle of fixed wing type as shown in figure 10. FWMAV was designed with a high leading edge sweep angle of 40 • . A dihedral angle of 2 • was proposed for lateral stability. Detailed design features of FWMAV is shown in authors' earlier work [50].

A. DESCRIPTION OF WIND TUNNEL
In this research, a subsonic wind tunnel manufactured by Aerolab, USA was used to conduct wind tunnel experimentation (refer figure 8). The wind tunnel is a closed circuit, continuous return, rectangular test section and horizontal type tunnel which is capable of generating linear velocity of 110 m/s at an atmospheric pressure. Maximum Reynolds number is 8.32 × 10 6 based on 1.00 m of characteristics length. Subsonic wind tunnel has a rectangular test section of 2ft × 3ft × 6ft and equipped with state of the art digital data acquisition system (DAQ) for measurements of aerodynamic forces and moments caused by these forces. During experimental testing of FWMAV, accurate results from wind tunnel are obtained, because no scaling effects are encountered during data reduction.Test section of subsonic wind tunnel accommodated full scale FWMAV. Reynolds number and Mach number, both were matched with actual flight conditions which was 2 • pitch attitude which is corresponding to cruise flight of FWMAV. Motor drive, flow conditioning, test section, and control console constitute major parts of subsonic wind tunnel. Flow conditioning consists of stilling chamber where stagnation pressure is stored, diffuser section VOLUME 8, 2020 and contraction section. 150 HP electric motor drives variable pitch counterclockwise rotating fan, through which linear velocity can be changed and adjusted for wind tunnel tests. Fan speed up to 1500 RPM can be achieved which corresponds to 110 m/s of linear velocity. Honer comb screen is placed inside stilling chamber to achieve attenuation of flow turbulence. Fixed contour of converging section accelerates the flow from stilling chamber to test section. Since converging section of subsonic wind tunnel is fixed, therefore variation in linear velocity inside test section is achieved through pitch variation of CCW rotating fan which has four impellers. Linear velocity inside test section, pitch attitude and yaw attitude of test object is controlled by control console. Pitch attitude of ± 30 • and yaw attitude of ± 90 • can be achieved inside test section through control console.

B. WIND TUNNEL SUPPORT SYSTEM AND ITS EFFECTS
Subsonic wind tunnel has six strain gauge sting balance as well as pyramidal balance for installation of test models inside the test section. Extensive experimental calibration of pyramidal balance was carried out using loading in forces direction and loading in moments direction to proof load the balance and to ascertain component sensitivity. Balance calibration matrix, calibration slopes and interaction constants were evaluated and used to purify results reflected by data acquisition system. The details of complete calibration process along with calibration matrix values are shown in [50]. Since FWMAV was placed over pyramidal balance with the help of a supporting system, therefore, this affects overall drag measured by data acquisition system of subsonic wind tunnel (refer figure 9). Drag of supports is called ''tare drag'' whereas variation in air flow pattern due to the presence of these supports is called ''interference drag'' [51]. Correct estimation of tare and interference drags is a complicated process but it is explained here for clarity purposes. The process can be thought of a three step process. In first step, test object is attached to pyramidal balance in normal position and supports are attached to its lower surface. The drag D 1 is noted (refer equation 3). In second step of the process, the test object is attached to pyramidal balance in an inverted position and supports are attached to its upper surface only. The drag D 2 is noted then (refer equation 4). In third step, supports are added on lower surface of test object while keeping it in an inverted positions like in step two. The drag D 3 is noted (refer equation 5). After mathematical operations on these equations, drag on test object is estimated (refer equations 6 and 7). The drag on test object, D 5 is the drag which is free from tare and interference drag. Same procedure was adopted while evaluating drag on FWMAV for both types of wind tunnel tests conducted in this research, angle of attack sweep and tests with varying propeller RPMS at 2 • pitch attitude.
. Support System inside Wind Tunnel Test Section. The support system is used to position test objects at correct attitude for conduct of wind tunnel tests. Tare and Interference drag generated by these supporting structure is estimated and removed from overall drag of the object using equations 3 to 7.

C. WIND TUNNEL CORRECTIONS
Wind tunnel boundary corrections for three dimensional bodies in the form of horizontal buoyancy, solid blocking, wake blocking and turbulence factor were applied to obtain final results [52]. Total blockage factor denoted by was calculated and the following expressions were used to determine corrected parameters of velocity, dynamic pressure and Reynolds number, refer equations 8, 9 and 10. 'u' denotes uncorrected values and 'TF' denotes tunnel turbulence factor. Chordwise downwash correction factor denoted by τ and boundary correction factor denoted by δ were also used to apply correction for streamline curvature in the coefficient of forces and moments. As an example, corrected coefficient of lift denoted by C L is shown in equation 11. The values of δ and τ are taken from wind tunnel testing by Pope [52].
Spatially standardized air velocity inside the test section is a prerequisite for authentic subsonic wind tunnel testing. Several parameters need to be quantified like flow angularity, noise level inside the test section, vortices generated from walls, pressure variation along the test section floor. Among these, velocity fluctuations along streamwise direction was considered most critical and therefore measured. Axial velocity was measured in the center, 200 mm downstream of the test section with hot wire anemometer. Hot wire anemometery system provided by DANTEC R dynamics was used to determine turbulent intensity inside test section of subsonic wind tunnel. CTA probe was made from tungsten wire with gold plated ends. Data was recorded using 16-bit, 4-velocity channel NI 9215 A/D board with a sampling frequency of 200kHz. Turbulent intensity as recorded by single wire probe miniCTA 54T42 type 55P16 is shown in figure 11. Turbulence intensity is equivalent to the ratio of mean velocity fluctuation to mean velocity and recorded as percentage. It can also be defined as the ratio of root mean square of turbulent velocity fluctuations and mean velocity. Turbulence intensity of 0.25% was recorded for subsonic wind tunnel once it was operated at 20 m/s, refer figure 11. Experimental uncertainty assessment has been carried out to measure uncertainty in the measurement of forces and moments as a function of repeatability of results. Forces and moments are calculated for three times against one wind tunnel test and error was evaluated. Error estimates are as   A total of 27 wind tunnel tests were planned and conducted. Forward velocity, propeller rotational speed and propeller diameter variation was recorded during these tests. Windmilling tests were also included and considered as a baseline for calculation purposes. Detail description of wind tunnel tests along with their ID numbers is shown in table 4.

TABLE 4.
Wind tunnel test descriptions and settings. Three wind tunnel speeds, three propeller rotational speeds and three propeller diameters were selected for 27 advance ratios. 9 wind mill tests were also planned for every propeller and every wind tunnel speed. A total of 36 wind tunnel tests were planned and conducted. 12 tests were conducted for complete angle of attack sweep, whereas remaining tests were conducted at longitudinal trim angle of attack of 2 • as shown against each. Longitudinal trim angle of attack of 2 • was determined experimentally in early research of authors [50]

VI. RESULTS AND DISCUSSIONS
The coefficient of lift, C L , was plotted against angle of attack for wind tunnel tests WM 5−20, T 3, T 6, T 9, WM 6−20, T 12, T 15, T 18, WM 7−20, T 21, T 24, T 27 and shown in figure 15. It was found that lift increases with a decrease in advance ratio. The decrease in advance ratio was obtained through an increase in propeller diameter and propeller rotational speed. Wind milling lift coefficient was found lesser in all tests invariable to advance ratios. Post stall region was found with more spread for the cases of 6 inch and 7 inch propellers, refer figures 15(b) and 15(c). No appreciable change with advance ratio was noticed for 5 inch propeller in post stall region, refer figure 15(a).
Like coefficient of lift, the coefficient drag also followed conventional trend of variation with angle of attack. The coefficient of drag was found higher at lesser advance ratio, refer figure 16. It was determined that higher propeller rotations and larger propeller diameter increased drag significantly. Pyramidal balance calibration graph. This graph shows calibration slope between applied load in lift direction (LF) and lift force reading (LFR) from data acquisition system of wind tunnel. Relation of voltage from strain gauge as shown by system with applied load is also shown. Further details on calibration matrix and interaction constants is shown in [50].
The coefficient of drag for wind milling cases was found least in numbers. A slight irregularity was noticed for 7 inch propeller at high angles of attack once tested at higher propeller rotational speeds, refer figure 16(c). It was determined that high swirl velocity and strong helical flow structure could be the most probable cause of this irregularity specially at higher angles of attack.
Coefficient of lift to drag ratio, C L C D , was found maximum for wind milling case. Maximum L D decreased with an increase in propeller diameter. Maximum L D value was noticed at the angle of attack of 3 • for all propeller diameters, refer 17.
Since FWMAV was trimmed at longitudinal angle of attack of 2 • during wind tunnel testing as shown in [50], therefore aerodynamic forces and moments due to forces were extracted at angle of attack of 2 • and plotted against advance ratios, refer figures 18 and 19. It was noticed that all coefficients of forces decreased with an increase in advance ratio. The gradient of lift coefficient was found larger as compared to the gradient of drag coefficient and side force coefficient. It was highlighted, that decrease in advance ratio was achieved either by an increase in propeller diameter or by an increase in the propeller rotational speed. Pitching moment coefficient varied linearly with advance ratio with a slope of −0.03177, figure 19. Yawing and rolling moment coefficients reflected nonlinear behavior with advance ratio, figure 19. Their dependency on advance ratio was calculated by taking derivative of polynomial equation and putting desired value of J. For J = 0.15, ∂C n ∂J becomes −0.1178 and ∂C l ∂J becomes −0.1425.
Aerodynamic coefficients of forces and moments were plotted for varying propeller diameter and propeller rota-      gradient. For moments, there was a negligible difference with propeller diameters of 4 inch and 5 inch, However, all three moment coefficients increased sharply when 7 inch propeller was used, refer figure 22(a), 22(b) and 22(c). This reflected that the larger propeller induced higher rolling and yawing moments, whereas the effect of propeller diameter on pitching moment was insignificant. Pitching moment coefficient was found more dependent on propeller rotational speed rather than propeller diameter, refer figure 22(a). Maximum of force and moment aerodynamic coefficients was found with no definitive trend, however wind tunnel tests ID T 9, T 18 and T 27 were identified for maximum rolling and yawing moment coefficients. It was determined that increased Aerodynamic efficiency parameter, denoted by L D was plotted to identify the best combination of propeller diameter and propeller rotational speed in figure 23. It was noticed, that aerodynamic efficiency initially increased from 5 inch propeller to 6 inch propeller and then decreased once 7 inch propeller was used. Maximum of L D was noticed at propeller diameter to wingspan ratio, D b , of 0.338, which corresponded to 6 inch propeller. At this propeller diameter, aerodynamic efficiency was found to decrease with propeller rotational speed. The Maximum of L D was noticed at 50% throttle instead of 100% throttle. The maximum value of L D was noticed for wind tunnel tests ID T 11, refer table 5. With larger propeller diameter to wingspan ratio, coefficients of lift and drag, both increased but the gradient of drag increment with propeller diameter was more significant. Therefore maximum L D was observed at forward velocity of 15 m/s, 0.338 of propeller diameter to wingspan and minimum throttle setting, which corresponded to wind tunnel test ID T − 11 and had an advance ratio of J = 0.120. L D was found to be the nonlinear function of propeller diameter and wingspan ratio and propeller rotational speed in form of RPM, refer figure 23 and equation 15.

VII. PROPOSED 6-DEGREE OF FREEDOM MODEL WITH INCLUSION OF ADVANCE RATIO TERMS
Six degree of freedom model presents mathematical expressions for modeling flight dynamics of a vehicle. These are set of six nonlinear differential equations which represents linear line) and gravity force acts at center of mass of FWMAV. Since propeller rotational axis is aligned with FRL, therefore, thrust force did not generate any moment around any axis. Six equations for 6-DOF model are stated as 16 till 21, refer [53].
For low angles of attack, X-force can be estimated with − C D Sq ∞ and Z force can be estimated with − C L Sq ∞ . These forces are proposed with an additional factor concerning advance ratio J , which were determined experimentally in this research. It is highlighted that X-force and Z-force consist of aerodynamic and propulsive forces, both [54]. where Here, C DJ determines change in the coefficient of drag with change in the advance ratio, C DJ = ∂C D ∂J . Similarly, VOLUME 8, 2020 Longitudinal coefficients like C L , C D and C m are dependent upon angle of attack, α, where as later-directional coefficients are dependent upon side slip angle, β. In literature these coefficients are known as static coefficients. Static coefficients are calculated experimentally using wind tunnel tests while keeping test model stationary. Other derivatives are estimated when test model maneuvers either with slow rate of change of attitude known as "rate derivatives" OR when unsteady aerodynamic effects become significant as "acceleration derivatives". Most prominent rate derivatives are C lp , C mq and C nr where as notable acceleration derivative is translation acceleration derivative, C mα .
Wind tunnel tests were conducted for very low free stream velocities, therefore no appreciable change in aerodynamic coefficients was noticed once forward velocity was changed. Although change in forces magnitude and moments magnitude was registered during wind tunnel tests, but due to increase in propeller induced flow velocity, aerodynamic coefficients did not reflect any significant change, as explained in equations 29 and 30. Force coefficients with proposed advanced ratio derivative are shown in equations 23, 24 and 25. Proposed moments coefficients with the addition of change in moment due to change in advance ratio terms are shown as equations 26, 27 and 28.
A. ACCELERATION DERIVATIVES Acceleration derivatives are important parameters specially for vehicles which undertake maneuvers with high maneuverability rates or when unsteady aerodynamics has a prominent role to play like flutter. Once a dynamical system like an aircraft is written in state space format, mass matrix is modified to include acceleration derivatives. Overall effect of acceleration derivative is to alter the apparent mass and inertia properties of the dynamical system. For this reason, acceleration derivatives are usually referred to as virtual or apparent mass. Physically, whenever any flying object moves through an air, the surrounding air mass is entrained which moves with the moving object. The mass and inertia of entrained air modifies inertia of airplane. Acceleration derivative quantifies this change [53]. For small airplanes, mass of entrained air is small, therefore acceleration derivatives become insignificant. For most of the modern fighter aircraft which can change their orientation very fast, translational acceleration derivative may be significant. Second physical effect is the transient disturbance in the wing downwash field which interacts with horizontal tail of aircraft with some time delay and causes a change in incidence angle of tailplane, a short time later. This particular phenomenon is known as downwash lag effect. For this reason, in many text books, the translational acceleration derivative is assumed to be zero for tailless airplane, refer [55]. Third physical explanation ofα derivatives is the settlement of pressure distribution to an equilibrium value on wing and on tail once angle of attack of the airplane is changed suddenly. Calculation of this physical effect comes under the purview of Unsteady aerodynamics as acceleration derivatives. All other derivatives of airplane can be calculated on the basis of steady state aerodynamics. Flutter phenomenon of large aspect ratio wings is a classical case where acceleration derivatives becomes so important that they can not be ignored [56].
Translation acceleration derivative can be obtained either by flight tests, wind tunnel testing or computational fluid dynamic analysis, refer [57]. Wind tunnel testing for determination of translational acceleration derivatives is time consuming and a very costly solution. Subsonic wind tunnel used in this research is not equipped with mechanism of generating step or impulse acceleration on test object; therefore wind tunnel tests could not be performed. With current advancements of computational power, CFD can be the only viable solution for evaluation of acceleration derivatives. Computational evaluation of dynamic derivative is a cumbersome process and current research is not directly focused towards their evaluation. Readers are requested to consult [58]- [62] for more details on estimation and evaluation of acceleration derivatives through application of CFD. In this research however, it is assumed that C mα has negligible effect on over all flight dynamics of FWMAV as compared to other factors, and therefore, regarded as zero in 6DOF modeling. It is however appreciated by all the authors of this research, that flying machines designed even without horizontal tail do not have zero translation acceleration derivative as explained by Baigang M in [63].

B. RATE DERIVATIVES
Several methods exist for the evaluation of aircraft rate and acceleration derivatives [64] and [57]. Potential solvers is one way which solve Laplace equation valid for inviscid, incompressible and irrotational flow [65]. AVL (Athena Vortex Lattice) [66], PANAIR (Panel Aerodynamics) [67], Tornado [68] and XFLR [69] are few renowned software out of many. In this research, subsonic wind tunnel is not equipped with dynamic pyramidal balance, therefore important rate derivatives (C Lq , C mq , C lp , C lr , C nr and C np ) were estimated using potential flow solvers (XFLR and Tornado) along with existing theoretical correlations with moderate to medium level of confidence (figure 25) [70].
Potential flow solvers use Vortex Lattice Method (VLM) for calculations of forces, moments and aerodynamic coefficients for airplane at low speeds and low angles of attack [71]- [73]. VLM is capable of solving complete flight dynamic problem with full linearization of the aerodynamic model about any trim point. VLM methodology works on vortex filament with vortex circulation strength, . Lifting surfaces are modeled using vortex filaments using Helmholtz's first, second and third theorems, [74] and Kelvin's circulation theorem (refer equation 31) [75]. VLM is implemented on 3 dimensional lifting surfaces for which interaction with arbitrary point is described by Biot-Savart Law (refer equation 32) [76]. Kutta-Joukowski theorem relates circulation produce by closed 2D body and aerodynamic force that body experiences per unit span (refer equation 33). In VLM lattice of panels define lifting surface where every panel has a horseshoe vortex and a control point. Bound vortex on every panel is placed at c 4 whereas control point is defined as 3c 4 of panel chord length. Influence coefficient matrix is found by the influence of these vortices on all the control points of every panel. Flow tangency condition is forced to satisfy and strength of vortices of all panels is obtained (refer equation 34). Using vortices strength, forces on all panels are calculated which could be converted into pressure using unsteady Bernoulli's equation (35).
However, in this research, potential solvers used Unsteady Vortex Lattice Method (UVLM) to evaluate rate derivatives [77] and [78]. The UVLM code provides a medium-fidelity tool for the estimation of dynamic aerodynamic loading. UVLM has been found very beneficial in quick estimation of derivatives where unsteady aerodynamics of lifting surfaces undergoes complex kinematics [77], [79]. In UVLM, Wake effects are added in the form of free wake or prescribed wake. Mathematical details along with formulation of wake models can be found in [80]- [83].
XFLR was developed by Mark Drela with an earlier version of XFOIL [69] whereas Tornado was improved for propeller effects by Tomas Melin [84]. XFLR did not include propwash effects and therefore software like QBLADE needed to be used for propeller effects independently [85]. In order to estimate rate derivatives of FWMAV with propeller effects, Tornado was selected as a primary software. Tornado has been tested and found to give satisfactory results for propwash affected vehicles [86]. Vortex Ring Modeling (VRM) and the Helical Wake Modeling (HWM) are the two approaches which this software can employ to couple propwash effects on vehicle aerodynamics, refer [87], [88] and [89]. In this work, Helical Wake Modeling approach was used to evaluate rate derivatives of FWMAV inclusive of propwash effects. A propeller diameter of 5 inch, RPM of 7800 and a free stream velocity of 20 m/s was used as propeller advance ratio conditions.
From the geometrical details of FWMAV, elevons tail volume ratio in elevator role, V H was 0.340, while tail volume ratio for winglets was 0.0411 (refer figure 25). For an estimation of elevon lift curve slope, wing lift curve slope of was used which was determined experimentally through wind tunnel tests [50]. For an estimation of winglets lift curve slope, an inviscid airfoil lift curve slope of 2π was corrected for winglets leading edge sweep of 30 • and aspect ratio of 0.98 as proposed by Helmbold's equation 36, taken from [90]. In equation 36, AR is aspect ratio, LE is leading edge sweep angle of winglets. Using analytical expression shown in equations 37 till 42, rate derivatives were calculated and shown in tabletablerates. While comparing analytical results with potential flow solvers in table 6, few differences are observed. These differences could be for many reasons, however few reasons are mentioned below which could affect results obtained from UVLM code; • Thickness of lifting surface is ignored. • Small angle approximations. VOLUME 8, 2020 • Low aspect ratio effects. • VLM valid for potential flows only. • Prediction of vortices from wing leading edge. • Formulation to capture unsteady aerodynamics effects. • Formulation to capture correct propwash effects These derivatives were utilized in proposed 6DOF modeling while appreciating known differences and shortcomings. However for more trusted and dependable values of rate and acceleration derivatives, new computational procedure as proposed by Bai-gang M may be adopted [63].
After implementation of small perturbation theory at trim point (all staes zero except u velocity and pitch angle), equations of motion for small perturbation were written. Forces

VIII. SUMMARY OF RESULTS
Following results are summarized from wind tunnel testing conducted in this research.
1) Lift coefficient and drag coefficient both increased with a decrease in advance ratio at all angles of attack. The trend remained same for every propeller diameter (figures 15 and 16). The probable cause is an increase in propeller tangential velocity and induced velocity components with respect to propeller RPMs. Since flow over the wings increased with propeller rotational speed (reducing advance ratio, J), therefore lift and drag coefficients, increased with angle of attack by decreasing advance ratio. 2) Post stall region of lift coefficient reflected variation with advance ratio for 6 inch and 7 inch propellers, whereas this trend was not observed for 5 inch propeller, refer figures 15(b) and 15(c). The spread in the curve in post stall region is caused by aggressive flow separation at high angles of attack once large propellers were used. 3) Drag coefficient followed a conventional trend of variation with angle of attack. However, slight irregularity was observed at higher angles of attack, (refer figure 16). The most probable cause is the increase in wake region which consequently increased pressure drag, once angle of attack was increased. 4) Maximum L D was found at angle of attack of 3 • for all tested advance ratios. Wind tunnel tests conducted in wind mill conditions reflected maximum value of L D . Most probable cause of increase in L D is reduced drag coefficient during windmill cases. 5) Among three different sizes of propellers, the smallest propeller of 5 inch diameter generated maximum value of L D as compared to 6inch and 7 inch propellers (refer figure 17). 6) For wind tunnel tests conducted at longitudinal trim angle of attack of 2 • , it was found that • All force coefficients reflected linear variation with advance ratio, where lift coefficient reflected major change (refer figure 18).
• All moment coefficients reflected variation with advance ratio. Pitching moment coefficient reflected linear variation whereas yawing and rolling moment coefficients reflected non-linear variation with advance ratio. Significant amounts of negative rolling moment and negative yawing moments were recorded at low advance ratios, which corresponded to high propeller rotations with large propeller diameter (refer figure 19).
• L D was found to be a non-linear function of advance ratio, D b and propeller RPM (refer figure 23). Maximum L D was found for wind tunnel test ID T 11, which utilized 6 inch propeller ( D b =0.338) with 50% propeller speed at 15m/s of forward velocity (refer figure 20).
• Side force coefficient and pitching moment coefficient reflected least dependency on advance ratio. Since the propeller was installed on FRL therefore, propeller slipstream was found symmetric in x-y plane of FWMAV, and hence did not change pitching moment to a considerable extent. No concrete justification of minor side force coefficient change with advance ratio could be found (refer figures 18 and 19). VOLUME 8, 2020 • A major contribution of forces and moments caused by propeller slipstream is proposed in the form of C DJ , C LJ , C YJ , C lJ , C mJ and C nJ . These values against trim angle of attack of 2 • is shown in table 7. 7) All above findings were mathematically captured and an improved 6-DOF model has been proposed which has contribution of aerodynamic coefficients variation with advance ratio. The set of six improved 6DOF equations of motion are shown as equations 43, 44, 45, 46, 47 and 48. Contribution of advance ratio is found to be significant and can not be ignored for realistic 6 DOF modeling.

IX. CONCLUSION
Aerodynamic characteristics of flying wing micro aerial vehicle under propeller slipstream were investigated through wind tunnel experimentation at various advance ratios. Left turning tendencies for right-handed fixed pitch propellers were quantified through exhaustive wind tunnel experimentation to relate aerodynamic coefficients variations with advance ratio. It was found, that force coefficients increased with a decrease in an advance ratio. Aerodynamic efficiency parameter, L D , reflected maximum values for wind mill cases which decreased with propeller rotations. The FWMAV was also tested in subsonic wind tunnel at longitudinal trim angle of attack of 2 • . It was observed, that lift coefficient and drag coefficient showed high dependency on advance ratio as compared to side force coefficient. Pitching moment coefficient though varied linearly with negative gradient, but variation was not found significant. Rolling and yawing moment coefficients, both showed a non-linear behavior. A significant amount of negative rolling moment and negative yawing moment was observed in right-handed propellers at low advance ratios. It was determined that left turning dependencies for right-handed propellers are related to helical nature of flow in prop wash. The dependencies become significant at low advance ratio that occur at high propeller rotations and with large propeller diameters. Research concluded, that aerodynamic force and moment coefficients vary significantly with change in advance ratio specially, for vehicles with large propeller diameter to wingspan ratio, D b . The research proposed an improved 6-DOF model with included effect of propeller slipstream on low aspect ratio wing for better estimation of flight instabilities.

X. RECOMMENDATIONS
In this research wind tunnel experiments were limited by angle of attack sweep, however it is recommended that effect of side slip angle may also be investigated. FWMAV was designed with large winglets to provide lateral stability. Effect of lateral derivatives on side slip angle would be significant. It is therefore recommended that exhaustive wind tunnel experimentation may be planned and dependencies on side slip angle may also be determined. Further, in this research rates and acceleration derivatives were estimated using poten-tial flow solvers which presented errors in comparison with analytical results. It is recommended, that CFD analysis may be carried out to determine exact values of rate and acceleration derivatives for FWMAV by adopting new CFD methods as proposed in [63] and [57].

XI. FUTURE WORK
This research discussed variation of L D with propeller diameter to wingspan ratio, D b by varying propeller diameters only while keeping wingspan constant. However, it is suggested that analysis may also be carried out to study variation of L D by varying D b through variation of wingspan for fixed propeller diameter. propeller induced velocity δ boundary correction factor δ e elevon deflection angle τ chordwise downwash correction factor C lα airfoil section lift coefficient variation with α C Lα lift coefficient variation with β C Dα drag coefficient variation with α C Y β side force coefficient variation with β C lβ rolling moment coefficient variation with β C mα pitching moment coefficient variation with β C nβ yawing moment coefficient variation with β C Lq lift coefficient variation with dimensionless pitch rate C mq pitching moment coefficient variation with pitch rate C Lδe lift coefficient variation with δ e C Dδe drag coefficient variation with δ e C Y δe side force coefficient variation with δ e C lδe rolling moment coefficient variation with δ e C mδe pitching moment coefficient variation with δ e C nδe yawing moment coefficient variation with δ e C lp rolling moment coefficient variation with roll rate C lr rolling moment coefficient variation with yaw rate C np yawing moment coefficient variation with roll rate C nr yawing moment coefficient variation with roll rate TAIMUR ALI SHAMS was born in Lahore, Pakistan, in 1979. He received the B.E. degree in aerospace engineering from the National University of Sciences and Technology, Islamabad, Pakistan, and the M.S. degree in fluid dynamics from Air University, Islamabad. He is currently pursuing the Ph.D. degree with the College of Aeronautical Engineering, NUST, and working on flight instabilities of fixed wing micro aerial vehicles. His industrial experience includes two years of managing organizational level and five years at intermediate level aircraft maintenance. He has a vast experience of carrying out Research and Development projects concerning modifications on aircraft. His research interests include experimental aerodynamics, wind tunnel testings, flight dynamics, 6DOF modeling, and evaluation of flight test data.

APPENDIX COMPLETE LIST OF NOMENCLATURE
SYED IRTIZA ALI SHAH graduated from the Georgia Institute of Technology, USA. He received the B.E. degree from NED University and the M.S. degree in aerospace engineering from NUST. He then proceeded abroad to Georgia Tech, from where he completed his M.S. degree in flight mechanics and controls, his M.S. degree in machine vision and robot controls, and his Ph.D. degree in aerial robotics. Since then, he has been teaching at the National University of Sciences and Technology, Pakistan. He has authored more than 100 research articles and had supervised 15 Ph.D. and M.S. students and a number of B.E. students from aerospace engineering, mechanical engineering, manufacturing engineering, robotics, intelligent machine engineering, and biomedical engineering disciplines. He also had been the Head of NUST community service program from 2011 to 2015. He designed a course on Community Service Learning for B.E. and M.S. students of NUST, which has now been made compulsory part of NUST curriculum. During his tenure, he raised the NUST Community Service Club from 200 volunteers to 4500 volunteers and taught, trained and supervised over 5000 students in community service.
AAMER SHAHZAD received the bachelor's degree from CAE, NUST, Risalpur, in 2004, with a specialization of aerospace engineering, the master's degree in computational science and engineering (CS and E) from the Research Center for Modeling and Simulation (RCMS), NUST, Islamabad, Pakistan, and the Ph.D. degree in mechanical engineering from the University of New South Wales (UNSW), Canberra, Australia. He is currently an Assistant Professor of aerodynamics with the Department of Aerospace Engineering, CAE, NUST. His industrial experience includes five years of managing organizational level and Intermediate level aircraft maintenance and two years of managing depot-level maintenance in an Engine Overhaul facility. Since 2010, he has been involved in the research of bio-inspired and bio-mimicking flapping wings, with an emphasis on improving aerodynamic performance in hovering flight. The primary focus of his research has been to study the effects of wing geometric parameters, kinematics, and wing corrugations and flexibility on the aerodynamic performance of flapping wings.
ALI JAVED received the bachelor's degree in aerospace engineering from CAE, NUST, Islamabad, Pakistan, in 2001, the master's degree in solid mechanics from Air University, Islamabad, and the Ph.D. degree in fluid structure interaction from the University of Southampton, UK. He is currently an Associate Professor of astronautics with the Department of Aerospace Engineering, CAE, NUST. He is involved in active research and development activities relating to solid mechanics, structural dynamics, radial basis functions, mesh free particle methods, and computational modeling of flow induced vibrations.