AERODYNAMIC INVESTIGATION OF THE DEFORMABLE MEMBRANE AIRFOIL WITH EXCESS LENGTH

Original scientific paper The deformable membrane airfoils (DMA) are considered in the present paper, which could be potentially used on a new class of modern Micro Aerial Vehicles for their compact and lightweight form. Numerical simulations of low Reynolds number airflow around the DMA with excess length were performed and compared with the results obtained in the experiment. The aerodynamic performance was studied using a fluid-structure interaction between the DMA with excess length and the surrounding fluid flow. Fluid part was solved using Reynolds averaged Navier – Stokes (RANS) equations and applying the resulting pressure distribution to the airfoil. Experiment was performed in the low-turbulence wind tunnel with a model of an airfoil made from a deformable membrane with excess length, fixed at the leading and trailing edge. Measurements were obtained for the predetermined range of excess length ratio (0,025 to 0,150) and angle of incidence (0° to 15°), and showed a reasonable agreement with the numerical results with smallest discrepancies at 1,8 %. Investigation of the excess length shows that it could increase longitudinal stability of the airfoil in the non-oscillating regime, compared to standard rigid airfoils.


Introduction
Deformable airfoils can often be seen in nature: birds, insects, bats, etc. Controlled flexibility of the deformable airfoil enables it, to have the best aerodynamic characteristics at current conditions which can also be applied in aerial vehicles.This can be accomplished with the use of rigid structural elements, which provide stiffness, and flexible membrane, which adapts to the shape that provides desired aerodynamic characteristics, for the chosen flight regime.Wings could be implemented as compact and easily stowable, which due to its small size, would be heavily affected by wind gusts and significant complexities associated to their flight mechanics.To supplement the bulk knowledge in the field of Micro Aerial Vehicle aero mechanics, the aerodynamic investigation of the DMA with excess length membrane is performed with the use of the fluid-structure interaction between the DMA and the airflow.
Pertinent initial research of the flow field around sails was based on analytical methods as shown in [1] and did not consider airflow separation on the upper part of the airfoil.These methods did not produce results in agreement with experimental outcomes.In [2] authors solved the problem of nonviscous flow interactions around a membrane using complex functions, the Kutta-Joukovski theorem, and the small displacement linear theory of solids.Similar procedure was used in [3] and [4], to design the shape of a DMA (sail) with finite difference numerical methods.This approach did not take into account viscosity and airflow separation, and was limited to the calculation of airflow around membrane airfoils, with a small profile camber and angle of attack.A number of measurements on the DMA for small camber at small angle of attack were made in [5], where discrepancies with the results of before mentioned work were found.To obtain improved calculations of the airflow around deformable airfoils, artificial vortices in the upper part of the airfoil were introduced in [6], intending to simulate airflow separation.Authors assumed constant pressure inside the separation vortex.
To properly compute flow separation around DMAs Navier-Stokes equations have to be solved as shown in [7], while to accurately predict the shape of the membrane, the non-linear theory with large displacement must be considered.Due to the numerical complexity, the fluid flow and membrane simulation need to be solved numerically as shown in [8].Because the fluid equations of conservation and the structural equations of state are time dependent, the numerical problem of fluid-structure interaction is solved sequentially, in numerical iteration cycles [9].In this method, the flow field around the membrane airfoil and pressure distribution on a rigid airfoil, is calculated first.Pressure is then applied to the deformable airfoil and a new equilibrium state and geometry is calculated.This approach assumes stationary airflow around the membrane, with small displacement velocity of the airfoil.In [10] authors used indirect numerical iterative methods to calculate the flow around a DMA.They numerically solved the Reynolds Averaged Navier-Stokes (RANS) equations with the use of the - turbulence model.They observed the airflow as twodimensional and the membrane as non-extendable, with no bending resistance.Their research was limited to small camber deformable airfoils, because of flow separation issues and large displacements of the membrane.
In this paper, which will complement the work presented in [11] and [12], a numerical algorithm for the coupled fluid structure interaction is presented, which is made to account for excess length of the membrane that results in large camber.Influence of the excess length is shown with regards to the aerodynamic lift coefficient and pressure distribution.In the first part of the iteration step, a two-dimensional numerical simulation of turbulent airflow around a rigid membrane airfoil was performed, using the - turbulence model.In the second part of the iteration, calculated air pressure around the airfoil from the previous step, was applied on the DMA and solved using large non-linear displacement theory.
For the comparison of the numerical results, experiment is designed and measurements are performed in the wind tunnel.All the predefined ranges of the angle of incidence and the excess length ratio are also tested, for the non-oscillating DMA configurations.

Governing equations of the fluid-structure interaction
The airflow surrounding the DMA provides a 2 way coupled fluid-structure interaction.Because of the highly deformable thin membrane, which is used as an airfoil, the expected change in shape requires the use of a deformable membrane structure with small deformations and large displacements.Since the integral quantities in the flow are of interest, a stationary, non-compressible, viscous airflow is assumed.
For the calculation of the turbulent fluid flow, Reynolds averaged Navier -Stokes equations, using Boussinesq approximation were used.For closure of the system of equations, two equation RNG (Re-Normalisation Group) - model was used [13].
where   is mean fluid velocity,  f is kinematic viscosity of the fluid, and  T is the kinematic turbulent viscosity Constants in (4) ÷ ( 6) are   = 1,0,   = 1,3,  2 T = 1,92,   = 0,09.The former constant of the standard - model  1 T , changes to where  ∞ = 4,38 and  is the renormalization coefficient which is calculated as Sign <> is a sum of squared tensor components by summation indexes  and .Compared to the standard - model, RNG model has additional Eq. ( 4), which attempts to account effects of more than one scale of motion in determining eddy viscosity.Eq. ( 1) ÷ (5) close the system for the solution of the following variables: pressure, turbulent kinetic energy, energy dissipation and velocity vector.These equations are solved iteratively, using the Petrov-Galerkin upwind corrected, finite element method.
Structural strain was calculated using Hook's linear rheological law and incorporated geometric nonlinearities as in [14].The momentum conservation equation for infinitesimal control volume in Lagrange formulation according to [15], is where   is the displacement vector,  is density,  is Young's modulus of elasticity, and  s is the Poisson's ratio.Although the rheological law is linear in the case of large displacements and small deformations, the partial differential equations of momentum conservation are nonlinear and are solved with non-linear iterative numerical methods.
Interactions between fluid and structure are observed on the interface-common border Γ fs of the fluid domain Ω f and the structural domain Ω s as shown on Fig. 1.

Figure 1 Interactions between the fluid and structural domain
On the interface, the magnitudes of the scalar products of the fluid and structural stress tensor and the common border normal are equal, but have opposite direction.
where   f and   s are respectively stress tensor of fluid and solid.The normals   f and   s , have the same orientation, but opposite directions   f = −  s .Despite that the symmetrical stress tensors of fluid and solid are not equal in general, because three additional scalar equations are needed for mathematical problem closure.
For incompressible fluid flow and structures with small deformations, Eq. ( 7) is written as: The stress tensor of the fluid consists of a hydrostatic pressure component  and viscous component.For DMA the pressure part of fluid stress tensor on the fluidstructure interface and the normal components of structural stress tensor are dominant.

Geometrical characteristics of the deformable membrane airfoil with excess length
Since the airfoil in this research was not actively controlled, the flow conditions of the surrounding airflow determined the shape of the DMA through the equilibrium of aerodynamic forces on the membrane with excess length.Therefore the final shape was determined as a final result of a fluid-structure interaction, whereas our starting assumption of the shape was an arc between support points with constant radius of curvature.The initial shape of the DMA attached to rigid support points is presented in Fig. 2. The initial shape of the airfoil is determined by the arc of radius  and membrane with excess length .Airfoil chord  is the distance between the leading and trailing airfoil edges.In the case of a deformable airfoil with excess length, chord  represents the distance between support points.Angle of attack  is defined as the angle between the airfoil chord and the airflow direction.To define the airfoil flexibility,  is defined as a measure of the airfoil excess length ratio of the deformable membrane Flat airfoil has the value of excess length ratio  = 0.The membrane arc angle  is calculated as: Eq. ( 10) is transcendental, and is solved numerically using e.g.Newton's iteration method.Specific camber of the deformable membrane airfoil / is: The basic elements which define the aerodynamic characteristics of rigid airfoils are: angle of attack, airfoil thickness, chord length and camber, whereas the initial shape of the deformable membrane airfoil is defined by: angle of attack , chord length , and excess length ratio .

Numerical modelling of fluid-structure interaction
The displacements of the deformable airfoil are large and non-linear, and direct solution of solid displacement has to be obtained by iterative process.The simulation of the fluid-structure interaction has to be done in sequential iteration steps [16].The fluid flow over a rigid structure is determined first, and the new equilibrium state of the structure is calculated second, calculated from the fluid forces acting on the structure.Fig. 3 describes the numerical algorithm used to solve fluid-structure interactions on a DMA.
As the DMA's initial shape is determined as shown in Fig. 2, the parameters sufficient to determine the initial geometrical shape, are calculated by using: the membrane length , excess length ratio , angle of incidence  and membrane thickness.The finite element mesh is generated based on the initial airfoil shape.The stationary fluid flow is then calculated based on fluid properties and boundary conditions.The pressure distribution around the airfoil from the previous step is applied as a boundary condition for the displacement calculation.
Subsequent forces calculated at the support points were compared, to determine the convergence.If the relative difference between the currently calculated reaction force and the mean value from reaction forces calculated in the previous two iterations is smaller than 2 %, the interaction is complete.For a possible additional iteration, a new mesh is generated around the deformed membrane airfoil, where the new flow conditions are evaluated.Due to large displacements, partial mesh remeshing is necessary for proper flow and structure modelling.Since the solution of the fluid flow is complex and slow, two dimensional modelling is chosen.Finite element mesh used around the DMA, is shown in Fig. 4. Mesh is constructed with two sub domains: a fixed domain in which nodes are fixed within the calculation, and a remeshing domain, where mesh is regenerated by each fluid-structure iteration step.The length of the entire mesh domain is 20 membrane lengths and the height is 4 membrane lengths which are in accordance with the experimental setup.On the upper left part of Fig. 4, the cross-section of the rods on the leading and trailing edge, which support the DMA in the experiment, represents the gap in the mesh.
The fluid domain mesh is composed of approximately 77 000 four node bilinear finite elements where 42 000 of these elements were used in the remeshing domain.At the far distance from the DMA mesh, the characteristic dimensions of the elements are 1/50 of the airfoil length.In proximity to the DMA, the dimensions of the elements are 1/600 of the airfoil length.At the membrane gap and at the upper and lower boarders of the channel, no-slip wall boundary conditions with turbulent wall functions of Van Driest [17] were defined.On the left side of the mesh, the inlet velocity of  0 = 13,2 m/s was defined, and on the right side of the mesh, the outlet pressure  0 = 0 Pa was defined.The comparison of numerical gained aerodynamic forces acting on the DMA, was performed with the measurements of the lifting force on an experimental model in the wind tunnel.The experiment was performed in the low turbulence wind tunnel with test section dimensions 0,355 × 0,407 m, maximum velocity of 25 m/s and turbulence intensity lower than 0,1 %.The location of wind tunnel is approximately 300 m over the sea level and is located at the Faculty of Mechanical engineering in the University of Ljubljana.A deformable membrane airfoil model was mounted in the test section and connected to a measurement chain designed to measure lift coefficient   .
The DMA was mounted in the wind tunnel at the specific excess length ratio and angle of attack.The airflow was then run through the tunnel, which formed the airfoil.When the final shape of the DMA was reached and the membrane was not oscillating, the measurements were performed.Afterwards, the airflow was stopped and the model was tuned for the next measurement configuration -changed excess length ratio or angle of attack.This was done for the previously defined range of the excess length ratios and angles of attack that were simulated numerically.

Deformable membrane airfoil model
Membrane used in the experiment was made out of thin polymer sheet with thickness  = 0,2 mm and Young's modulus  = 3,0 GPa and of a rectangular shape  = 100 mm and span  = 280 mm.The leading and the trailing edge of the membrane were wrapped around a steel rod of diameter   = 4 mm as shown on Fig. 5.This type of fixation does not increase the boundary layer separation area on the upper side of the leading edge.Other two edges of the foil were not attached and were left to move freely.ensure rigidness of the model, two binding rods with 5 mm diameter were added at the bottom of the plate, and one on the top of the plate on the downwind side.Supports did not interfere with the flow regime around the airfoil.
Angle of attack  and excess length ratio  were adjusted using adjustment screw, with the lead of 0,7 mm.The adjustment screw for the setting of the vertical position was 30 mm long, and the one for the horizontal position, was 40 mm long.

Measurement chain
Mean air velocity was measured with a static Pitot tube  = 4 mm, connected to a differential low-pressure transmitter Omega PX655-02DI, with working range of 0 ÷ 500 Pa.For computation of air density, measurements of temperature and relative humidity were performed, using Omega XV-11V transmitter, and measurements of absolute air pressure using Omega PX139-030A4V transducer.Electric signals were digitalized through a 16 channel analogue to digital converter National Instruments PCI-6251 with 16 bit accuracy and refresh frequency of 100 Hz.For the purpose of displaying values of current mean velocity and its progress in time, graphical user interface was created.
Mean air velocity as a function of measured quantities is expressed from the Bernoulli's principle and equation of density calculation for moist air [18] where △  represents the pressure difference on the Pitot tube and  represents the air density, which is expressed with ambient temperature , ambient pressure  0 , specific gas constant for dry air   , specific gas constant for water vapour   , relative humidity  and saturation vapour pressure   , according to [18].For measurements of forces acting on the model, a support system was designed which had 3 support points.On every support point there was one dimensional force sensor with full Wheatstone bridge and carrier amplifier.

Results of numerical simulations and measurements
The chosen pre-determined range for the aerodynamic investigation of the DMA with excess length, was set for a range of airfoil excess length ratios:  = 0,025 to  = 0,150 in steps of 0,025, and angle of attack  = 0, 0 ∘ to  = 15, 0 ∘ in steps of 1, 5 ∘ .In Tab. 1, the membrane arc angle , and specific camber / are calculated, based on the excess length ratio .The flow was simulated around the deformable PVC membrane airfoil with the following parameters: arc length  = 100 mm, thickness 0,2 mm, Young elasticity modulus  = 3,0 GPa, and Poisson number  s = 0,38.
For airflow at standard atmosphere at an altitude of 300 m, the density is  = 1,17 kg/m 3 and the kinematic viscosity is  f = 15,4 × 10 −6 m 2 /s.The air velocity at the inlet was  0 = 13,2 m/s, and the calculated Reynolds number based on the deformable membrane airfoil arc length was Re = 85,7 × 10 3 .
The average time for the calculation of one case (specific  and ), was approximately 49 core hours, on the computer with 12 processors with frequency of 2,93 GHz and 24 GB of RAM.For complete calculation of the previously defined range of excess length ratios and angles of attack, computation took approximately 2500 core hours.
For each measurement case, the pressure difference, temperature, relative humidity, absolute pressure and the aerodynamic forces were measured.From the obtained data, the air density and velocity was calculated.Due to ambient condition changes during the measurements, the parameters were not fixed.The density was 1,155 ± 0,005 kg/m 3 , inlet velocity 13,5 ± 0,2 m/s, and chord Reynolds number 86 × 10 3 ± 2 × 10 3 .

Calculated distribution of the pressure coefficient on the DMA
For the pressure distribution comparison, the pressure coefficient   was used.It is calculated as: where  denotes the static pressure on the surface of the airfoil,  0 and  0 are static pressure and fluid velocity at the far distance from the airfoil respectively.In Fig. 6 the calculated distributions of pressure coefficient on the airfoil, for the excess length ratios  = 0,025 and  = 0,150, are presented.Distributions of the pressure coefficients for the following parameter ranges are presented: angle of attack from 3, 0 ∘ to 15, 0 ∘ at excess length ratio  = 0,025, and angle of attack from 9, 0 ∘ to 15, 0 ∘ at excess length ratio  = 0,150.The vertical axis is inverted, to express the fact that a negative value of   is associated with the positive lift.The lines on the upper side of the graph represent the   on the upper side of the airfoil.On the horizontal axis, the current position on the deformable membrane airfoil  (shown on Fig. 2) is shown, as a percentage of airfoil arc length.The first peak of the pressure on the upper surface occurs after the stagnation point, due to the flow around the leading edge, where the velocity rises.It is larger at smaller excess length ratios because of better flow alignment with the leading edge and therefore higher velocity.At the trailing edge on the lover side of the membrane is another peak, which is a consequence of a local velocity increase, due to the flow reattachment at the trailing edge.At smaller angles of attack ( < 7, 5 ∘ ), for excess length ratio  = 0,025, the   on the upper side of the membrane decreases to a minimum value between 35 % to 55% of the airfoil length.Afterwards it increases steadily, which ensures a stable boundary layer without the airflow separation on the upper surface.This was also observed in the experiment.At higher angles of attack, the separation area with the constant   expands from the trailing edge up to 60% of the chord.Flow on the lower surface separates at the leading edge, and reattaches at the trailing edge.
Only at higher angles of attack  ≥ 9, 0 ∘ at  = 0,150, is the DMA shown to be non-oscillating.Compared to the previous excess length ratio, the leading edge pressure drops are smaller and the significant drops are between 30 ÷ 40 % on the upper edge.Beyond this area, there is a region of strong positive gradient of   , which rises up to a constant value.There airflow separation occurs and area of separated flow expands from the trailing edge up to 60 % of the chord.The flow on the lower side is separated from the leading edge, to the trailing edge, where it reattaches.

Lift coefficients of the DMA
Lift coefficients are used to evaluate the aerodynamic properties of the airfoil.They are calculated as: where  represents the reference area of the membrane airfoil, and   is lift force respectively.The area of a flat airfoil was used as the reference area.In Fig. 7(a) ÷ 7(f), lift coefficients as functions of the excess length ratio and angle of attack are presented.Next to them, the ideal rigid airfoil lift coefficient line is presented, with the slope of 2π/rad, as it is derived with the classical thin airfoil theory.
The calculated curves do not have constant gradients of lift coefficients relating to angle of attack, which would be typical for rigid airfoils at small angle of attack.In the range of angles of attack, where the   slope is greater than the one for the ideal rigid airfoil, the higher lift gradient is a consequence of change of shape, because of the deformability and the excess length.Since membrane was not pre-strained, the change in shape is substantial and greatly influences the lift curves.At higher angles of attack the separation area increases, which leads to smaller lift gradient.At higher  (Fig. 7(f)) the separation area is constant, and leads to constant lift gradient, since the separation occurs at the trailing edge and slowly advances toward the leading edge.
Maximum lift coefficients for all excess length ratios that were simulated and measured are presented in Tab. 2. The maximum coefficient of lift is calculated at the excess length ratio  = 0,050 and angle of attack   max,num = 15, 0 ∘ .The angle of attack for maximum measured lift coefficient is   max,exp = 13, 5 ∘ and excess length ratio  = 0,075.At  > 0,075, the  max decreases with  and reaches the smallest value at  = 0,150.

Position of the centre of pressure
Centre of pressure is the point, where the total sum of the pressure fields act, with no moment about that point.The calculated relative position of the centre of pressure on the DMA (/) (shown on Fig. 2), as a function of angle of attack at a specific excess length ratio, is presented in Fig. 8. Calculated centre of pressure lies between 35,7 % and 53,4 % of the chord length .The maximum change of centre of pressure is at excess length ratio 0,025 which has the greatest range of validated angles of attack.The mean rate of change, of the centre of pressure, is 1,3 % per degree.At excess length ratio of 0,150, the mean rate of change is 0,55 % per degree.This implies that excess length could increase longitudinal stability of the airfoil.The centre of pressure at the minimum angle of attack in simulations is decreasing from 53,4 % for  = 0,025 at  = 0, 0 ∘ , to 42,5 % for  = 0,025 at  = 9, 0 ∘ .The dashed curve on Fig. 8 represents the stability curve and connects the points at the minimum angle of attack (for specific ), from which the membrane was not oscillating.The calculated position of the centre of pressure is shifted downwind, compared to standard rigid airfoils, where the centre of pressure lies between 30 % and 35 % at higher angles of attack.Maximum curvature point of the DMA is near the middle of the airfoil, as opposed to rigid airfoils, where the maximum curvature lies on the upper airfoil surface, near the leading edge.Because of these factors, a DMA has a flight envelope that is shifted towards higher lift coefficients as compared to that of a standard airfoil.

Conclusion
In the present paper, a comprehensive numerical and experimental analysis of the airflow around the deformable membrane airfoil (DMA) was presented with the focus on the influence of the excess length ratio on the aerodynamic characteristics.Numerically obtained lift coefficients were compared with experimental results and showed similar trend with small discrepancies.Excess length of the membrane without pre-tension, proved to have a large influence on the flow separation, which has an impact on the lift and pressure coefficients.
To investigate further into the MAV design with the use of DMA, the airfoil should be transformed into wing and implemented and tested, with different material properties.Due to oscillations, transient modelling of the unfolding of the wing could be considered, along with different types of folding mechanisms, to take the advantage of deformable membranes.Obtained results of the stationary airflow simulations show that DMA allows for lower take-off speeds, but the narrow non-oscillating range of angles of incidence has to be considered.Because of small change in the location of centre of pressure, the size of the flight control surfaces can be decreased and provide better roll stability during stall conditions.
It provides additional transport equations for turbulent kinetic energy  T

Figure 2
Figure 2 Geometry of the initial shape of the deformable membrane airfoil with excess length

Figure 3
Figure 3 Numerical simulation algorithm of fluid-structure interaction used for calculation of a DMA

Figure 4
Figure 4 Finite element mesh of fluid flow calculation domain: entire mesh, remeshing domain and membrane gap 5 Measurements of aerodynamic forces acting on the DMA

Figure 5
Figure 5 Experimental model For measurements of airfoil characteristics, endplates of dimensions 200 × 150 × 5 mm were added perpendicular to the rod.They eliminate the tip effects and enable two-dimensional configuration for the flow over airfoil.Distance between them was 281 mm, leaving a gap of 1 mm between the airfoil and the plate.To

Figure 6 Figure 7
Figure 6 Distribution of the pressure coefficient for excess length ratios =0,025 (left) and =0,150 (right) for specified ranges of angles of attack

Figure 8
Figure 8 Calculated position of the centre of pressure, as a function of angle of attack, at a specific excess length ratio

Table 1
Arc angle and specific camber calculated from the excess length ratio

Table 2
Calculated (num) and measured (exp) maximum lift coefficient  max , at the angle of attack   max at excess length ratio