Creating a Code to determine Lift Curve Slope
조회 수: 104 (최근 30일)
이전 댓글 표시
For an Aerodynamics couse, we need to create a code that computes lifting surfaces (wing) aerodynamics. I cannot get the graph to display a nice CL vs AoA curve, but only a single value. I don't know how I can accomplish this. As I am a novice to MAtlab, I would appreciate some pointers or hints. Thank you.
clc;
% Default airfoil or not...
Default_AF = 0;
% A/F Characteristics:
% Lift coefficient [1/rad]:
Cl_alpha = 5.79;
Cl_alphaInRadians = deg2rad(Cl_alpha);
% Sectional lift curve slope [deg]:
alpha_0 = -2.0;
% Pitching moment [1/rad]:
Cm_alpha = 0;
Cm_alphaInRadians = deg2rad(Cm_alpha);
% 2D Pitching moment:
Cm_0 = -0.025;
% Coefficient in series expansion of circulation distribution [deg]:
alpha_nl = 8;
% Wing planform Geometry [inches]:
% [root-LE,root-TE,tip-LE,tip-TE]
x1 = 0;
y1 = 0;
x2 = 200;
y2 = 0;
x3 = 0;
y3 = 1000;
x4 = 200;
y4 = 1000;
% Wing twist angle [deg]:
theta_t = 0;
% Number of spanwise coefficients:
N = 100;
% Flight Conditions:
% Free-stream Velocity [ft/sec]
V = 100;
% Density [slugs/ft^3]:
rho = 0.002344;
% Viscous Force [lb.s/ft^2]
mu = 0.0000003737;
% Angle of Attack [deg]
alpha = 5;
% Wing span:
b = sqrt((x2-x1)^2 + (y2-y1)^2);
% Wing chord:
c = sqrt((x3-x1)^2 + (y3-y1)^2);
% Wing area:
S = b*c;
% Wing semi-span:
s = b/2;
% Aspect Ratio:
AR = b^2./S;
% Applying segment length of the wings:
alpha_segment = zeros(1,N);
theta_segment = zeros(1,N);
% Taper Ratio and its range within the wing (which is 0 to 1):
for lambda_r = 0:0.25:1
% Root Chord:
c_root = (2*S) / ((1+lambda_r)*b);
% Tip Chord:
c_tip = ((2*S)/((1+lambda_r)*b)) * (1-((2*(1-lambda_r))/b)*(b/2));
% Mean Aerodynamic Chord:
MAC = (2/3)*(c_root+c_tip-(c_root*c_tip)/(c_root+c_tip));
% Taper ratio:
lambda = c_tip / c_root;
% Providing segment limitation within range:
alpha_segment(1) = alpha;
theta_segment(1) = pi/2;
% Completing the Lifting Line Theory to determine the
% matrix value of the segments:
for N = 2:N
c(N) = c_root+(c_tip-c_root)/N.*N;
alpha_segment(N) = alpha+(theta_t)/N.*N;
theta_segment = pi/(2*N):pi/(2*N):pi/2;
end
% Lifting line theory construction:
xb = 4*b./(pi*c);
% Completing the [B] matrix:
for i = 1:N
for j = 1:N
B(i,j) = (sin(j.*theta_segment(i))).*(xb(i)+j/(sin(theta_segment(i))));
end
slope(i) = (alpha_segment(i)-alpha_0)*(pi/180);
end
A = B\transpose(slope);
end
% Leading edge sweep angle
for fai = 1:N
delta_LE(fai) = fai*(A(fai)/A(1))^2; %#ok<*SAGROW>
end
% Section Lift Coefficient of Airfoil
cl = 0.5*cos(pi/b);
% Wing Lift Coefficient:
CL = pi*AR*A(1);
% Span Efficiency:
delta = sum(delta_LE);
CD_0 = 1/(1+delta);
% Induced drag coefficient:
CD_i = CL.^2 / (pi*CD_0*AR);
% Speed of sound (assuming 20 degree dry air) [ft/sec]:
C = 1125.33;
% Mach Number:
M = V / C;
% Dynamic Pressure:
q = (0.5*rho*V^2);
% Pitching Moment Coefficient:
CM = M / q*S*c;
% Subsonic Lift:
alpha_inf = 2*pi*cos(delta_LE);
% 3D Lift Curve Slope of airfoil:
CL_alpha = 2*pi*A;
% Geometric Angle of Attack of Wing:
alpha = (alpha_inf)/(1+(alpha_inf)/(pi*AR));
%----------------------------------------------------------
% Wing lift coefficient (CL) vs Angle of Attack (alpha):
figure(1);
plot(alpha,CL,'-o')
grid on
title('Wing lift coefficient (CL) vs Angle of Attack (alpha)')
ylabel('Wing Lift Coefficient, CL')
xlabel('Angle of Attack, alpha [deg]')
% Pitching moment coefficient(CM) vs Angle of Attack (alpha):
figure(2);
plot(alpha,CM,'-o')
grid on
title('Pitching moment coefficient(CM) vs Angle of Attack (alpha)')
ylabel('Pitching Moment Coefficient, CM')
xlabel('Angle of Attack, alpha [deg]')
% Wing lift coefficient (CL) vs Induced drag coefficient (CD_i):
figure(3);
plot(CD_i,CL,'-o')
grid on
title('Wing lift coefficient (CL) vs Induced drag coefficient (CD_i)')
ylabel('Wing Lift Coefficient, CL')
xlabel('Induced Drag Coefficient, CD_i')
댓글 수: 0
채택된 답변
Star Strider
2019년 12월 3일
There are several problems with your code.
You are only saving the last iteration of ‘A’, and since it does not appear that ‘A’ is calculated recursively, this change in the ‘lambda_r’ loop will save all of them:
lambda_rv = 0:0.25:1;
for k = 1:numel(lambda_rv)
lambda_r = lambda_rv(k);
% Root Chord:
c_root = (2*S) / ((1+lambda_r)*b);
% Tip Chord:
c_tip = ((2*S)/((1+lambda_r)*b)) * (1-((2*(1-lambda_r))/b)*(b/2));
% Mean Aerodynamic Chord:
MAC = (2/3)*(c_root+c_tip-(c_root*c_tip)/(c_root+c_tip));
% Taper ratio:
lambda = c_tip / c_root;
% Providing segment limitation within range:
alpha_segment(1) = alpha;
theta_segment(1) = pi/2;
% Completing the Lifting Line Theory to determine the
% matrix value of the segments:
for N = 2:N
c(N) = c_root+(c_tip-c_root)/N.*N;
alpha_segment(N) = alpha+(theta_t)/N.*N;
theta_segment = pi/(2*N):pi/(2*N):pi/2;
end
% Lifting line theory construction:
xb = 4*b./(pi*c);
% Completing the [B] matrix:
for i = 1:N
for j = 1:N
B(i,j) = (sin(j.*theta_segment(i))).*(xb(i)+j/(sin(theta_segment(i))));
end
slope(i) = (alpha_segment(i)-alpha_0)*(pi/180);
end
A(:,k) = B\transpose(slope);
end
You can eliminate the ‘fai’ loop with:
fai = 1:N;
delta_LE = fai(:).*(A/A(1)).^2; %#ok<*SAGROW>
Note that ‘alpha’ was then calculated as a (100x100) matrix. Since I believe you may want it to be a (100x5) matrix, do element-wise division:
alpha = (alpha_inf)./(1+(alpha_inf)/(pi*AR));
The most serious problem is that ‘alpha’ and many of the rest of your values are all NaN. This is the result of 0/0, Inf/Inf, or existing NaN values. NaN values do not plot.
I defer to you to solve the NaN problems.
댓글 수: 3
kaleab Alamnew
2022년 5월 18일
Hello Srat Strider,
Hope you are doing well. Can you keep continue to solve the problem please? Thank you a lot.
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Airfoil tools에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!