3
$\begingroup$

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:

enter image description here

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:

enter image description here

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

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

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;
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.