Code:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* this is a simplified code using BEM equations */
int main()
{
int i; // counter for referring to rows of array
int Nb = 3; // number of blades
float pi=3.141592654; // pi
float R=1; // radius of turbine (m)
float tsr=7; // tip speed ratio
float rho=1.225; // density of air (kg/m^3)
float U=10; // free stream velocity of wind (m/s)
float maxchord=0.25; // maximum allowed chord length (m)
float aguess = 0; // initial guess for axial induction factor
float aprimeguess = 0; // initial guess for tangential induction factor
float diff1 = 1; // variable showing difference between axial induction factor and current guess
float diff2 = 1; // variable showing difference between tangential induction factor and current guess
float delta_r = 0.01; // this is the width of one blade element (m)
float opt_deltaQ = 0; // optimum values of deltaQ
float alpha_deg, alpha_rad, Cl, Cd, Vr, chord, idealchord, phi, beta_rad, beta_deg, PTLfactor, a, aprime, solidity, deltaQ, omega, opt_alpharad, opt_phi, opt_chord, dimRi;
float array[19][3]; // this declares an array
FILE * pFile; // declares file
char mystring [1000]; // declares string to write file into
int j=0; // counter for incrementing row of array to copy string into
pFile = fopen ("c:\\Users\\Shaun\\My Documents\\flat plate info.txt" , "r"); // this opens the specified file
if (pFile == NULL) perror ("Error opening the file"); // gives error message if file can't be opened
else {
while(fgetc(pFile) != EOF){ // runs loop until end of file
if ( fgets (mystring , 300 , pFile) != NULL ) // gets string from the file provided there is no error
sscanf (mystring, "%f %f %f",&array[j][0], &array[j][1], &array[j][2]); // takes values from the string and puts them into the array
j++; // increments to next row of array
}
fclose (pFile); // closes file after use
}
for( dimRi=0.1; dimRi<1.01; dimRi=dimRi+0.05 ) //runs the code for various points along the radius
{
for (i=0; i<19; i++ ) // this runs through each set of values for alpha, Cl and Cd
{
alpha_deg = array[i][0]; // angle of attack, alpha, in degrees
alpha_rad=alpha_deg*pi/180; // angle of attack, alpha, in radians
Cl = array[i][1]; // coefficient of lift
Cd = array[i][2]; // coefficient of drag
while(diff1 > 0.0001 || diff2 >0.0001) // this runs the loop until the values for a and aprime have converged
{
omega=tsr*U*(1-aguess)/R; // angular velocity of blade (rad/s)
phi = atan((U*(1-aguess))/(omega*dimRi*R*(1+aprimeguess))); // incoming flow angle (rad)
Vr = (U*(1-aguess))/sin(phi); // resultant wind velocity (m/s)
idealchord = (16*pi*U)/(9*Nb*Cl*tsr*dimRi*Vr); // optimum chord length (m/s)
if(idealchord < maxchord) chord = idealchord; // this limits the chord length if it exceeds the maximum
else chord = maxchord; // allowed length defined by the variable maxchord
solidity = chord*Nb/(2*pi*dimRi*R); // solidity is ratio of blade area to swept area
PTLfactor = (2/pi)*acos(exp((Nb*(dimRi - 1))/(2*dimRi*sin(phi)))); // Prandtl tip loss corrector factor
a = 1/((4*PTLfactor*sin(phi)*sin(phi))/(solidity*(Cl*cos(phi) + Cd*sin(phi))) + 1); // axial induction factor
aprime = 1/((4*PTLfactor*sin(phi)*cos(phi))/(solidity*(Cl*sin(phi) - Cd*cos(phi))) - 1); // tangential induction factor
diff1 = fabs(aguess -a); // returns absolute difference between aguess and a
diff2 = fabs(aprimeguess - aprime); // returns absolute difference between aprimeguess and aprime
/*printf("ag= %1.4f a= %1.4f apg= %1.4f ap= %1.4f Vr= %3.3f ch= %1.3f phi= %1.4f \n", aguess, a, aprimeguess, aprime, Vr, chord, phi); */ //prints variables
aguess = a; // redefines aguess for next iteration
aprimeguess =aprime; // redefines aprimeguess for next iteration
}
deltaQ = 0.5*rho*U*(1-a)*omega*dimRi*dimRi*R*R*(1+aprime)*Nb*(Cl*sin(phi)-Cd*cos(phi))*chord*delta_r/(sin(phi)*cos(phi)); // torque produced by element
if(deltaQ > opt_deltaQ) // this loop records the current variable values if they produce a higher deltaQ than any previous angle of attack
{
opt_deltaQ = deltaQ;
opt_alpharad = alpha_rad;
opt_phi = phi;
opt_chord = chord;
}
}
beta_rad= opt_phi - opt_alpharad; // setting angle beta, in radians
beta_deg= beta_rad*180/pi; // setting angle beta, in degrees
printf("\n\nopt deltaQ = %3.2f for beta= %2.2fdeg, alpha= %2.1fdeg, chord= %1.2fm at R= %1.2fm\n", opt_deltaQ, beta_deg, opt_alpharad*180/pi, opt_chord, dimRi*R);
}
return 0;
}