recieving error while using atand command

Hi everyon,

when V0 <10 the code is working without problem.

but when V0 >= 10

Phi is giving NaN values and so a and a_prime are giving missing value.

could you please help me where is the error in my code, or why Phi is giving NaN

Please kindlly find the attachments for excel files also attched

clear all;
clc

%%% Fixed Parameters

B = 2;                      % Number of Blades
R = 10.06;                  % Blade Radius in meter
Rho = 1.225;                % Density of air kg/m3
%V0 =10;                     % Wind Speed  m/s       I Have problem %when I recomment this line
V0 =5;                     % Wind Speed  m/s
w = 7.501;                  % Rotational Speed rad/s
Theta_P = -3;               % Tip Pitch Angle in degree
ac = 0.2;                   % Glauert correction
eps = 1.e-5;                % Error tolerance

%%% Given Blade Data
ns=26;
sections = readmatrix('chord and twist variations.xlsx');
r=sections(:,2);                % Radial Distance in meter
Chord=sections(:,4);            % Chord length in meter
Beta=sections(:,5);             % Twist  in degree


%%%  Airfoil Date 
afoil = readmatrix('airfoil data.xlsx');
Aoa = afoil(:,1);
Cl = afoil(:,2);
Cd = afoil(:,3);

for i = 8:length(r)
    
sigma(i) = (B*Chord(i))/(2*pi*r(i)); 

a=0;
a_prime=0;

alist(i)=[0];
alistprime(i)=[0];

disp(['Section ',num2str(i)])

for n=1:1000
    
Phi =  atand(((1-a)*V0)/((1+a_prime)*w*r(i)));  % Flow Angle     

Theta = Beta(i) + Theta_P ;

aoa = Phi- Theta;                  % local angle of attack in degree

f = (B/2)* ((R-r(i))/(r(i)*sind(Phi)));
F = (2/pi)*acos(exp(-f));           % Prandtl correction factor 

cl=interp1(Aoa,Cl,aoa);
cd=interp1(Aoa,Cd,aoa);

Cn = cl*cosd(Phi)+cd*sind(Phi);
Ct =cl*sind(Phi)-cd*cosd(Phi);

K=(4*F*(sind(Phi))^2)/(sigma(i)*Cn);

%%%%%%%%%% a and ac condition
    
a = 1/(((4*F*(sind(Phi)^2))/(sigma(i)*Cn))+1);

if a < ac || a==ac
a=a;
else
a = 0.5*(2+K*(1-2*ac)-sqrt(((K*(1-2*ac)+2)^2)+4*((K*(ac^2))-1)));
end

a_prime =1/(((4*F*(sind(Phi)*cosd(Phi)))/(sigma(i)*Ct))-1);

alist(i,n+1)=a;           
alistprime(i,n+1)=a_prime;


disp([" a ",a,' a_old ',alist(i,n)," a_prime ",a_prime,' a_Prime_old ',alistprime(i,n)])

%%% condition for epsilon or tolerance

if a-alist(i,n) < eps && a_prime-alistprime(i,n) < eps
    break;
end

end

%%% V relative 
v_rel(i) = sqrt( (V0*(1-alist(end)))^2 + (w*r(i)*(1+alistprime(end)))^2 ); 

%%%  The tangential force per length Pt
Pt(i) = 0.5*Rho*(v_rel(i)^2)*Chord(i)*Ct;

end

Answer 1

[Simon Chan]

The value of variable 'aoa' is outside the interpolation range 'Aoa' and hence it gives NaN as the result.

So the root cause is not coming from Phi.

You may use extrapolation and select your method of doing that.

The following is using linear extrapolation as the method.

cl=interp1(Aoa,Cl,aoa,'linear','extrap');
cd=interp1(Aoa,Cd,aoa,'linear','extrap');

[mehmet salihi]

thank you very much for your help.

this new command helped me.

% I used for my problem
'pchip'