Progress Report, Mesicopter Project

Research and Development Work in the Area of Ultra-Low Reynolds Number Aerodynamics

Introduction:

The operating regime of the mesi-scale helicopter poses certain difficulties for aerodynamic analysis and design. Current sizing and motor parameters result in a rotor tip Reynolds number of approximately 5000. It is unclear to what extent classical airfoil and finite wing analysis and design methods are applicable in this flow regime. The highly viscous nature of the flow field, large increases in the boundary layer thickness, and the potential for large regions of separated flow, all create the potential for large discrepancies in performance from what might be expected from experience at higher Reynolds numbers.

Very little experimental or computation work exists to date for aerodynamic lifting surfaces operating at such low Reynolds numbers. Much of the work completed thus far under the current program has gone towards a comparison of the analysis tools currently available, investigation of the unique characteristics of two-dimensional airfoils operating in this regime, and the assessment of the potential performance benefits of detailed section and planform design. The last point is of particular interest for it forms a strong motive and basis for future research and provides insight into the complexity and depth of the available design space.

Current Tools:

Current computational analysis tools fit into two categories, full viscous flow field solvers working with some formulation of the Navier-Stokes equations, and methods which divide the flow field into an outer inviscid flow region and a viscous boundary layer.

Two different Navier-Stokes codes have been examined for application to this work. The first is the FLO103 code [1], a RANS solver for the full Navier-Stokes equations. This solver was not intended for low Mach number and ultra low Reynolds number analysis and has had very limited success in this current application. Fully laminar converged solutions have been obtained for Mach Numbers as low as 0.3 and Reynolds numbers as low as 20,000, but this is still significantly out of the current operating regime. Further improvements in convergence may be possible through careful grid refinement, but this tool appears to be only marginally applicable to this problem. It is important to note that this is not a shortcoming of this particular code, but of this type of flow solver in general.

The second Navier-Stokes solver considered is the INS2D code, developed by Stuart Rogers at NASA Ames [2]. This code solves an incompressible formulation of the Navier-Stokes equations. It is much better suited to analysis at ultra-low Reynolds numbers and has been used for the majority of two-dimensional section analysis to date.

The inner/outer flow field methods available solve the outer flow field using either a potential flow formulation via a classical panel method, or by solving the Euler equations for inviscid, compressible flow. The viscous boundary layer solution is generally some variant of an integral formulation of the Prantdl boundary layer equations.

One inner/outer flow field code that has been implemented, with limited success, is the MSES program developed by Mark Drela at MIT [3]. This is a two-dimensional Euler solver coupled with an integral boundary layer formulation. It appears to give reasonable drag predictions over a narrow range of angle of attack, but the limitations of the boundary layer formulation cause the solution to diverge if substantial regions of separated flow exist. Unfortunately, this is often the case with these airfoils at even moderate angles of attack.

The much simpler and faster panel methods can be unsuited for analyzing very thin sections, below 1% to 2% thick. Depending on the implementation of the surface boundary conditions and the types and placement of singularity elements, the resulting linear system may become nearly singular. The current section thickness is approximately 2% thick due to manufacturing constraints. This does allow for limited analysis using any of the commonly available airfoil design and analysis codes such as Xfoil [4] or the Eppler code [5]. These codes face the same severe angle of attack limitations as MSES, since all share similar boundary layer formulations. For incompressible flows with adequate section thickness, panel methods are faster and appear to be considerably more robust. For very thin sections or higher Mach numbers, MSES would be more applicable.

Comparison of MSES and INS2D Solutions:

MSES and INS2D fully laminar, steady-state solutions have been compared for two airfoils at a Re 1000. The airfoils sections are the NACA 4402 and 4404 profiles. The resulting lift curve and drag polar plots are provided in Figure (1) and Figure (2). Each data set is complete up to the point at which a given code failed to converge.

Several trends are apparent from the figures. The most obvious is that MSES fails to converge significantly earlier than the INS2D code, reaching only 50% of the INS2D predicted C_lmax . The angle of attack at which it fails roughly corresponds to the development of substantial stable eddy regions along the upper surface in the INS2D solution. The coupled Euler and integral boundary layer formulation has only limited capability to deal with laminar separation bubbles and small regions of separation. The INS2D solution continues to converge past C_lmax until the behavior of the flow field becomes unsteady at high angles of attack. There is the capability to run the INS2D code for an unsteady, time varying solution. This would likely be incorporated into future work.

The other visible trends are consistent shifts in both the lift curve plots and the drag polars. The is a 0.05 constant shift in the lift coefficient corresponding to a 15% to 25% variation between the two solutions. The constant shift in drag for a fixed lift is approximately 50 counts, corresponding to a 5% shift in Cd_o .

This variation can be partially attributed to the governing equations of the integral boundary layer formulation. The Prandtl boundary layer equations assumes pressure is constant along the perpendicular from the surface to the edge of the boundary layer. This can introduces a significant error for ultra-low Reynolds numbers. Under these conditions the pressure is not necessarily constant through the boundary layer. The variation in pressure coefficient through the boundary layer is shown in Figure (3) for the upper surface 50% chord location of a NACA 4402. The section is operating at one degree angle of attack and Re 1000. There is variation of 18% in the pressure coefficient between the surface of the airfoil and the outer edge of the boundary layer. The outer edge of the boundary layer is defined as the location of maximum velocity along the normal to the surface. This variation will be related to the thickness of the boundary layer, with less variation forward of this station and more aft, but it is clearly non-zero. Neglecting this variation introduces an error in both lift and drag calculations.

Considering the added complexity and computational effort required for the Navier-Stokes solution, MSES performs very well over its much smaller range of convergence. It is the inability to converge at higher angle of attack that firmly establishes the need for a higher fidelity solver, capable of dealing with the large regions of separation that form. Reaching only 50 % of maximum lift and 65% of the maximum lift to drag ratio is inadequate. Although the lift and drag of the simpler method agrees well with the Navier-Stokes solution at low lift coefficients, the operational requirements of this very small flight vehicle necessitates maximizing the performance of the airfoil section.

General Performance Characteristics of Airfoils at Ultra-Low Reynolds Numbers:

The primary goal of this work has been to establish the effects of section camber, thickness, and operational Reynolds number on two-dimensional airfoils. All sections are members of the NACA 4-digit family. All analyses utilized the INS2D flow solver. All solutions are steady state and fully laminar. The steady state assumption is supported by the adequate convergence histories of the flow solutions. The fully laminar assumption is reasonable and in fact necessary at the ultra-low Reynolds numbers, 20,000 and lower, being considered. Imposing a transition location would likely be an artificial mechanism that could not be physically realized. Even if sufficiently large artificial trips could be implemented without causing complete laminar separation, it is likely that the transition length would exceed the chord of the airfoil.

Basic Camber and Thickness Effects:

The initial analyses considered the general effects of camber and thickness at a Re 1000. Three sections were considered: NACA 0002, 4402, and 4404. The resulting drag polars are provided in Figure (4). Prior to this analysis is was unclear if the overall effect of camber would be the same at very low Reynolds numbers. It was thought that the decambering effect of the considerably thicker boundary layers might effectively wash out the camber completely. The INS2D results for the 0002 and the 4402 demonstrate that camber is still a very effective parameter. The 4% cambered section exhibits the expected upward shift, providing lower drag at higher lift coefficients, but the addition of camber also significantly increases the predicted C_lmax of the section from 0.55 to 0.74. These two effects combine to increase the maximum lift to drag ratio from 4.5 to 5.5. The C_l for maximum L/D also increases from 0.54 to 0.71. The effect of section thickness on drag is also visible in Figure (4). The increase in thickness from 2% to 4% results in a 15% drag penalty around maximum L/D and a similar reduction in the magnitude of the maximum L/D.

Reynolds Number and Additional Thickness Effects:

Using the 4402 and 4404 sections, additional runs were completed to explore the effect of Reynolds number variations and to gain further insight into thickness effects.

The lift curves and drag polars for the two sections over a range of Reynolds numbers are provided in Figure (5) and Figure (6). The lift curves show the general trend of increasing C_lmax with decreasing Reynolds number. This may seem counter intuitive, but the flow field appears to become much less susceptible to the onset of unsteady separated flow with decreasing Reynolds number. Lower Reynolds number implies a reduction in the ratio of inertial to viscous effects in the flow field. The increase in the viscous effects can be seen in Figure (6) as the substantial increase in drag with decreasing Reynolds number. The effect of the reduction in the inertial effects is seen in the increase in C_lmax. As Reynolds number is decreased, large stable eddies develop on the upper surface of the airfoil. The inertially damped flow field is less likely to separate catastrophically and exhibit unsteady behavior.

This lift behavior is opposite, and the drag behavior is consistent, with what is seen as normal behavior for airfoils as Reynolds number is decreased from high values down below 100,000. In this regime, the increase in the viscous effects, the increase in boundary layer thickness, etc., must be more destabilizing then the reduction in inertial effects is stabilizing. Then at some point above the range of Re considered here, there is a crossover, and the reduction in flow inertia begins to win out, delaying the onset of unsteady separation and increasing the maximum section lift coefficient.

The trends in C_lmax and maximum L/D for the 4402 and 4404 are provided in Figure (7). As the Reynolds number approaches zero, the C_lmax increases from rather low values towards values approaching one. The C_lmax recovers to a value quite similar to those generally attainable for Reynolds numbers of order 100,000 and greater. The increase in lift does not compensate for the rapid increase in drag when considering L/D , but even as the Re approaches zero maximum L/D stays well above one, approaching 2.5 in the limit. Similar trends have been observed in limited experimental data utilizing a 14% thick airfoil [6]. This experiment found maximum lift coefficients in the range of 0.8 to 0.5 in the Reynolds number range of 100 to 1000. The reported drag is substantially higher, but this is likely due to the large increase in section thickness.

The trends depicted in Figure (7) may appear erroneous in that the maximum L/D appears to asymptotically approach values between 10 and 14. Typical sections at much higher Reynolds numbers exhibit values of order 100. The reason for the plotted behavior is that the thin, relatively sharp-nosed sections, become less efficient as the Reynolds number is increased. At the Reynolds numbers being considered, small magnitude changes correspond to substantial changes in the behavior of the flow field. In order to see the expected trend, a series of different sections, each optimized for a small range of operation, would have to be superimposed.

Also visible in Figure (5) and Figure (6) is a lateral and upward shift of the lift curve and drag polar with decreasing section thickness. This effect is similar to that associated with an increase in camber. The increased thickness of the 4404 section appears to increase the decambering effect of the boundary layer.

Parametric Study of Camberline Variations:

The previous results have shown that section geometry variations have a considerable effect on section performance. Previous analyses considered only one camber location and magnitude. In an attempt to further improve section performance and explore the value of a detailed design effort, a nine point parametric study of section performance has been completed. The varied parameters are the camber magnitude and the chordwise location of maximum camber. All sections are 2% thick with an NACA 4-digit camberline and thickness distribution. All sections were analyzed at Re 12,000. The test points and a summary of the results are shown in Figure (8) and Tables 1, and 2.

	Max.	Camber
Max. Camber Location	2%	4%	6%
30%	14.67	13.47	10.35
50%	14.78	14.09	11.87
70%	15.40	16.20	14.78

Table 1 Maximum L/D vs. Camber and Camber Location

	Max.	Camber
Max. Camber Location	2%	4%	6%
30%	0.46	0.45	0.38
50%	0.47	0.51	0.47
70%	0.49	0.64	0.71

Table 2 CL for Maximum L/D

The performance of the sections varies considerably from point to point, supporting the premise that detailed section geometry design is a worthwhile endeavor. Somewhere within the design space exists at least one optimal design. The substantial variations in performance across the test matrix imply that an optimal design may offer significant performance gains over the limited number of sections considered thus far.

The only strong trend visible in the results is a performance advantage in conjunction with an aft shift in the maximum camber location. For this limited test matrix, the 4702 section exhibits the best performance in terms of high maximum L/D occurring at a high lift coefficient. This section will likely be implemented in the next iteration of rotor design, replacing the current 4402. Figure (9) incorporates the 4702 drag polars with those of the 4402 and 4404 for several Reynolds numbers. The much higher C_lmax of the section results in considerable gains in L/D, but the stall characteristics of the section at the higher Reynolds numbers appear to be more abrupt and less forgiving. In application, it would be prudent to back away somewhat from the maximum L/D in order to provide some stall margin.

Areas for Future Work

The two-dimensional work to date has established that he the ultra-low Reynolds number regime offers a fruitful area for research and an operationally feasible regime for micro air vehicles, more specifically the mesicopter. The lack of quantitative historical research into the aerodynamics of this regime, coupled with the preliminary results of the research associated with this program, hold promise for a significant and relevant contribution to the state of the art. To bring this flight vehicle to fruition and to adequately develop our understanding of the aerodynamics, there are still many areas for continued effort.

Two-Dimensional Analysis

The INS2D analyses completed thus far indicate that as the airfoil sections approach C_lmax, relatively large stable eddies form along significant portions of the upper surface. These regions appear to grow in strength and size with angle of attack until the flow becomes unsteady and the steady-state analysis fails. This gradual growth process is seen in Figure (1) as a soft stall behavior, which diminishes with increasing Reynolds number. The character of the stall region is of particular interest in this case because the location for maximum L/D is very close to C_lmax. If a soft stall is present, it allows one to maximize the performance of the airfoil with a larger margin of safety and confidence. A hard stall behavior forces the design operating point to fall below the maximum performance point in order to have some reasonable operating margin.

In order to provide some validation of the analyses and to expand our current understanding, two areas of two-dimensional work would be useful, one analytical and one experimental. The analytical work would extend the two-dimensional analyses with a fully unsteady, time accurate method. The currently used code, INS2D, has most of this capability and would require only minor modifications. This would allow for analysis beyond the current analytical limitation and well into the post stall region. It would also help to validate the existence and stability of any large stable separation regions prior to reaching C_lmax.

Two-Dimensional Experimental Testing:

Following further computational work, it would be beneficial and appropriate to develop some form of experimental validation. A two-dimensional airfoil test would provide both quantitative validation of predicted lift and drag forces, as well as an improved qualitative understanding of the flow using visualization techniques. The primary difficulty is in locating and obtaining access to suitable test apparatus. At the Reynolds numbers of interest some fluid other than air must likely be used in order to achieve reasonable physical model size for fabrication and instrumentation, and in order to bring the dynamic pressure and the forces up to a level where accurate measurements are possible. Some form of oil is the most likely candidate. For a Reynolds number of 5000 per foot, flow velocities on the order of one foot per second would be necessary, generating dynamic pressure on the order of one pound per square foot.

Instrumentation would depend on the capabilities of the facility, but ideally would utilize surface pressure transducers for the integration of lift and wake survey capabilities for drag computation. Two-dimensional airfoil section testing using force balances can be problematic from the standpoint of integrating the model with the balance. The model must be completely isolated from the channel walls, introducing gaps, or a portion of the wall must be made metric with the balance. In either case there is considerable added complexity and the effects of the wall boundary layer must be accounted for. An advantage of a balance would be in the ability the measure drag forces beyond C_lmax, where accurate wake integration would become problematic. In spite of any difficulties, experimental validation of the computational analyses would be of high value in refining and improving the accuracy of current analysis and design tools.

Two-Dimensional Geometry Optimization:

Analysis of a simple 9-point test matrix of airfoil geometry variations has shown that considerable variations in performance are possible. These results are displayed in Figure (8) with the variation of maximum L/D provided in Table (1). The results do show certain trends, but the complexity of the design space, which incorporates highly non-linear flow behavior, makes the selection of an optimal configuration non-intuitive.

Preliminary work has begun on implementing numerical optimization of the airfoil geometry to maximize section performance. A numerical optimization code will be coupled with an automatic grid generator to create candidate designs. The analysis of candidate designs will be done with the INS2D solver. The airfoil will be represented by a 2% thick NACA 4-digit thickness distribution fitted over a spline-generated camberline. The free design parameters will be the location of a finite number of interior knots. These knots, along with the fixed endpoints at the leading and trailing edges, uniquely determine the camberline shape. Each knot will have only one degree of freedom perpendicular to the chordline and will be constrained to move only a small fraction of the chord length in either direction. This limits the complexity of the problem and restricts the optimizer from attempting highly non-feasible designs.

The angle of attack could also be incorporated as an additional design parameter or iteration over the angle of attack could occur within each design iteration. Including it as a design variable is tremendously more efficient, but may pose difficulties if the optimizer requests an angle of attack for which the analysis fails to converge. Continuing on with the optimization after this occurs could be problematic. Excluding angle of attack as a design variable for the optimizer is safer, but much more costly. The latter option will likely be implemented first to assure a result, but the problem handling non-feasible candidate designs is not trivial and is worth investigation.

Three-Dimensional Analysis:

The transition from two-dimensional sectional properties to the efficient design of a three-dimensional planform introduces many factors into the design analysis which require further exploration. One issue is the stability of the large separation regions, predicted by two-dimensional analysis, in the presence of the three-dimensional flow field of a finite wing. It is possible that additional factors, such as cross flow and finite-wing wake effects, could destabilize these regions and substantially reduce the maximum lift and overall performance. As in the airfoil section analyses, both steady and unsteady calculations will be performed. This problem will be investigated using the INS3D code, a three-dimensional variant of the currently used INS2D code.

For the mesicopter rotor blade application, blade/wake interaction effects may be of great significance. The low section Reynolds number results in a viscous wake thickness on the order of the half-chord. The impingement of the wake of the leading blade upon the trailing blade would certainly have an effect on performance. Investigation of this area would require a time-accurate unsteady rotor code. Such codes do exist and are normally used for large-scale rotor applications, but could likely be adapted for this problem with minimal modifications.

REFERERNCES:

1. Martinelli, L., Jameson, A., "Validation of a Multigrid Method for the Reynolds Averaged Equations," AIAA Paper 88-0414, AIAA 26^th Aerospace Sciences Meeting, Reno, NV, January 1988.
2. Rogers, S. E., Kwak, D., "An Upwind Differencing Scheme for the Time Accurate Incompressible Navier-Stokes Equations," AIAA Paper 88-2583, June 1988.
3. Drela, M., Giles, M. B., "ISES - A Two-Dimensional Viscous Aerodynamic Design and Analysis Code," AIAA Paper 87-0424, AIAA 25th Aerospace Sciences Meeting, Reno, NV, January 1987.
4. Drela, M., "XFOIL – An Analysis and Design Program for Low Reynolds Number Airfoils," Low Reynolds number aerodynamics; Proceedings of the Conference, Notre Dame, IN, June 1989, Springer- Verlag, Berlin and New York, 1989, p. 1-12.
5. Eppler, R., Somers, D. M., "A Computer Program for the Design and Analysis of Low-Speed Airfoils," NASA-TM-80210, August 1980.
6. Thom, A., Swart, P., "The Forces on an Aerofoil at Very Low Speeds," Journal of the Royal Aeronautical Society, Vol. 44, pp. 761-770, 1940.