Creating a Code to determine Lift Curve Slope

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

Answer 1

[Star Strider]

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.

[John Shackleford]

TO Star Strider

Thank you very much! I'll do my best to configure the NaN problems.

[Star Strider]

As always, my pleasure!

[kaleab Alamnew]

Hello Srat Strider,

Hope you are doing well. Can you keep continue to solve the problem please? Thank you a lot.