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')

## Best Answer