If-else statement not working??

Hi,

I am trying to calculate values of Ct before finding axial induction values for a 9x1 array using the following equations:

C_T=(sigma .* ((1-axial_induction).^2) .* ((Cl.*cos(relative_wind))+(Cd.*sin(relative_wind)))) ./ ((sin(relative_wind)).^2);

if C_T<=0.96
    axial_induction=1 ./ (1+(4.*F.*(sin(relative_wind).^2)) ./ (sigma.*Cl.*cos(relative_wind)));  
else
    if C_T>0.96
    axial_induction=1 ./ (((4.*F.*cos(relative_wind)) ./ (sigma.*Cl))-1); 
end
end

However I am getting a warning which states that the variable CT might be set by a non scalar operator. I am unsure how to perform the if else statement for all entities inside the 9x1 vector? I researched a function call ismember but am unsure how to use it? Could anybody please help with this?

[Jasmine]

Ct depends on the value of axial_induction and axial_induction depends on the value of Ct, so how can they evaluate? Does axial_induction have an initial value or something?

[Kevin]

Thank you for your reply. Yes axial induction has an initial value. This is my entire code:

% Inputs
R=0.4;                                          % Radius of Rotor
B=3;                                            % Number of blades
V=2;                                            % Fluid velocity
Rho=998;                                        % Fluid Density
N=9;                                            % Number of Blade Elements
Cp_estimate=0.5;                                % Estimate power coefficient
Alpha_design=4;                                 % Design alpha
Cl_design=1;                                    % Design lift coefficient

% Variables
TSR=1;                                          % Initial tip speed ratio
Cp=0;                                           % Initial power coefficient
i=1;                                            % Counter
alpha_new=0;                                    % Initial value for alpha new
axial_induction=0;                              % Initial axial induction factor
tolerance=0.01;                                 % Tolerance Value
Check=1;                                        % Initial check value
axial_induction_old=0;                          % Initial value for old axial induction factor

Cl=[1.3; 1.1; 1; 0.9; 0.86; 0.83; 0.8; 0.75; 0.5];                                          % Lift Coefficients
Cd=[0.027; 0.024; 0.02; 0.019; 0.018; 0.016; 0.013; 0.012; 0.01];                           % Drag Coefficients

r_local=R/N*(1:9)';
r_over_R=r_local / R;

for TSR=1:10                                     % TSR from 1 to 10
    disp(TSR)                                                                    
    Check=1;
    Cp=0;   
   while abs(Check)>=tolerance

TSR_local=r_over_R .* TSR;

Phi=(2/3)*atan(1./TSR_local);
C=((8*pi.*r_local) ./ (B.*Cl_design)).*(1-cos(Phi));
sigma=(B*C) ./ (pi.*r_local.*2);

axial_induction= 1 ./ (((4.*(sin(Phi).^2)) ./ (sigma.*Cl_design.*cos(Phi)))+1);
angular_induction= (1-(3*axial_induction)) ./ ((4.*axial_induction)-1);

axial_induction_old = axial_induction;        
TSR_local = TSR .* (r_local./R);                                                                 % Local Tip Speed Ratio
Phi = (2/3) .* atan(1./TSR_local);                                                               % Angle of Relative Fluid

relative_wind = atan((1-axial_induction) ./ ((1+angular_induction) .* TSR));

F=(2/pi) .* acos(exp(-(((B/2) .* (1-(r_over_R))) ./ ((r_over_R) .* sin(relative_wind)))));      % Tip Loss Factor

C_T=(sigma .* ((1-axial_induction).^2) .* ((Cl.*cos(relative_wind))+(Cd.*sin(relative_wind)))) ./ ((sin(relative_wind)).^2);

if C_T<=0.96
    axial_induction=1 ./ (1+(4.*F.*(sin(relative_wind).^2)) ./ (sigma.*Cl.*cos(relative_wind)));  
else
    if C_T>0.96
    axial_induction=1 ./ (((4.*F.*cos(relative_wind)) ./ (sigma.*Cl))-1); 
    end
end

D=(8./(TSR.*N)).*(F.*(sin(Phi).^2).*(cos(Phi)-((TSR_local).*(sin(Phi)))).*(sin(Phi)+((TSR_local).*(cos(Phi)))).*(1-(Cd./Cl).*atan(Phi)).*(TSR_local.^2));
Cp=sum(D);

Diff=axial_induction-axial_induction_old;
Check=max(Diff(:));

     end

store_Phi(:,TSR)=Phi;
store_TSR_local(:,TSR)=TSR_local;
store_axial_induction(:,TSR)=axial_induction;
store_angular_induction(:,TSR)=angular_induction;
store_relative_wind(:,TSR)=relative_wind;
store_Check(:,TSR)=Check;
store_Diff(:,TSR)=Diff;
store_Cp(:,TSR)=Cp;
store_TSR(:,TSR)=TSR;

end

figure(1)
plot(store_TSR,store_Cp)
hold all
title('Cp vs Tip Speed Ratio')
xlabel('TSR')
ylabel('Cp')

[Jasmine]

I think the warning is there because you compare Ct, a vector, to 0.96, a scalar. Are you saying if all values of Ct are <= to 0.96, or do you want that bit of code to execute if any value within Ct is less than 0.96?

[Kevin]

I want to execute for each individual value inside the vector

[Jasmine]

So you want the if-statement to execute 9 times for each time Ct is changed? So that it can check each value in Ct against 0.96 and then set axial_induction accordingly? Does this mean that axial_induction should be a vector of length 9 as well?

[Kevin]

Yes thats exactly what I want. I changed my code to try solve for CT but i'm still having no luck with it. Here's my altered code:

% Inputs
R=0.4;                                          % Radius of Rotor
B=3;                                            % Number of blades
V=2;                                            % Fluid velocity
Rho=998;                                        % Fluid Density
N=9;                                            % Number of Blade Elements
Cp_estimate=0.5;                                % Estimate power coefficient
Alpha_design=4;                                 % Design alpha
Cl_design=1;                                    % Design lift coefficient

% Variables
TSR=1;                                          % Initial tip speed ratio
Cp=0;                                           % Initial power coefficient
i=1;                                            % Counter
alpha_new=0;                                    % Initial value for alpha new
tolerance=0.01;                                 % Tolerance Value
axial_induction=[0,0,0,0,0,0,0,0,0];
Check=1;                                        % Initial check value
axial_induction_old=0;                          % Initial value for old axial induction factor

Cl=[1.3; 1.1; 1; 0.9; 0.86; 0.83; 0.8; 0.75; 0.5];                                          % Lift Coefficients
Cd=[0.027; 0.024; 0.02; 0.019; 0.018; 0.016; 0.013; 0.012; 0.01];                           % Drag Coefficients

r_local=R/N*(1:9)';
r_over_R=r_local / R;

for TSR=1:10                                     % TSR from 1 to 10
    disp(TSR)                                                                    
    Check=1;
    Cp=0;

TSR_local=r_over_R .* TSR;    
Phi=(2/3)*atan(1./TSR_local);
C=((8*pi.*r_local) ./ (B.*Cl_design)).*(1-cos(Phi));
sigma=(B*C) ./ (pi.*r_local.*2);

axial_induction= 1 ./ (((4.*(sin(Phi).^2)) ./ (sigma.*Cl_design.*cos(Phi)))+1);
angular_induction= (1-(3*axial_induction)) ./ ((4.*axial_induction)-1);

relative_wind = atan((1-axial_induction) ./ ((1+angular_induction) .* TSR));

F=(2/pi) .* acos(exp(-(((B/2) .* (1-(r_over_R))) ./ ((r_over_R) .* sin(relative_wind)))));      % Tip Loss Factor

C_T=(sigma .* ((1-axial_induction).^2) .* ((Cl.*cos(relative_wind))+(Cd.*sin(relative_wind)))) ./ ((sin(relative_wind)).^2);

     while abs(Check)>=tolerance

axial_induction_old = axial_induction;        
TSR_local = TSR .* (r_local./R);                                                                 % Local Tip Speed Ratio
Phi = (2/3) .* atan(1./TSR_local);                                                               % Angle of Relative Fluid

for i=1:length(C_T)
   if C_T(i) <= 0.96
      axial_induction(i) = 1 / (1+(4*F(i)*(sin(relative_wind(i))^2)) / (sigma(i)*Cl(i)*cos(relative_wind(i)))); 
   else C_T(i) > 0.96
      axial_induction(i) = 1 / (((4*F(i)*cos(relative_wind(i))) / (sigma(i)*Cl(i)))-1);
   end;
end;

D=(8./(TSR.*N)).*(F.*(sin(Phi).^2).*(cos(Phi)-((TSR_local).*(sin(Phi)))).*(sin(Phi)+((TSR_local).*(cos(Phi)))).*(1-(Cd./Cl).*atan(Phi)).*(TSR_local.^2));
Cp=sum(D);

Diff=axial_induction-axial_induction_old;
Check=max(Diff(:));

     end

store_Phi(:,TSR)=Phi;
store_TSR_local(:,TSR)=TSR_local;
store_axial_induction(:,TSR)=axial_induction;
store_angular_induction(:,TSR)=angular_induction;
store_relative_wind(:,TSR)=relative_wind;
store_Check(:,TSR)=Check;
store_Diff(:,TSR)=Diff;
store_Cp(:,TSR)=Cp;
store_TSR(:,TSR)=TSR;

end

figure(1)
plot(store_TSR,store_Cp)
hold all
title('Cp vs Tip Speed Ratio')
xlabel('TSR')
ylabel('Cp')