function coursework2009_v01
% this function performs calculation of flapwise stability of a 3-bladed
% horizontal axis wind turbine with constant chord length
%
% governing equations from 'Wind Energy Explained' page 181
tic
disp(' ')
disp('>>> Begin flapping angle calculation')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculation parameters
%
finalN  = 50;        % number of cycles
psi0    = 0;        % initial azimuth angle
beta0   = 5;        % initial flapping angle
alpha0  = 0;        % initial flapping velocity

% Operation parameter
tsr     = 5;        % tip speed ratio
gamma   = 01;       % [deg] - yaw error angle
q       = 0.01;     % [deg/sec] - yaw rate
a       = 1/3;      % axial induction factor
thetap  = 9.0;      % [deg] - pitch angle of untwisted blade

% Atmospheric parameters
rho     = 1.165;    % [kg/m^3] - air density
V       = 10;      % [m/sec] - mean air speed
g       = 9.81;     % [m/sec^2] - acceleration due to gravity
Kvs     = 0.00;     % vertical wind shear linear constant

% Structural parameters
c       = 0.30;     % [m] - chord length
R       = 2.50;     % [m] - rotor RADIUS
mb      = 10;       % [kg] - mass of blade
omegaR  = 325;      % [Hz] - ROTATING spring natural frequency
omegaNR = 0.9*omegaR;       % [Hz] - NONROTATING spring natural frequency
dyaw    = 0.25;     % [m] - yaw moment arm
Dhub    = 0.40;     % [m] - hub diameter


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculated aeroelastic system parameters

dbar    = dyaw/R;               % non-dimensional yaw moment arm
gamma   = gamma*pi/180;         % [rad] - yaw error angle
V1      = V*cos(gamma);         % [m/sec] - axial wind velocity
V0      = V*sin(gamma);         % [m/sec] - cross wind due to yaw error
OMEGA   = tsr*V1/R;             % [rad/sec] - blade angular velocity
q       = q*pi/180;             % [rad/sec] - yaw rate
qbar    = q/OMEGA;              % non-dimensional yaw rate
V0bar   = V0/(OMEGA*R);         % non-dimensional cross flow due to yaw error
V1bar   = V1/(OMEGA*R);         % non-dimensional axial flow
z       = (omegaR^2-omegaNR^2)/OMEGA^2;
e       = 2*(z-1)/(3+2*(z-1));  % non-dimensional offset
Ib      = mb*((R^3)/3)*(1-e)^3;     % area moment of inertia of hinged blade
Kb      = Ib*omegaNR^2;         % torsional spring stiffness
epsilon = 3*e/(2*(1-e));        % offset term
K       = 1+epsilon+Kb/(Ib*OMEGA^2); % Flapping inertia natural frequency
rg      = 0.5*Dhub + 0.5*(R-Dhub/2); % radial distance to blade CG
G       = g*mb*rg/Ib;
B       = G/(2*OMEGA^2);        % gravity term
LAMBDA  = V1*(1-a)/(OMEGA*R);   % non-dimensional inflow
A       = (LAMBDA/3)-(thetap/4);    % Axisymmetric flow term #1
A3      = (LAMBDA/3)-(2*thetap/3);  % Axisymmetric flow term #3

% Calculated integration parameters
finalT  = finalN*2*pi;          % range of azimuth angle rotation


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parsing parameters to ODE solver
inputs.gamma    = gamma;    inputs.V0bar    = V0bar;    inputs.qbar = qbar;
inputs.dbar     = dbar;     inputs.K        = K;        inputs.B    = B;
inputs.A        = A;        inputs.A3       = A3;       inputs.Kvs  = Kvs;      
inputs.V1bar    = V1bar;    inputs.thetap   = thetap;

% Solve ODEs using ODE45
options = odeset('RelTol',1e-09,'AbsTol',1e-09);
[psi,X] = ode45(@flapping_eqn, [0 finalT],[alpha0 beta0], options, inputs);

alpha   = X(:,1)*180/pi;    % [deg/(azimuth angle in rad)] - flapping speed
beta    = X(:,2)*180/pi;    % [deg] - flapping angle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% plotting results

RSTARTREV	= 2500;		% When does steady state start

Rstop   = length(psi);
if RSTARTREV < finalN
	psiR	= psi/(2*pi);
	Rstart	= min(find(psiR>RSTARTREV));
else
	Rstart  = floor(0.20*Rstop);
end
range   = Rstart:Rstop;


figure(1)
clf
plot(psi(range)/(2*pi),beta(range))
xlabel('Rotor revolution')
ylabel('Flapping angle (deg.)')

figure(2)
clf
plot(alpha(range),beta(range))
xlabel('Flapping speed')
ylabel('Flapping angle')

toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function dy = flapping_eqn(psi,y,p)

K(1,1)  = -(p.gamma/8)*(1-(4/3)*cos(psi)*(p.V0bar+p.qbar*p.dbar));
K(1,2)  = -(p.K+2*p.B*cos(psi)+(p.gamma/6)*p.V0bar*sin(psi));
K(2,1)  = 1;
K(2,2)  = 0;
F(1,1)  = p.gamma*p.A/2 - (p.gamma*p.qbar/8)*sin(psi) - ((2*p.qbar + ...
    (p.gamma/2)*(p.A3*(p.V0bar+p.qbar*p.dbar)+(p.Kvs*p.V1bar/4)))*cos(psi));
F(2,1)  = 0;

dy  = K*y + F;