I need help with elemental wise operation for a first point iteration method

I have created a BEM (Blade Element Momentum) program in MATLAB, I am not getting the correct results when I run the program.

My Program:

%% Input Variable

Uinf  = 11;             % freestream veloctiy (m/s)
R     = 6.25;           % Radius of Rotor (m)                               % Input 1
N     = 90;             % Rotor Speed (rpm)
N_rad = 2*pi*N/60;      % Rotor Speed (rad/s)

%% Blade Segment

r       = [1.25; 1.56];     % segments (m)                                  % Input 2
c       = [0.70; 0.63];     % chord length of the airfoil at segement (m)   % Input 3
theta_t = [17.46; 14.22];   % twist angle (deg)                             % Input 4

%% Initialization Vairable

a     = zeros(length(r),1); % axial induction factor inital assumption
at    = zeros(length(r),1); % angular induction factor
anew  = zeros(length(r),1);
atnew = zeros(length(r),1);
B     = 3;                  % no of blades
rho   = 1.225;              % density of air @ 25 degCel

%% Setup Details

n = 400;    % no of iterations
e = 1e-03;  % tolerance limit

%% Problem Solution
% For Loop for cycling through each blade element
for j = 1:length(r)
    % For Loop for no of iterations
    for i= 1:n
        %% Disk Velocity and wake Velocity
        Vd = (1-a(j)).*Uinf;
        Vw = (1-2*a(j)).*Uinf;
        
        %% relative velocity and angle of attack determination
        t       = (Uinf/(N_rad.*r(j)))*((1-a(j))./(1+at(j)));
        phi     = atan(t);          % rad
        phi_d   = 180/pi()*phi;     % deg
        alpha_d = phi_d - theta_t;  % deg
        alpha   = alpha_d*pi()/180; % rad
        
        %% Coefficient of Lift and Drag,  and Relative Velocity
        Cl = -0.0091*(alpha_d).^2 + 0.2695*(alpha_d) - 0.1205;
        Cd = -0.0001*(alpha_d).^2 - 0.0051*alpha_d + 0.0527;
        Vr = sqrt(((Uinf.*(1 - a(j))).^2) + ((N_rad*r(j)*(1 + at(j))).^2)); % relative velocity
        
        %% Normal Force Coefficient and Tangential Force Coefficient
        Cn = Cl*cos(phi) + Cd*sin(phi);     % normal force coefficient
        Ct = Cl*sin(phi) - Cd*cos(phi);     % tangential force coefficient
        
        %% Pradntl Tip Loss Factor
        f = (B*(R-r(j)))./(2.*r(j)*sin(phi));
        F = (2/pi)*(acos(exp(-f)));
        
        %% Solidity Constant
        sigmar = (B*c(j))/(2*pi.*r(j));
        Za     = sigmar*Cn + 4*F*(sin(phi)).^2;
        Zt     = N_rad.*r(j)*Cn + 8*F*pi.*r(j)*(sin(phi)).^2;
        
        %% Axial and Tangential Induction Factor re-calculation
        anew   = (sigmar*Cn)/Za;
        atnew  = (Uinf*Ct)/Zt;
        
        %fprintf('a%d = %.4f\n',i,anew)
        %fprintf('at%d = %.4f\n',i,atnew)
        %% Test Tolerance Limit for Loop Break
        if abs(anew(j)-a(j))<e & abs(atnew(j)-at(j))<e
            break;
        end
        %% Reassigning Induction factor for next iteration of loop    
        a(j)  =  anew(j);
        at(j) = atnew(j);
    end
end

%% OUTPUTS
a

a = 2×1

0.0671 0

at

at = 2×1

0.3556 0

OUTPUT (I NEED):

(a matrix) 0.322

0.3115

(at matrix) 0.3222

0.2404

OUTPUT (I GET):

(a matrix) 0.06712

0

(at matrix) 0.3556

0

This is an iterative approach like fixed point iteration method to determine the convergence for the variable a and at but all I need is to save the final result for a and at in matrix form. I need help debugging the the mistake in the code.

Answer 1

[Torsten]

%% Input Variable

Uinf  = 11;             % freestream veloctiy (m/s)
R     = 6.25;           % Radius of Rotor (m)                               % Input 1
N     = 90;             % Rotor Speed (rpm)
N_rad = 2*pi*N/60;      % Rotor Speed (rad/s)

%% Blade Segment

r       = [1.25; 1.56];     % segments (m)                                  % Input 2
c       = [0.70; 0.63];     % chord length of the airfoil at segement (m)   % Input 3
theta_t = [17.46; 14.22];   % twist angle (deg)                             % Input 4

%% Initialization Vairable

a     = zeros(length(r),1); % axial induction factor inital assumption
at    = zeros(length(r),1); % angular induction factor
B     = 3;                  % no of blades
rho   = 1.225;              % density of air @ 25 degCel

%% Setup Details

n = 400;    % no of iterations
e = 1e-03;  % tolerance limit

%% Problem Solution
% For Loop for cycling through each blade element
for j = 1:length(r)
    % For Loop for no of iterations
    for i= 1:n
        %% Disk Velocity and wake Velocity
        Vd = (1-a(j)).*Uinf;
        Vw = (1-2*a(j)).*Uinf;
        
        %% relative velocity and angle of attack determination
        t       = (Uinf/(N_rad.*r(j)))*((1-a(j))./(1+at(j)));
        phi     = atan(t);          % rad
        phi_d   = 180/pi()*phi;     % deg
        alpha_d = phi_d - theta_t(j);  % deg
        alpha   = alpha_d*pi()/180; % rad

        %% Coefficient of Lift and Drag,  and Relative Velocity
        Cl = -0.0091*(alpha_d).^2 + 0.2695*(alpha_d) - 0.1205;
        Cd = -0.0001*(alpha_d).^2 - 0.0051*alpha_d + 0.0527;
        Vr = sqrt(((Uinf.*(1 - a(j))).^2) + ((N_rad*r(j)*(1 + at(j))).^2)); % relative velocity
        
        %% Normal Force Coefficient and Tangential Force Coefficient
        Cn = Cl*cos(phi) + Cd*sin(phi);     % normal force coefficient
        Ct = Cl*sin(phi) - Cd*cos(phi);     % tangential force coefficient
        
        %% Pradntl Tip Loss Factor
        f = (B*(R-r(j)))./(2.*r(j)*sin(phi));
        F = (2/pi)*(acos(exp(-f)));
        
        %% Solidity Constant
        sigmar = (B*c(j))/(2*pi.*r(j));
        Za     = sigmar*Cn + 4*F*(sin(phi)).^2;
        Zt     = N_rad.*r(j)*Cn + 8*F*pi.*r(j)*(sin(phi)).^2;
        
        %% Axial and Tangential Induction Factor re-calculation
        anew   = (sigmar*Cn)/Za;
        atnew  = (Uinf*Ct)/Zt;
        
        %fprintf('a%d = %.4f\n',i,anew)
        %fprintf('at%d = %.4f\n',i,atnew)
        %% Test Tolerance Limit for Loop Break
        if abs(anew-a(j))<e & abs(atnew-at(j))<e
            break;
        end
        %% Reassigning Induction factor for next iteration of loop    
        a(j)  =  anew;
        at(j) = atnew;
    end
end
a

a = 2×1

0.3222 0.3113

at

at = 2×1

0.3228 0.2410

[Mohamed Farmaan]

Thank you very much, can you tell me where I went wrong, and how it affects the code.

[Torsten]

alpha_d = phi_d - theta_t(j);  % deg

instead of

alpha_d = phi_d - theta_t;  % deg

Removal of

anew  = zeros(length(r),1);
atnew = zeros(length(r),1);

and use of the two variables as scalars instead of arrays.