Calculating lift with lifting-line theory in MATLAB

Question

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

Answer 1

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:
AoA = deg2rad(0);
AoA_0 = deg2rad(-6.4);

% 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;