Convergence Error in Matlab Code for Fluid Flow Simulation

[HACRS4]

Homework Statement

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).

Homework Equations

Code based on the one presented in the following document:

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

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

end

end

A=B\transpose(LHS);

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

end
end

% 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
plot(y_s,CL1,'-o')
elseif Lambda == 0.5
plot(y_s,CL1,'-*')
elseif Lambda == 0.75
plot(y_s,CL1,'-s')
else
plot(y_s,CL1,'-d')
end
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;

end

delta1 = sum(delta);

e = 1 / (1 + delta1)

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

end


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

end

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')

 

[jedishrfu]

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.

 

[HACRS4]

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.

 

[prada]

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

 

[SteamKing]

Apologies, but I'm not sure how to go about doing that. I am a Matlab, and coding in general, novice.

It's pretty simple to add

 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.

 

[AlephZero]

Don't hold your breath waiting for a response. Check the dates of the previous posts

 

[SteamKing]

Ah, the necropost!

 

[gsb]

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

end

plot((1:t), q);

end

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.
thanks