MATLAB-Final Independent Project

The flight performance is to be calculated at 2000 m above sea level, with pressure and density as per the I.S.A. [3]. for the vehicle's velocity range from 30 to 130 m/s.

Before getting into the 3D(wings') analysis, the behaviour of the airfoil and its data has been taken from an online airfoil simulating platform[4]. For this stage, the 2D lift curve slope of an airfoil, 'a0' is calculated from 0 - 9 degree angles, which has a linear nature. Subsequently, the lift curve slope for the high for a finite wing is calculated, assuming airfoil's efficiency 0.8 as per the following formula[5],

`a = (a0)/(1 + ((57.3*(a0))/(pi*e*AR)))`

This CL(lift coefficient), for the actual wing, can then be obtained by subtracting the selected angle(0 degree alpha) from selected 4 degree, and multipying with 'a'. Hence, a range of CL, with corresponding alphas can be obtained by using 'linspace' command, i.e. (from -15 to +15 with 100 intervals), plotted in Fig(1).

Being in a steady flight, Weight has to be balanced by the Lift Force, therefore here, to retain the flight at low velocity, a higher 'CL' will be required, and with a proportionality of v^2, at high speeds, CL required will be acute. The variation of Cl with velocities can be seen in below chart,

The flight is subsonic and puts a little step in next regime because Mach Number (M) is getting above 0.3, and somekind of inverse nature is observed from the left sided curve, at some point, at low Mach Numbers, the Drag coefficient is low, but this relationship ends after CD rises from 0.1 again, and a proportionality exists.

The below graph indicates the variation in Reynolds Number (Re), on the wing distance. To be clear, near the Root Chord, the Reynolds Number is quite high, while being directed toward the tip section, the value of (Reynolds Number) is decreasing, as if the wing is divided into number of section, at Root, the MAC is more as compared to the tip section.

Secondly, the total drag acting on the aircraft is something that cannot be neglected, as this single parameter is enough to manipulate the desired performance. The Net drag comprises of 1. Parasite drag, acting just due to motion, but not lift, 2. Induced drag when vehicle is in flight condition and 3. Drag arising because of remaining elements, like fuselage, stabilizers, rudders etc.

Here the Parasite drag acting on whole design is calculated in the OpenVSP software[6], using following equations,

For Laminar skin friction coefficient, using Blasius equation,

`Cf = 1.32824 / (Re^(0.5))`,

For Turbulent skin friction coefficient, Power Law Prandtl Medium Re,

`Cf = 0.0315/(Re^(1/7))`

'Re' demotes the Reynold's number, which here is approximately 600000. The obtained averaged Cf, is taken into account with Form Factor (assuming here as 1), with Parasite drag and the ratio of wing's Reference area and Wetter area (Swet).

So, Equivalent Skin Friction, `Cfe = FF * CD0*((Sref)/(Swet))`, the Wetted Surface area for the wing is 49.25 m^2. Additionally the drag from other parts including fuselage is assumed to be 1.2xCDi. The Total drag results by summation of all subcomponents of drag, and its variation is plotted againt alpha and velocity is shown below. The noticeable thing here is that Drag being a function of Lift coefficient (for induced one), with rise proportionately and also direclty proportional to velocity's square with increase too.

Coming the propulsion parameters, as Thrust should be more than Drag for maintaining the flight. Also, the Power required is an important factor to be considered as to keep the fuel consumption efficicient for the given aircraft. The Power, however largely depends on the Drag Coefficient and aircraft weight. The matrix of Drag coefficient obtained previosly is taken, and its recorded values with a function of varying Cl (i.e. CDi) from starting to end, is divided in 100 intervals, and is integrated with the empricial relation of Power, as a function for the given configuration. So Power (in Newtons) is given as,

`P = 1/2*rho*v^3*Sref*CD`

The calculated flight parameters are a prediction of the vehicle's performane with variation in velocity. But, it is helpful in determining the results with input conditions.

% Specifications of MQ-9A UCAV

Wnet = 4760;                 
MAC = 1.2470;         
Sref = 22.48;         
AR = 17.80;
wl = linspace(10,1,100);

% Velocity range in m/s
v = [31:130];            % Max. limit is upto 133 m/s 

% Atmospheric variables at 2000 m under I.S.A.
rho = 0.8217;         % Density in kg/m^3
p = 0.7950            % Pressure in kPa
vu = 1.715e-5;        % Kinematic viscosity in m^2/s
Re = (rho*v*MAC)/vu;  % Reynold's number as per velocity range
A = 332.5;            % Sound speed in m/s
M = v/A;              % Mach number variation


% The 3D Lift curve slope, or wing's is
e = 0.8;                          % efficiency same for airfoil and wing
a0 = 0.0844;                      % lcs of MH114
a = a0/(1 + (57.3*a0/(pi*e*AR))); % lcs of MH114 wing                            
alpha0= -7.25;                    % angle @Cl_MH114 = 0
  
% To find Cl(WING) vs Alpha
alpha = linspace(-15,15,100);
CL=a*(alpha-alpha0);

% PLot of CL(wing) vs Alpha
figure(1);
plot(alpha ,CL ,'linewidth',2, 'color','r');
xlabel('Angle of attack')
ylabel('Wing Lift Coefficient')
title('Lift coefficient vs Alpha')
grid on

% Required Lift coefficient with angle
figure(2);
C_L = 2*Wnet./(rho*Sref*v.^2);
plot(v ,C_L ,'linewidth',2, 'color','r');
xlabel('Velocity');
ylabel('Cl to maintain L=W');
title('Required Lift coefficient with velocity')
grid on

figure(3)
plot(CD, M ,'linewidth',2,'color', 'm')
xlabel('Drag coefficient');
ylabel('Mach Number');
title('Drag variation with Mach Number')
grid on

figure(4)
plot(wl, Re ,'linewidth',2,'color', '[0.4660 0.6740 0.1880]')
xlabel('Distance from root to tip chord');
ylabel('Reynolds Number');
title('Reynolds Number variation on wing')
grid on

% Drag calculation and plotting
CD0 = 0.08866;              % Parasite drag coefficient overall
CDi = CL.^2/(pi*AR*e);      % Induced drag coefficient
CDf = 1.2*CDi               % Fuselage drag estimation
CD = CD0 + CDi + CDf;       % Total drag coefficient
 
figure(5)
plot(alpha ,CD ,'linewidth',2, 'color','r')
xlabel('Angle of attack');
ylabel('Total Drag Coefficient');
title('Drag coefficient vs Alpha')
grid on

figure(6) 
plot(v, CD, 'linewidth', 2,'color','b')
xlabel('Velocity')
ylabel('Drag Coefficient')
title('Drag variation with velocity')

grid on

% To ovrecome Drag force,
D = 0.5*rho*v.^2*Sref.*CD;
T =D;

figure(7)
plot(v, T,'linewidth',2,'color', 'y')
xlabel('Velocity');
ylabel('Thrust required(Newtons)');
title('Thrust vs velocity')
grid on

% Net Power Required, for given configuration
PREQ = power_req(4760 ,0.8217, 22.48)

figure(8)
plot(v, PREQ,'linewidth',2,'color', 'g')
xlabel('Velocity');
ylabel('Power required(Watts)');
title('Required Power vs velocity')
grid on
 

function [Power] = power_req(Wnet,rho, Sref)

CD = linspace(0.105795270359314,0.229896708133328, 100);
v = [31:130]
Power = 0.5*rho*v.^2*Sref.*CD.*v;

end

[1]. https://en.wikipedia.org/wiki/General_Atomics_MQ-9_Reaper

[2]. http://airfoiltools.com/airfoil/details?airfoil=mh114-il

[3]. https://www.engineeringtoolbox.com/international-standard-atmosphere-d_985.html

[4]. http://webfoil.engin.umich.edu/

[5]. http://www.dept.aoe.vt.edu/~lutze/AOE3104/airfoilwings.pdf

[6]. http://openvsp.org/wiki/doku.php?id=parasitedrag