Numerical Nonlinear Lifting Line Theory in MA

[meltems]

1. Homework Statement
Hello all.
It is not a homework actually. I just didn't know at which forum I should post. I am working on a MATLAB code solving the finite wing properties iteratively by using the Anderson's Numerical Lifting Line Method. However, I got some wrong results. The circulation (c) values do not seem right and the code found the induced angle of attacks have negative values. I checked all equations. I think I may have mistakes at the structure of the code or while defining the initial values.

Many thanks in advance for any help given.

2. Homework Equations
The paper (3.44690) including the algorithm I follow is attached.

The Attempt at a Solution

The MATLAB code (nll) is also attached.

function nll
format longg
span=2;
ym=-span/2; %%%general lift distribution
yp=span/2;
noseg=101;
dy=(yp-ym)/(noseg-1);
vel=30;
locc=1000000; %%bu neden 1000000
% cdif=1;
alfa=0.1;

for i=1:noseg
y(i)=ym+(i-1)*dy;
c(i)=10*sqrt(1-(2*y(i)/span)^2); %%% c dediğimiz circulation. 10 bizim c0 ımız

end
y';
c';
%%
for i=1:(noseg-1)/2.+1
cbar(i)=0.3+(0.004*(i-1));
%cbar(i)=0.4;
end

for i=1:noseg
cbar(noseg+1-i)=cbar(i);
%cbar(i)=0.1;
end
%%%%%%%%%%%%%
cbar';

iter=0;
while (iter<=200)
iter=iter+1;
% cdif=0;
 
totint=0;
for j=1:noseg
 
 derc=0;
 down=0;
 intpart=0;
 
for i=1:noseg
if (i==1)
derc(i)=(c(i+1)-c(i))/(y(i+1)-y(i));
down(i)=y(j)-y(i);
intpart(i)=derc(i)/down(i);
elseif (i==noseg)
derc(i)=(c(i)-c(i-1))/(y(i)-y(i-1));
down(i)=y(j)-y(i);
intpart(i)=derc(i)/down(i);
    else
derc(i)=(c(i+1)-c(i-1))/(y(i+1)-y(i-1));
down(i)=y(j)-y(i);
intpart(i)=derc(i)/down(i);   
    end    
% locc;
if(abs(down(i))<(dy/10))
locc=i;
end
% locc;

end
% j;
% locc;

if (locc==1)
intpart(locc)=intpart(locc+1);
elseif(locc==noseg)
    intpart(locc)=intpart(locc-1);
else
        intpart(locc)=(intpart(locc-1)+intpart(locc+1))/2.;
end

% locc=10000000;
totint(j)=0;
for i=2:2:noseg-1;
totint(j)=totint(j)+dy/3*(intpart(i-1)+4*intpart(i)+intpart(i+1));
end

% totint(j);
aind(j)=totint(j)/(4.*acos(-1.)*vel);
aeff(j)=alfa-aind(j);
cl(j)=6.917*(aeff(j));

cnew(j)=0.5*vel*cbar(j)*cl(j);
j;
end
cold=c;
cold';
for i=1:noseg
c(i)=cold(i)+0.05*(cnew(i)-cold(i));
end
c';
aind';
cdif=1;
end

aind';
Lift=0;
indrag=0;
for i=2:2:noseg-1
Lift=Lift+dy*(c(i-1)+4*c(i)+c(i+1))/3.;
% indrag(2)=dy*(c(2-1)*aind(2-1)+4*c(2)*aind(2)+c(2+1)*aind(2+1))/3.;
indrag=indrag+dy*(c(i-1)*aind(i-1)+4*c(i)*aind(i)+c(i+1)*aind(i+1))/3.;
% indrag2(i)=indrag(i)*1.225*vel
end
Lift=Lift*1.225*vel;
indrag=indrag*1.225*vel
S=(cbar(1)+cbar((noseg-1)/2+1))/2.*span;
AR=span^2/S;
CL=Lift*2/(1.225*vel^2*S);
% cdi=indrag*2/(1.225*vel^2*S);
% oswald=CL^2/cdi/acos(-1.)/AR;

end

Added code tags jmc

 

[jim mcnamara]

Code HTML tags added to improve readability.

 

[magoo]

I couldn't read the paper other than the references. The rest was blank.

 

[meltems]

I couldn't read the paper other than the references. The rest was blank.

Oh, I can see the rest of it. Here is the algorithm I followed.

 

[meltems]

Code HTML tags added to improve readability.

Thank you