BEM_GPU

%% Exercise 3
% TEP4175
% Design of a Wind Turbine

% Wind turbine calculations using BEM

clear all;
close all;
format long
load airfoildata
clc;

% Fixed parameters
rho=1.225;
V=10;
z=2;

% Variable parameters
Cp_0=0.45;
Rtip=0.45;
Rinner=0.045;
n=1000;
TSR=7.5;

% Airfoil
ALPHA=gpuArray(airfoildata(:,1));
CL=gpuArray(airfoildata(:,2));
CD=gpuArray(airfoildata(:,3));



%Iteration tolerances and limits
itmax_Cp = 30;
tol_Cp = 1e-4;
tol_a_ap = 1e-9;


% Run
profile on
[Cp, Power, Thrust, Torque, r, dT, dM, Lc] = calculate(rho, V, z, TSR,...
    Cp_0, Rtip, Rinner, n, CL, CD, itmax_Cp, tol_Cp, tol_a_ap);
profile off
profile viewer

%1. Maximum power coefficient, Cp [1]
Cp
%2. Maximum power output, P [W]
Power
%3. Thrust force, T [N]
Thrust

%4. Plots
figure;
plot(r/Rtip,dT,r/Rtip,dM)
grid
title('Blade loads')
xlabel('r/R')
ylabel('SI units')
legend('Thrust [N]','Torque [Nm]', 'Location', 'Best')
%5. Plot showing the chord length versus radius
figure;
plot(r/Rtip,Lc)
title('Chord lengths')
grid
xlabel('r/R_{tip}')
ylabel('Chord length [m]')

function [Cp, Power, Thrust, Torque, r, dT, dM, Lc] = calculate(rho, V, z, TSR,...
    Cp_0, Rtip, Rinner, n, CL, CD, itmax_Cp, tol_Cp, tol_a_ap)

    dr=(Rtip-Rinner)/n;
    r=(Rinner+dr/2:dr:Rtip);

    omega=TSR*V/Rtip;
    rpm=30*omega/pi;
    ac = 0.2;

    alphaopt=ALPHA(CL./CD==max(CL./CD));
    Cl0=CL(ALPHA==alphaopt);
    Cd0=CD(ALPHA==alphaopt);

    dA=pi*(((r+dr/2).^2 - (r-dr/2).^2))/z;
    dP=0.5*rho*V.^3*dA*Cp_0;
    dT=dP/omega;
    dFt=dT./r;
    U=omega*r;
    Vrel=sqrt(U.^2+V.^2);
    phi0=(360/(2*pi))*(atan(V./U));
    Lc=(2*dFt)./(rho*dr*Vrel.^2.*(Cl0*sin(((2*pi)/360)*(phi0))-Cd0*cos(((2*pi)/360)*(phi0))));

    it=0;
    Cp=1;
    Cp_prev = 1.1;

%     figure
%     hold on
%     title('C_p convergence')
%     xlabel('Iteration number')
%     ylabel('C_p')
%     grid

    while (abs(Cp-Cp_prev)>tol_Cp)

        it=it+1;

        if it==itmax_Cp
            %disp('No Cp convergence')
            break;
        end

        [a,ap,F,Ca,Cr,phi,Cl,Cd]=a_ap(V,omega,r,CD,CL, ac,n,alphaopt,z,Rtip,Lc, tol_a_ap);

        Bep = (a/(1+ap))*4.*sin((phi));
        Lc= Bep.*(pi*V./(Cl*omega));

        dT=z.*Ca*0.5*rho.*Lc.*((V^2.*(1-a).^2)./sin(phi).^2).*F*dr;
        dM=z.*Cr*0.5*rho.*Lc.*(V.*(1-a)./sin(phi)).*(omega.*r.*(1+ap)./cos(phi)).*F*dr;

        Torque=sum(dM.*r);
        Power=Torque*omega;
        Cp_prev = Cp;
        Cp=Power/(0.5*rho*V^3*pi*Rtip^2);
        Thrust = sum(dT);

%         datavec(it,1) = Cp;
%         datavec(it,2) = abs(Cp-Cp_prev);
%         error = abs(Cp-Cp_prev);
%         pause(0.0001)
%         plot(datavec(:,1), '-r' )

    end
end


function [avec,apvec,Fvec,Cavec,Crvec,phivec,Clvec,Cdvec,thetavec]=...
        a_ap(V,omega,r,CD,CL,ac,n,alphaopt,z,Rtip,Lc, tol_a_ap)


    a0=0;
    ap0=1;
    sigma=z.*Lc./(2*pi*r);

    avec = zeros(1,n);
    apvec = zeros(1,n);
    Fvec = zeros(1,n);
    Cavec = zeros(1,n);
    Crvec = zeros(1,n);
    phivec = zeros(1,n);
    Clvec = zeros(1,n);
    Cdvec = zeros(1,n);
    thetavec = zeros(1,n);

    for i=1:n

        a=1/3;
        ap=0;

        while (abs(a0-a)>tol_a_ap) && (abs(ap0-ap)>tol_a_ap)

            phi=atan((1-a)*V/((1+ap)*omega*r(i)));

            theta=(360/(2*pi))*(phi)-alphaopt;
            alpha = (360/(2*pi))*(phi)- theta;

            Cl=interp1q(ALPHA,CL,alpha);
            Cd=interp1q(ALPHA,CD,alpha);


            Ca=Cl*cos(phi)+Cd*sin(phi);
            Cr=Cl*sin(phi)-Cd*cos(phi);

            F=2/pi*acos(exp(-z*(Rtip-r(i))/(2*r(i)*sin(phi))));
            K=4*F*sin(phi)^2/(Ca*sigma(i));

            a0=a;
            ap0=ap;
            if a<ac
                a=1/(K+1);
            else
                a=0.5*(2+K*(1-2*ac)-sqrt((K*(1-2*ac)+2)^2+4*(K*ac^2-1)));

            end
            ap=1/(4*F*sin(phi)*cos(phi)/(Cr*sigma(i))-1);
        end

        avec(i) = a;
        apvec(i) = ap;
        Fvec(i) = F;
        Cavec(i) = Ca;
        Crvec(i) = Cr;
        phivec(i) = phi;
        Clvec(i) = Cl;
        Cdvec(i) = Cd;
        thetavec(i) = theta;

    end
end