Aerodynamic Analysis of a Light Aircraft at Different Design Stages

During the evolution of an airplane aerodynamic design, proper calculation methods and software tools should be utilized, which correspond to the airplane category and project development level. In case of light aircraft, the general trend is the application of analytical and semiempirical methods at the initial stages, combined with simplified inviscid CFD computational models, and fairly complex viscous CFD analyses at higher design levels. At the present stage of light aviation development, it is assumed that the contemporary design tools for each of those steps should be appropriate enough, so that they actually verify and additionally fine-tune each other's results. This paper describes the calculation tools and methods applied during the aerodynamic analyses of a new light aircraft at different development stages, and compares the results obtained by them, with the aim to verify and support the above statement, considering light aircraft aerodynamic design.


INTRODUCTION
The Innovation Center of the Faculty of Mechanical Engineering, University of Belgrade, has recently been involved in the development of a new light aircraft (NLA -see Fig. 1).The posted design requirements were: a single-engine airplane with two seats, which satisfies the Aerobatic and Utility requirements of EASA CS 23 regulations, with metal primary structure and plastic composite engine cowling, wing and tail tips and fillets, and with fixed tricycle landing gear, aimed for primary and advanced VFR and IFR flight training, sport flying, glider towing, etc.
Within the design process of the NLA, the initial aerodynamic calculations have been performed by Datcom method [1], well known and recognized in aviation industry.At the same time, the 3D vortex lattice method (VLM) has been tested and calibrated, according to the wind tunnel test data of an existing aircraft, the Utva 75.Based on the inviscid CFD concept, the VLM analyses can give results for lift and moment coefficients in their linear domains, and only lift-induced drag, for the entire airplane configuration (without and with flaps and control surfaces deflections).Datcom and VLM gave good agreements for lift and moment curves for the NLA.For better determination of the NLA's polar (lift-drag) curves, a hybrid method has been established, using parasite drag calculus from Datcom, and full-configuration induced drag from VLM.Finally, the viscous CFD calculations have been performed, using the RANS k- SST turbulent model.Although the NLA's cruising Mach number is of the order of 0.15 (compressibility influence is small), for higher accuracy of the results, all CFD settings have been applied as for fully compressible flow analyses.Comparison of the results obtained by Datcom, VLM and hybrid method, and viscous CFD have shown fair agreements for practical engineering purposes.Differences were observed at higher angles of attack, where viscous CFD generally gave slightly larger maximum lift coefficients and lower critical angles of attack, than values obtained by Datcom.

INITIAL AERODYNAMIC CALCULATIONS
At the initial stages of computational analyses, the authors have used a fairly simple CFD model based on a 3D vortex lattice method.Such calculation models [2] inherently neglect viscosity, and thus the boundary layer influence and separation effects are not taken into account.So the effectiveness of the flaps and control surfaces are overestimated at moderate and higher deflection angles, only induced drag can be calculated but parasite drag can not, while the obtained lift and moment curves are confined to their linear domains.

Vortex lattice analyses
In vortex lattice modeling, the airplane is initially segmented into a system of panels.The number of panels applied for the NLA was about 3800 (Fig. 2).Each of the panels (Fig. 3) is represented by a horseshoe vortex of the strength  n (which consists of a bound vortex positioned at the quarter panel chord position, and two semi-infinite trailing vortices), by the control point at 3/4 of the panel chord at its mid section, and the appropriate normal vector n (Fig. 4).

Figure 4. Representation of a single panel
The velocities induced by all horseshoe vortices in a control point of a given panel are calculated using the law of Biot-Savart [2], [3].The summation is performed for all control points, and a set of linear algebraic equations for the calculation of horseshoe vortex strengths is generated.Their  n values are determined by satisfying the boundary condition of "no flow through the surface", i.e. that the total velocity component (free stream plus induced) in the direction of normal vector n must be equal to zero for each of the panels.
Although the applied vortex lattice method is based on the planar presentation of the airplane configuration, the influence of actual mean surface cambers, incidences, dihedral and twist angles, deflections of control surfaces and flaps, must be taken into account.For example, coordinates of NACA 652-415 airfoil [4] applied for the wing had to be assigned, while for fuselage, its side shape had to be defined (see Fig. 5).

Figure 5. Definition of control surfaces and characteristic contours
The vortex strengths are related to the lifting surface circulation, and the pressure differential between the upper and lower surface pressure coefficients is defined as (the convention which was applied in the calculations -see example in Fig. 6).
The pressure differentials are integrated to yield the total forces and moments, and their appropriate coefficients.More details about the performed VLM calculations can be found in [5] and [6].
In this paper only symmetrical airflow cases have been considered, namely, variations of angle of attack without and with different combinations of elevator and flaps deflections (without sideslip or roll, aileron or rudder deflections).After performing initial calibrations for wing and elevator incidence angles (details can be found in [6]), first set of cases analyzed for Utva 75 was the configuration with elevator deflections in range  e = −30 o (up) ÷ +20 o (down), with 10 o steps, without the deflection of flaps.Neglecting the boundary layer effects, VLM inherently overestimates control surface effectiveness.In this case the obtained gradients (slopes) of lift coefficient C L and moment coefficient C M diagrams were good, but shifts, specially of moment curves, were quite large (Fig. 7.).
The use of some generalized corrections, such as suggested in reference [4], failed to give any satisfactory matches with wind tunnel data.Through repeated tests, finally the elevator calibration factor  e values that gave good match were  e = 0.6; 0.75; 0.9; 1.0; 0.82 and 0. Utva 75 has flaps of single slotted type [8], and the same type of flaps has been applied for the NLA.On the other hand, the VLM inherently treats flaps in the same way as the control surfaces -as simple plain flaps.Principal difference between the two types is that slotted flaps rotate about an axis which is below the wing structure.This way, a convergent gap between the slotted flap and the wing structure appears (Fig. 8), and this flap type generates higher lift (or "circulation" in mathematical sense) than the plain flap.Utva's wind tunnel data exist for flap deflections of  = 25 o and  = 45 o (also in combinations with same elevator deflections as before).Uncorrected VLM results for C M were quite good, but C L showed both gradient and shift discrepancies (Fig. 9).After a number of tests, initially for the case  e = 0 o , simultaneous combinations of flaps efficiency  f factor (affecting shift) and circulation calibration Cirw f factor (affecting gradient) have been defined, which gave very good agreements with experiment.The "fine tuning" of C M also required the introduction of the third parameter, the horizontal tail circulation calibration Cirh f due to flaps deflection which, combined with other two, gave precise agreements considering the moment coefficient slopes.Then, tests have also been repeated for all elevator deflections, and  e for  = 25 o and  = 45 o have been successfully determined.
Using interpolation methods, appropriate calibration factors have been recalculated for NLA, taking into account differences of wing areas in flaps domains, and different reference flap deflections ( = 20 o and  = 30 o ).
Details considering this procedure can be found in [9].Calibration diagrams obtained this way for NLA's vortex lattice analyses are shown in Fig. 10.

VLM results and comparisons with Datcom
Beside the VLM analysis, at the initial design stages of the NLA, parallel calculations have also been done by well recognized Datcom method, based on completely different calculation approach (details of this method exceed the scope of the paper, and can be found in [1]).
When placed on the same diagrams, the VLM calculations have overlapped very well (naturally, in linear domains) with Datcom results, verifying the accuracy of predicted calibrations.Within the same VLM calculations, the induced drag polars have been obtained, and their application will be explained in details in the following chapter.

DRAG CALCULATIONS USING HYBRID METHOD
Considering analytical drag calculations, a general formula which describes total airplane drag polar can be written as [10]: and often, the second and third member are grouped as: where "e" represents the Oswald's efficiency factor.For the purpose of here presented calculations, we shall use its initial form.
The first two members on the right side of (1) define the parasite drag, not generated by lift.The min D C represents its minimum value for the entire airplane, while 2 L k C  represents its "position" component, dependant on the angle of attack (and for the ease of calculations, it is formally expressed in terms of the lift coefficient).The third member Di C represents the liftinduced drag, and only this component can be calculated by VLM, because the first two are the consequence of viscous effects.Theory often treats the induced drag of the wing as the only relevant component contributing to it, defined as: where AR represents the wing's aspect ratio ("relative" wing span), u is the correction for non-elliptical wing planform shape, and s defines the influence of fuselage on induced drag.Taking (3) literally, when wing lift coefficient is equal to zero, the induced drag does not exist, which is acceptable for simple preliminary analytical calculations.
On the other hand, this is not quite true.Namely, with flaps and elevator deflected, the total lift of the entire airplane configuration at a certain angle of attack will be zero, but induced drag will not, since the local lift distributions around lifting surfaces, flaps and elevator will generate local lift-induced drags, and thus 0 Di C  .The capability of VLM to calculate the actual induced drag for any given configuration represents a great advantage over analytical methods.Because of that, for aerodynamic calculations of the NLA, a hybrid method approach has been applied.The first two members in (1) have been determined by reliable analytical methods, while the induced drag has been calculated by VLM: The value of NLA's minimum drag coefficient min 0.03737 , without flaps and elevator deflections, has been determined using Datcom.The increase of the minimum drag coefficient due to flaps deflections was calculated by DVL method (referenced and described in [8]), giving min 0.000748 2 VLM 0.03812 0.0145 2 VLM 0.04150 0.0158 Induced drag polars calculated by VLM, and total drag polar curves, determined using ( 5) ÷ ( 7 This way the analytical parabolic polars without elevator deflection, for flaps deflections  = 0 o ,  = 20 o and  = 30 o respectively, have been obtained as: 2 0 0.03737 0.06421 2 20 0.04935 0.06175 2 30 0.05854 0.06248 Polars for other elevator deflections have been obtained using the same principle.Manufacturers of many single-engine light airplanes of the NLA's size and category, with metal structures (which is the case with NLA as well), officially state that maximum lift/drag ratios of their airplanes are of order of ≈12 ÷ 13 (some examples are given in [11]).On the other hand, many pilots report that, during operational use, maximum On the other hand, presented methods have some inherent shortcomings: (A) As already mentioned, the VLM (also sometimes called the "inviscid CFD" method), neglects the viscosity, and can not predict aerodynamic characteristics at high angles of attack, require derivation of control surface calibration functions, etc.
(B) Parabolic polars, such as ( 9) ÷ ( 11), are often used in technical literature and applied in the design practice, but they are not fully accurate.In most cases, airplane minimum drag is not achieved at zero lift, but at some small positive value of the lift coefficient.Also, at high, around-critical angles of attack, drag increases very rapidly, and can not be expressed by simple mathematical formulas, such as parabolic function.In that sense, parabolic polars underestimate drag at very small and very large lift coefficients.
Due to that, at higher design stages, aerodynamic analysis required the application of the CFD model with included viscous effects, taking into account both micro and macro-vorticity influences on the flowfield around the airplane.

Description of the calculation method
For the purpose of the CFD analyses, the CAD model (Fig. 1  In cases with flaps deflected, the accurate 3D geometry of the convergent slot between the flap and the wing structure has also been modeled.On the NLA, flaps extend from the ailerons to the trapezoidal centerplane segments of the wing (Fig. 21).The selection of the meshing method and the applied physical model complexity was based on an optimum compromise dictated by the available hardware resources.After a certain number of test runs, the optimum choice, which gave satisfactory outcomes, was based on the following: (A) Application of a mesh with reasonably small number of elements, in this case of the order of 1.000.000tetrahedral elements for half-model calculations -see Fig. 22, the exact numbers vary slightly for different flaps configurations.(For very large and complex configurations, the fine meshes can exceed tens of millions of elements, but even with very powerful hardware, the CPU time for a single run might be measured in hundreds of hours.)(B) Application of a very sophisticated calculation method, capable of dealing with separation effects at high angles of attack, fully taking into account compressibility effects (for example, as for transonic 3D flow calculations, although compressibility influence at the NLA's cruising Mach number of M = 0.15 is rather small), etc. Calculation of angles of attack around and beyond critical has become possible after the RANS (Reynolds-Averaged Navier-Stokes) k- SST (Shear-Stress Transport) turbulent model [13], [14] has been applied, so this model has been adopted as standard for all calculations.
Since only symmetrical cases have been analyzed, another reasonable simplification was the analysis of half-models, instead of full (also see Fig. 22), while for the visual presentations, the images were generated using mirroring option with respect to the plane of symmetry.The solution convergence had been boosted by the initial optimum reordering of the mesh domain using the Reverse Cuthill-McKee method [15], the application of FMG -the Full Multi-Grid solution initialization at 4 levels [14], [15], and by active solution steering, applying the automatic optimization of Courant number for the achieved solution convergence stage.It had been assumed that the solution for the given angle of attack has converged when the solution monitors for lift, drag, and pitching moment coefficients show no change (constant values) within the last 100 iterations.
By this approach, with the available hardware, the required CPU time for the convergence of the solutions per one angle of attack had been reduced to only 1 ÷ 2 hours (longer time required for higher angles of attack).

Analysis of the results
Results obtained by CFD analyses have provided very useful, both qualitative and quantitative information, about the NLA's aerodynamic characteristics.
Qualitative information have been obtained trough visualization options.Figure 23 shows distribution of eddy viscosity behind the NLA, for the angle of attack  = 14 o .For example, this type of visualization enables quick and efficient verification of flowfield patterns around the tail surfaces.From the isometric view it is obvious that the central domain of the horizontal tail is (inevitably) shaded by fuselage vorticity, but dominant part of it is still in a reasonably "clean" flow, meaning that at the given angle of attack, horizontal tail will preserve its effectiveness and functionality.By making such visual verifications for the whole range of relevant angles of attack, from moderate negative to positive beyond critical (deep stall), it is possible to verify if the horizontal tail has been positioned properly.Figure 24 shows the example of pressure distribution over the NLA, calculated with respect to the ambient pressure at the adopted cruising altitude, taken as zero reference value (this way the underpressure is quantified as negative, and overpressure as positive value).The underpressure (bright yellow) domain over the upper wing surface overflows the fuselage in the domain of the canopy, by which this part of the fuselage additionally contributes to a certain extent to the overall lift.Known is the fact that, in the presence of wing, the fuselage generates more lift than the same fuselage as an isolated body.In the domain of profiled wing tips (which could be treated as very small winglets) the same low pressure extends all the way to the trailing edge, also improving lift in that part of the wing, etc.
Quantitative information obtained by CFD analyses can be categorized as global and detailed.Global information concern the variations of lift, pitching moment and drag coefficients for the entire airplane.Discrepancies are observed at high angles of attack, where CFD gave slightly higher maximum lift coefficient, lower critical angle of attack, and quite abrupt stall, compared with Datcom results.Polar curves also coincide well, except in the domains of very small and very large lift coefficients.CFD has identified large drag divergence in post-stall, and the lift coefficient value C L ≈ 0.2 at which the minimum drag is achieved.
The CFD results obtained for flaps at  = 30 o confirm all general trends of the three aerodynamic coefficient curves obtained by previous methods, but as expected, the differences are larger, and the reasons are quite obvious.In case of analytical methods, the generalized experimental data for single slotted flaps are applied, although for any of flap types, many local variations in actual geometry are possible, depending on the wing airfoil, flaps chord ratio, etc., which affect the accuracy of final results.In case of VLM, flaps calibrations had been derived from an existing design, and characteristics were then recalculated for the NLA's geometry.On the other hand, the CFD calculations have been performed with the real NLA's wing, and actual and very precisely defined 3D slot geometry between the wing and flap.
The results for  = 30 o case show good agreements for lift and moment curve slopes, but slightly lower zero lift angle of attack, and higher maximum lift coefficient.Also, CFD indicates the delayed stall behavior, due to the deflected flaps influence.The logical explanation could be that the convergent gap between the flap and wing accelerates and energizes flow at the separation domain, which is at the junction between the wing and the centerplane, at flaps root (see Fig. 29), thus delaying its propagation of separated flow forward along the chord with the increase of the angle of attack, and postpones the stall.Previous statement is supported by visualization shown in Fig. 29, where dark-red color represents very intensive vorticity, due to massive flow separation (it is interesting to note that, when case in Fig. 29(c) is viewed from bellow, the horizontal tail is almost completely outside of the red zone, contributing to the natural stall recovery tendency).The detailed quantitative information about the aerodynamic characteristics of the NLA have been obtained by grouping several characteristic airplane loft components into zones, or so called "named selections", such as flap, wing, fuselage, etc.One example of their contribution to total drag for a given angle of attack and configuration is shown in Table 1 (also see Fig. 29(a)).Similar tables can be given for lift and moment forces and coefficients as well, for any considered flaps deflection and angle of attack.
The great advantage of such available data is that the partial modifications on the airplane's design (such as adding profiled aerodynamic covers over the wheels, modifying the flaps channel geometry, wing tips shape, etc.), can be monitored both in details, and in the sense of the overall advantages achieved, including the mutual interaction with all other structure members.Such modifications made on the NLA, compared with its initial design, were the introduction of small leading edge strakes on centerplane, replacement of Hoerner wing tips with profiled tips shown in Fig. 21, etc.
Finally, these calculations have provided information which could be treated both as qualitative and quantitative.One of the most important is the spanvise position of the initial flow separation which, specially in case of the training airplanes, should not be in the domain of ailerons.
The initial NLA's wing was of the pure rectangular planform, which always gives initial separation at the wing root, which is very favorable for this category of airplanes.Later, a small wing strake had been considered as an option, as very common feature in case of manly low-wing light airplanes.It relaxes the wing loading at the root, thus reducing the wing-fuselage interference, since less lift is mapped on the fuselage side, and redistributes the wing load slightly outwards (and should be differed from very large wing strakes on jet fighters, which have a completely different role).
The initial strake analyses, considering local lift distribution, have been performed using VLM. Figure 32 shows one of the cases analyzed by this method, for the wing with flaps retracted.For all analyzed angles of attack, the maximum local lift coefficient obtained by VLM for wing with strakes, was at the wing-centerplane junction, as indicated in Fig. 32 (while for the wing without strakes, it was at the root).Although this method could not predict maximum lift coefficient, even at the time of the initial analyses, it had readily been assumed that the initial separation would occur exactly at the wing-centerplane junction, since the general profile of local lift distribution remains practically unchanged until the stall.
This assumption has been fully confirmed after the CFD analyses.Figure 33 shows the NLA with flaps retracted, at  = 15 o , where the initial stall, indicated by dark-red color, occurs exactly at the wing-centerplane junction.This spanvise initial separation position does not change with flaps extended to both of the analyzed angles, the example could be seen in Fig. 29.

CONCLUSION
The paper presents an overview of the aerodynamic methods and procedures applied in the analyses of a light aircraft, at different development stages.
In the early design steps, applied calculations had to be optimized to be time efficient, and enable quick estimates of the aerodynamic characteristics of possible design variations.On the other hand, these calculations had to be as reliable as possible as well, in order to properly predict the right directions in the design development.The analytical and VLM methods had been used simultaneously, first for the purpose of their mutual verification, in the domains where their applicability overlaps.The VLM calculations had been improved considering the flaps and elevator deflections, by introducing the calibration functions, which were derived from wind tunnel tests data of another light aircraft, and then recalculated for the NLA.For several analyzed cases, agreement of the analytical and VLM results was very good, confirming that the posted quality requirements have been satisfied.
At the next level, the two calculation models had been merged into a so called hybrid method for drag calculations.The parasite drag components were calculated by analytical method, since VLM does not take into account the viscosity, while the overall induced drag for the different airplane configurations was calculated by VLM.The results were then superimposed, thus giving a very powerful calculation tool, which utilized the best aspects of the two applied methods.
Finally, the aerodynamic analyses were performed using the CFD method, which took into account viscous, and also compressibility effects (no matter how small).The calculation model utilized mesh with a reasonably small number of elements, and a very sophisticated physical model, based on the RANS k- SST equations for turbulence.By this, and applying the half-model analysis, the CPU time for calculations of different angles of attack had been substantially reduced, while the range of angles of attack had been extended to the post-stall values.
Results obtained by these methods have shown fair agreements for practical engineering purposes, and the requirement in contemporary airplane design, that the calculation tools and methods applied at all different design levels should generally confirm, supplement and fine-tune each-other, has been satisfied.The FAA (Federal Aviation Administration of the USA) has recently posted an advice that, for light airplanes of the conventional design, the wind tunnel tests are not mandatory, if the contemporary calculation methods have been applied at all design stages.Thus the final verification of the presented results for the NLA will be performed through the oncoming flight tests.

Figure 6 .
Figure 6.Pressure differential distribution on the NLA in one of the test cases, with flaps and elevator deflections 7 for the above mentioned  e range, being a non-linear and asymmetrical function with respect to the  e = 0 o .Examples of "raw" VLM results and with elevator effectiveness calibration applied, for elevator deflection  e = −30 o are shown in Fig. 7, where they are compared with wind tunnel test results.

Figure 10 .
Figure 10.Calibration factors for the NLA Figures 11 and 12 show examples of lift curves for different flaps deflections and different elevator deflections, respectively.The next step in the application of VLM in the NLA aerodynamic analyses was determination of linear domains of lift coefficient C L and pitching moment coefficient C M for different elevator and flaps deflections (examples are given in Fig.'s 13 and 14).

Figure 15 .L
Figure 15.The NLA in a configuration and at an angle of attack when total airplane's lift is equal to zero, but induced drag is not An example is given in Fig. 15 where, at fuselage angle of attack  fus = −5.6 o , with flaps deflected to  = 30 o and elevator deflection  e = −30 o , lift coefficient of the NLA is 0 L C  , while 0.0469 0 Di C   .
 = 30 o .According to Douglas Aircraft Company[10], the parameter k from (the classical design, with unswept wings.This way, the general equations for total drag of NLA, with flaps deflections  = 0 o ,  = 20 o and  = 30 o respectively, have been obtained as: ), are shown in Fig's 16 and 17 for two characteristic cases  = 0 o , and  = 30 o .When the total configuration lift is zero, it is obvious that only for the case  = 0 o , and elevator deflections  e = 0 o ,  e = +10 o and  e = +20 o the induced drag coefficient is 0.0 Di C  , while for other cases it can be substantially larger than zero, at 0 L C  .

Figure 16 .
Figure 16.Induced drag (left) and total drag polars (right) for  = 0 o , and different elevator deflections

Figure 19 .
Figure 19.Three characteristic polars of the NLA, obtained by hybrid calculation method Using (9) ÷ (11), / L D C C curves versus L C have been drawn (Fig. 20).Manufacturers of many single-engine light airplanes of the NLA's size and category, with metal structures (which is the case with NLA as well), officially state that maximum lift/drag ratios of their airplanes are of order of ≈12 ÷ 13 (some examples are given in[11]).On the other hand, many pilots report that, during operational use, maximum /

Figure 20 .
Figure 20.The NLA's lift/drag curves for three typical flaps positions,  e = 0 o ) has been used to derive three lofts of the NLA, corresponding to flaps deflections  = 0 o ,  = 20 o and  = 30 o (without elevator deflections).All other undeflected control surface gaps on the airplane have been "sealed", in order to reduce the complexity of the analyses, without any substantial penalties considering the overall accuracy of the results.Example of the NLA's loft used for  = 30 o calculations is shown in Fig. 21.

Figure 21 .
Figure 21.The NLA loft with flaps deflected to 30 o

Figure 22 .
Figure 22.The mesh on the model and in its vicinity, with flaps deflected to 30 o

Figure 23 .
Figure 23.Visualization of eddy viscosity behind the NLA, flaps retracted, angle of attack  = 14 oIt should be emphasized that, for each angle of attack, a sufficient number of different views must be generated to obtain the proper information of the flow patterns around the whole airplane (for example, in Fig.23in side-view, the tip vortex from the left wing visually masks the flow in the tail domain).

Figure 24 .
Figure 24.Visualization of pressure distribution over the NLA loft, flaps retracted, angle of attack  = 14 o

Figure 25 .
Figure 25.Comparisons of lift and moment, case  = 0 o

Figure 29 .
Figure 29.The eddy viscosity around NLA, for case  = 30 o : (a) flow separation and initial stall at flaps roots, at  = 12 o ; (b) developed stall at  = 14 o ; (c) deep stall at  = 18 o The polar obtained by CFD for  = 30 o gives slightly smaller values of drag coefficient than those obtained by hybrid method for lift coefficients C L > 0.6, and larger below this value.These results are encouraging, because flaps are generally not used at lift coefficient values C L < 0.6 considering this category of airplanes, i.e. at higher flight speeds (one of the principles and values stated within the DVL method).Lift/drag ratios are shown in Fig.'s 30 and 31.The CFD calculations are also encouraging in this sense, giving larger maximum C L /C D values approximately by one unit, compared to the corresponding values calculated by previous methods.

Figure 32 .
Figure 32.Local lift coefficient distribution shape, obtained by VLM for the NLA wing, without and with strake