Lift Distributions on Low Aspect Ratio Wings at
Low Reynolds Numbers
by
Sagar Sanjeev Sathaye
A Thesis
Submitted to the Faculty of
WORCESTER POLYTECHNIC INSTITUTE
in partial fulfillment of the requirements for the
Degree of Master of Science
in
Mechanical Engineering
by
_____________________________
May 2004
Approved:
________________________________________________________
Professor David J Olinger, Thesis Advisor
________________________________________________________
Professor Hamid Johari, Committee Member
________________________________________________________
Professor William Durgin, Committee Member
________________________________________________________
Professor John Sullivan, Graduate Committee Representative
Abstract
The aerodynamic performance of low aspect ratio wings at low Reynolds
numbers applicable to micro air vehicle design was studied in this thesis. There is an
overall lack of data for this low Reynolds number range, particularly concerning details
of local flow behavior along the span. Experiments were conducted to measure the local
pressure distributions on a wing at various spanwise locations in a Reynolds number
range 3×104 YL1(i)& YL1(i+1)>YU1(i+1))|(YU1(i)==YL1(i))|(YU1(i)
YL1(i+1)
intsec(j)=X(i);
j=j+1;
end
end
r=j;
intsec(r)=X(n);
for j=1:r-1
k=j;
X1=[intsec(j):h:intsec(j+1)];
YU2=spline(xu,cpu,X1);
YL2=spline(xl,cpl,X1);
area1(k)=trapz(X1,YU2);%area integral- upper surface
area2(k)=trapz(X1,YL2);%area integral- lower surface
area3(k)=(area1(k)-area2(k));
j=j+1;
end
cl=sum(area3)%Lift Coefficient
110
%Program to plot lift coefficient against percent wing span for
measured data and elliptic distribution%
clear all;
clc;
clf;
h=0.001;% X co-ordinate accuracy
b1=(8*25.4)/1000;% Wing Span
b2=1;%normalised wing span
c=(8*25.4)/1000;% chord length
S=(64*((25.4/1000)^2));% Wing area
rho=1.2;%air density
x=[1,0.9375,0.8125,0.6875,0.5625,0.4375,0.3125,0.1875,0];% Z/(b/2) for the wing
cl= input('Enter the cl(coefficient of lift) matrix')% Cl Data
v= input('Enter the free stream velocity')% free stream velocity
X=[0:h:1];% Spline for x co ordinate
Y=spline(x,cl,X);% Spline for upper surface
plot(X,Y)% To plot cl vs Z/(b/2) curve
hold on;
area=trapz(X,Y);% To calculate the area under the curve, i.e. the total lift
ymax= (4*area/pi);%to calculate the maximum lift at the root
n=(1/h)+1;
b(1)=0;
y(1)=ymax;%Minor axis of the elliptical lift distribution
j=(1/h);
for i=2:n% To plot the elliptical lift distribution for the same amount of lift
b(i)=(h*(i-1));
y(i)=sqrt((ymax^2)*(1-(b(i)^2)));
i=i+1;
end
for k=1:j% For getting the ratios of the two lift coefficients
ch(k)=y(k)/Y(k);
k=k+1;
end
xc(1)=0;
for r=2:j % x co-ordinate for the plot of ratios
xc(r)=(r-1)*h;
r=r+1;
end
plot(b,y)
hold off;
figure(2)
plot(xc,ch)
CL=area
111
%Program to Calculate the Fourier Coefficients and span
efficiency factor based on the measured lift distribution and
elliptic distribution%
clear all;
clc;
clf;
format short;
h=0.001;% X co-ordinate accuracy
x=[1,0.9375,0.8125,0.6875,0.5625,0.4375,0.3125,0.1875,0];% Z/(b/2) for the wing
cl= input('Enter the cl(coefficient of lift) matrix');% Cl Data
%v= input('Enter the free stream velocity')% free stream velocity
X=[0:h:1];% Spline for x co ordinate
Y=spline(x,cl,X);% Spline for upper surface
figure(1)
plot(X,Y)% To plot cl vs Z/(b/2) curve
hold on;
area=trapz(X,Y);% To calculate the area under the curve, i.e. the total lift
ymax= (4*area/pi);%to calculate the maximum lift at the root
CL=area
n=(1/h)+1;
j=(1/h);
for i=1:n% To plot the elliptical lift distribution for the same amount of lift
b(i)=(h*(i-1));
y(i)=sqrt((ymax^2)*(1-(b(i)^2)));
b1(i)=-b(i);
i=i+1;
end
plot(b,y)
hold off;
for k=1:j% For getting the ratios of the two lift coefficients
ch(k)=y(k)/Y(k);
xc(k)=(k-1)*h;
k=k+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nn=40;
for i=1:nn
theta(i)=((i)/nn)*(pi/2);
end
for j=1:nn
xx(j)=cos(theta(j));
end
for ii=1:nn
for jj=1:n
if abs(b(jj)-xx(ii))<0.0005
kk(ii)=jj;
end
end
end
for i=1:length(kk)
m=kk(i);
gamma(i)=(1/4)*Y(m);
end
for p=1:length(kk)
for s=1:length(kk)
R(p,s)=sin((2*s-1)*theta(p));
end
end
Aact=R\gamma';
Aact(1:3);
for i=1:length(kk)
deltas1(i)=(2*i-1)*(((Aact(i)/Aact(1)))^2);
112
end
delta1=sum(deltas1);
eact=(delta1)^(-1)
%%%%%%%%%For elliptic distribution
for i=1:length(kk)
m=kk(i);
gamma2(i)=(1/4)*y(m);
end
Aellip=R\gamma2';
Aellip(1:3);
for i=1:length(kk)
deltas2(i)=(2*i-1)*(((Aellip(i)/Aellip(1)))^2);
end
delta2=sum(deltas2);
eellip=(delta2)^(-1)
113
APPENDIX C
114
For α = 15o
An Re = 30218 Re = 35966 Re = 43615 Re = 49345 Re = 84122
A1 0.1437 0.1363 0.1413 0.1293 0.1292
A2 0.0323 0.0355 0.0358 0.0273 0.0339
A3 0.0232 0.0163 0.0194 0.0093 0.0165
A4 0.0061 0.0091 0.0086 0.0025 0.0018
A5 -0.0054 -0.0067 -0.0059 -0.005 -0.005
A6 -0.0082 -0.009 -0.011 -0.0054 -0.0074
A7 -0.0116 -0.0085 -0.0088 -0.0071 -0.0079
A8 -0.0008 -0.002 -0.0009 -0.0002 -0.0005
A9 0.0028 0.0027 0.0019 0.0015 0.0014
A10 -0.0013 -0.0008 -0.0008 -0.0013 -0.0008
A11 -0.001 -0.0013 -0.0012 -0.0001 -0.0008
A12 -0.0011 -0.0008 -0.0008 -0.001 -0.0007
A13 0.0004 0.0004 0.0002 0.0003 0.0002
A14 -0.0001 -0.0002 0 -0.0001 0
A15 0 0.0001 0 0.0001 0
A16 -0.0002 -0.0002 -0.0001 -0.0002 -0.0001
A17 0.0001 0.0001 0 0.0002 0
A18 -0.0002 -0.0002 -0.0001 -0.0003 -0.0001
A19 0.0002 0.0001 0.0001 0.0002 0.0001
A20 -0.0002 -0.0002 -0.0001 -0.0002 -0.0001
A21 0 0 0 0.0001 0
A22 -0.0003 -0.0003 -0.0002 -0.0002 -0.0002
A23 0 0 -0.0001 0 -0.0001
A24 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001
A25 0.0001 0.0001 0 0.0001 0
A26 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001
A27 0.0001 0.0001 0 0.0001 0
A28 0 -0.0001 0 -0.0001 0
A29 0.0002 0.0002 0.0001 0.0002 0.0001
A30 0 0 0.0001 0 0
A31 0.0002 0.0001 0.0001 0.0002 0.0001
A32 0 0 0 -0.0001 0
A33 0.0001 0.0001 0 0.0001 0
A34 -0.0001 -0.0001 0 -0.0001 0
A35 0.0001 0.0001 0.0001 0.0001 0.0001
A36 0 0 0 -0.0001 0
A37 0.0001 0.0001 0.0001 0.0001 0.0001
A38 0 0 0 -0.0001 0
A39 0.0001 0.0001 0 0.0001 0
A40 -0.0001 -0.0001 0 -0.0001 0
115
For α = 6o
An Re = 30218 Re = 35966 Re = 43615 Re = 49345 Re = 84122
A1 0.0586 0.0583 0.0582 0.0567 0.0569
A2 0.0042 0.0018 0.0036 0.0036 0.0098
A3 0.0023 0.0024 0.003 0.0045 0.0069
A4 0.002 -0.0014 0.0006 0.0018 0.0023
A5 -0.0024 -0.0012 -0.0016 -0.0012 -0.0009
A6 -0.0034 -0.0009 -0.0019 -0.0026 -0.0031
A7 0.0002 -0.0011 -0.002 -0.0024 -0.0044
A8 -0.0011 -0.0013 -0.0001 -0.0005 -0.0004
A9 0 0.0006 0.0003 0.0006 0.0011
A10 0.0006 0.0006 -0.0002 0 -0.0003
A11 -0.001 -0.0011 -0.0002 -0.0006 -0.0006
A12 0.0003 0.0003 -0.0002 -0.0001 -0.0003
A13 -0.0001 0 0 0 0.0001
A14 0 0 0 0 0
A15 0 -0.0001 0 0 0
A16 0.0001 0.0001 0 0 0
A17 -0.0001 -0.0001 0 -0.0001 0
A18 0.0001 0.0001 0 0 -0.0001
A19 -0.0001 -0.0001 0 0 0
A20 0.0001 0.0001 0 0 0
A21 -0.0001 -0.0001 0 0 0
A22 0 0 0 0 -0.0001
A23 -0.0001 -0.0001 0 0 0
A24 0 0 0 0 0
A25 -0.0001 -0.0001 0 0 0
A26 0.0001 0.0001 0 0 0
A27 -0.0001 -0.0001 0 0 0
A28 0.0001 0.0001 0 0 0
A29 0 0 0 0 0
A30 0.0001 0.0001 0 0 0
A31 0 0 0 0 0
A32 0.0001 0.0001 0 0 0
A33 0 0 0 0 0
A34 0 0 0 0 0
A35 0 0 0 0 0
A36 0 0 0 0 0
A37 0 0 0 0 0
A38 0.0001 0.0001 0 0 0
A39 0 0 0 0 0
A40 0 0 0 0 0
116
APPENDIX D
117
Sample error calculations:
The error calculation in the pressure coefficients is based on RSS (Root sum of squares)
type uncertainty. A sample calculation and a sample table for errors in pressure
coefficient for Re = 30218 and α = 15o is shown here.
Error in Pressure coefficient:
For a specific pressure coefficient, Cp and a dynamic pressure, q we have
∆P
Cp =
q
where ∆P = P − Ps
where P = pressure being measured (upper or lower surface pressure)
Ps = static pressure
q = dynamic pressure
2 2
∂Cp ∂Cp
δ Cp = ⋅ δ ( ∆P ) + ⋅ δ (q)
∂ (∆P) ∂ (q)
Here δ(∆P) = δ(q) = 0.0005, which is the accuracy of pressure transducer. The table on
following page demonstrates a sample set of error calculated in pressure coefficients.
So for a particular pressure measurement we can obtain the error in pressure coefficient.
The error in pressure coefficient calculation introduces error in the local lift coefficient
calculation. This was based on the RSS type uncertainty as well. The error in cl is then
given as:
δ cl = ∑ (δ C pu − δ C pl )
So by knowing the error in Cp we can obtain the error in local lift coefficient. The code
for this calculation was done in MATLAB.
118
Error in pressure coefficients for Re = 30218 and α = 15o:
Error calculation in pressure coefficients
Upper Surface Lower Surface
Port # Cp,u Error in Cpu Port # Cp,l Error in Cpl
1 -0.08 0.042 1 -0.08 0.042
2 -0.75 0.052 13 0.58 0.048
3 -0.75 0.052 14 0.25 0.043
4 -0.75 0.052 15 0.17 0.042
5 -0.75 0.052 16 0.08 0.042
6 -0.75 0.052 17 0.08 0.042
7 -0.50 0.047 18 0.08 0.042
8 -0.17 0.042 19 0.08 0.042
9 -0.17 0.042 20 0.08 0.042
10 -0.08 0.042 21 0.08 0.042
11 -0.08 0.042 22 0.08 0.042
12 0.00 0.042 23 0.00 0.042
119