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