Calculating lift with lifting-line theory in MATLAB

I am attempting to calculate the lift generated by a model wing (~0.8$$m^2$$ wing area) using Prandtl's lifting-line theory in MATLAB. I am starting with the equation below from Anderson's 2005 book Fundamentals of Aerodynamics:

I have defined the $$\frac{2b}{\pi c(\theta)}$$ as '$$\mu$$' for simplicity, moved this back into the first summation and rearranged the above to:

$$\Sigma A_n\sin(n\theta)(\mu + \frac{n}{sin(\theta)}) = \alpha(\theta) - \alpha_0(\theta)$$

Presumably if I split the half-wing section (i.e from $$\theta = \pi/2$$ to $$\pi$$) into 10 segments I can generate a set of simultaneous equations from $$n = 1$$ to $$10$$ and $$\theta_1$$ to $$\theta_{10}$$. Solving that set of equations should then give me values of $$A_1$$ through to $$A_{10}$$. However, the values I'm generating for my coefficient matrix, $$B$$, (i.e. the $$\sin(n\theta)(\mu...)$$ bit) seem to be off, reaching 65 for $$\theta_1$$ and $$n =10$$. This is turn gives me incorrect values of the $$A$$ coefficients. Both are shown below:

This leads me to believe something must be off in either my re-arrangement above or my MATLAB code, which is shown below for completeness. It might be inefficient but please ignore that for now! Anyone able to assist with where I might be going wrong? Would be greatly appreciated!

MATLAB Code:

clear all
clc

% Wing segments:
N_seg = 10;
N_seg2 = N_seg;

% Flow:
U_0 = 15.433 % [m/s]
sweep = deg2rad(30); % Sweep angle
U_n = U_0 * cos(sweep); % Flow component travelling normal to wing sweep

% AoA and zero-lift AoA:
AoA = 0;
AoA_0 = -6.4;

% Wing twist, to be confirmed following further design:
twist = 0;

% Initiliase segment chord, theta and AoA matrix:
c = zeros(1,N_seg);
theta = zeros(1,N_seg);
AoA_seg = zeros(1,N_seg);

% Lambda range:
for Lambda = 0.25:0.25:0.5

% Wing area, span, chord etc.
S = 0.8; % [m^2], wing area
b = 1.56827; % Wing span
s = b/2; % Semi-wing span
c_root = (2*S) / ((1 + Lambda) * b); % root chord [m]

c(1) = c_root; % setting root chord as inner segment wing chord
theta(1) = pi/2;
AoA_seg(1) = AoA;

% tip chord [m]
c_tip = ((2*S) / ((1 + Lambda) * b)) * (1 - ((2* (1 - Lambda)) / b) * (b/2) );
% mean aerodynamic chord [m]
MAC = (2/3) * (c_root + c_tip - (c_root * c_tip) / (c_root + c_tip))

% Segment chord lengths and theta to wing tip:
for N_seg = 2:N_seg
% Chord length of each segment, linear interpolation
c(N_seg) = c_root + (c_tip - c_root)/N_seg2 .* N_seg
% Theta of each segment;- first segment must not be 0,
% otherwise singular matrix obtained
theta = pi/(2*N_seg):pi/(2*N_seg):pi/2;

% AoA of each segment , linear interpolation
AoA_seg(N_seg) = AoA + (twist)/N_seg2 .* N_seg
end
mu = 2*b./(pi*c);

% Setting up RHS (B) and LHS matrix:
for i = 1:N_seg
for j = 1:N_seg
B(i,j) = (sin(j .* theta(i))) .* (mu(i) + j/(sin(theta(i))))
end
LHS(i) = (AoA_seg(i) - AoA_0) * (pi/180);
end

% Finding A coefficients:

A = B\transpose(LHS)

end


• If you suspect a problem in your MatLab code you might do well to post that part of this question on Stack Overflow – CatchAsCatchCan Nov 13 at 22:34

Your general formulation seems fine. Two things to point out:

1. If you choose to sample only half the wing, then you cannot utilize the relation: $$C_L=A_1\pi A$$. You would have to do the integration: $$C_L=\frac{2}{V_\infty S}\int^{b/2}_{-b/2}{\Gamma(y)dy}$$. That being said, maybe you can simplify the expression and get a closed-form result? Haven't tried myself.

2. Your code doesn't seem to handle for the case where $$\sin(\theta)=0$$.

You can take a look at the code below for reference:

% Wing area, span, chord etc.
S = 0.8; % [m^2], wing area
b = 1.56827; % Wing span
s = b/2; % Semi-wing span
taper = 0.25;
Cr = S/s/(1+taper);
Ct = taper*Cr;
MAC = 2/3*Cr*(1+taper+taper^2)/(1+taper);
AR = b^2/S;

% AoA and zero-lift AoA:

% Discretization
num = 10;
[N,Ts] = meshgrid(1:num, linspace(0,pi,num));
ys = -s*cos(Ts);
Cs = (Ct-Cr)/s*abs(ys)+Cr;

% Construct linear system
B = 2*b/pi./Cs.*sin(N.*Ts)+N.*sin(N.*Ts)./sin(Ts);
if Ts(1)==0
% Singular case, use l'Hopital's rule
B(1,:) = N(1,:).*N(1,:);
end
RHS = AoA - ones(num,1)*AoA_0;

% Solve
A = B\RHS;
CL = A(1)*pi*AR;