1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Lifting Line Code in Matlab

  1. Nov 24, 2013 #1
    1. The problem statement, all variables and given/known data

    My apologies if the error is an obvious one and many thanks in advance for any help given:

    I am trying to create a Matlab code that simulates Lifting Line Theory in order to provide an estimate of the lift and drag of a 3D wing. My hope is to later use this as part of an optimization routine for the wing design.

    The plot showing the variation of the lift distribution with taper ratio appears to be the wrong way round. That is, it shows that the lift generated at the wing tip is greater for a smaller taper ratio. I thought that it should be the other way around (more taper, less lift generated at the wing tip).

    2. Relevant equations

    Code based on the one presented in the following document:

    http://faculty.dwc.edu/sadraey/Chapter 5. Wing Design.pdf

    3. The attempt at a solution

    % Lifting Line Code

    clc;clear all;close all

    % Input Wing Geometry

    for Lambda = 0.25:0.25:1

    % Number of wing segments
    N_seg = 30;

    % Wing Area [m^2]
    S = 65;

    % Aspect Ratio
    AR = 25;

    % Twist [deg]
    twist = -6;

    % Wing Setting Angle/ AoA [deg]
    set_ang = 2;

    % Wingspan [m]
    b = sqrt(AR*S);

    % Input Aerofoil Data

    % Lift Curve Slope [1/rad]
    slope = 6.9;

    % Mean Aerodynamic Chord
    MAC = S/b;

    % Zero Lift Angle of Attack
    zl_AoA = -6.56;

    % root chord (m)
    Root_Chord = (2*S) / ((1 + Lambda) * b);

    % tip chord (m)
    Tip_Chord = ((2*S) / ((1 + Lambda) * b)) * (1 - ((2* (1 - Lambda)) / b) * (b/2) );

    % MAC
    MAC = (2/3) * (Root_Chord + Tip_Chord - (Root_Chord*Tip_Chord) / (Root_Chord +Tip_Chord));

    % Lifting Line Algorithm

    theta = pi/(2*N_seg):pi/(2*N_seg):pi/2;

    % create vector containing each segment's angle of attack
    alpha = set_ang + twist: -twist/(N_seg-1):set_ang;

    z = (b/2)*cos(theta);

    % Mean Aerodynamic Chord at each segment (m)
    c = Root_Chord * (1 - (1-Lambda)*cos(theta));

    mu = c * slope / (4 * b);

    LHS = mu .* (alpha-zl_AoA) * (pi/180);

    % Determine Coefficients A(i) by Solving N_seg Equations:

    for i=1:N_seg

    for j=1:N_seg

    B(i,j) = sin((2*j-1) * theta(i)) * (1 + (mu(i) * (2*j-1)) / sin(theta(i)));




    for i = 1:N_seg

    sum1(i) = 0;
    sum2(i) = 0;

    for j = 1 : N_seg

    sum1(i) = sum1(i) + (2*j-1) * A(j)*sin((2*j-1)*theta(i));
    sum2(i) = sum2(i) + A(j)*sin((2*j-1)*theta(i));


    % Determine Lift Coefficient For Each Segment

    CL = 4*b*sum2 ./ c;

    % Plot Lift Distribution

    CL1=[0 CL];
    y_s=[b/2 z];

    if Lambda == 0.25
    elseif Lambda == 0.5
    elseif Lambda == 0.75
    hold on

    % Output Lift Coefficient for Wing

    CL_wing = pi * AR * A(1)

    % Output Oswald Efficiency Factor Induced Drag Coefficent

    for n = 2:30

    delta(n) = n * (A(n)/A(1))^2;


    delta1 = sum(delta);

    e = 1 / (1 + delta1)

    CD_induced = CL_wing^2 / (pi * e * AR)


    % Plot Elliptical Lift Distribution For Comparison

    y_ellipse = linspace(0,(b/2),100);

    for i = 1:length(y_ellipse)

    CL_ellipse(i) = CL1(end) * sqrt(1 - ((2*y_ellipse(i)) / b)^2 );


    plot (y_ellipse, CL_ellipse, 'r-')

    grid on
    title('Variation of Lift distribution with Taper Ratio')
    xlabel('Semi-Span Location [m]')
    ylabel ('C_L')
    legend('location','best','Lambda = 0.25','Lambda = 0.5','Lambda = 0.75','Lambda = 1','Elliptical Lift Distribution')
  2. jcsd
  3. Nov 24, 2013 #2


    Staff: Mentor

    could you please edit your post and bracket your code listing with code tags?

    It will look much nicer on the eyes and will preserve indentations.

    Also you haven't said what the error is.
  4. Nov 24, 2013 #3
    Apologies, but I'm not sure how to go about doing that. I am a Matlab, and coding in general, novice.

    The error is in the plot that is generated - higher taper ratios have a better efficiency (closer to elliptical lift distribution) which I don't think should be the case. Many thanks for any help you can provide.
    Last edited: Nov 24, 2013
  5. Aug 20, 2014 #4
    Hi HACRS4,
    I'm not sure that i perfectly understand your question. I have tried to use your code, and I noticed that the efficiency goes from a minimum to a maximum and then decreases again while lambda goes from 1 to 0.25. This is exactly what is suppose to happen. In fact for a fixed aspect ratio wing the maximum efficiency should happen at a tapper ratio of approximately 0.4.
    Don't know if it was of any help or if it was too late. But I have actually used your code in my master thesis to do some simple calculations, and I also added the functionality to calculate the best linear wash-out geometric twist. Unfortunately I can't make a proper reference to your work. If you are so kind to tell me who to reference to I will.
    Thanks you in advance
  6. Aug 20, 2014 #5


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    It's pretty simple to add
    Code (Text):
     tags to your reply.

    You can hit the 'Go Advanced' tab at the bottom of the page and select the icon at the top which looks like a sheet of paper with a '#' on it.

    To check what you've done, hit the 'Preview Post' button and check the post before hitting the 'Submit Reply' button.
  7. Aug 20, 2014 #6


    User Avatar
    Science Advisor
    Homework Helper

    Don't hold your breath waiting for a response. Check the dates of the previous posts :smile:
  8. Aug 20, 2014 #7


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Ah, the necropost!
  9. Aug 26, 2014 #8


    User Avatar

    hi, i am having some problem with a code in matlab. my code is:

    t = 30;
    H= 30;
    h= 3;
    n = H/h;
    qqq= 10;
    a= 3;
    b = .2;
    k= 0.04;
    Lo= 100;
    g= 9.81;

    for i = 1:t
    q(i) = 0;

    for j = 1:n
    nl = qqq + (j*h)/(a + b*(j*h));
    nlt = nl*(0.0172*log(i) + 1);
    q(i) = q(i) + 2*k*(Lo/100)*(nlt/3)*g*exp(-k*(i));


    plot((1:t), q);


    this is the code and i keep getting this error- "Attempted to access q(2); index out of bounds because numel(q)=1."

    help would be greatly appreciated.
    Last edited: Aug 26, 2014
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted