Multiple streamtubes model: Code, sample run, & output
Matlab code file “predictcpcurve_sm.m”
%Single multiple streamtubes model for calculating Cp curve
%>for fixed pitch straight bladed Darrieus wind turbine
%>with momentum equation between 0 < a < 0.4
%>and Glauert empirical formula between 0.4 < a < 1.0
%>Include Reynolds number effect
%>Not including dynamic stall airfoil data modification
%>NOTE: A_min = 0.49 so start at tsr 1.5 when data until 90deg only
%
%INPUT
%fileCtang = filename for tangential cof table from "makedatatable.m"
%fileCnorm = filename for normal coefficients data table
%sol = single value of rotor solidity Nc/D
%tsr = vector of tsr values
%uRe = vector of airfoil Re numbers at stationary condition
%>the vector of Reynolds number correspond to the tsr values
%
%OUTPUT
%tsr_Cp_Ct_meanA_minA_maxA_meanRe_minRe_maxRe_maxAlp =
%> table of [tsr Cp Ct meanA minA maxA meanRe minRe maxRe maxAlp]
function tsr_Cp_Ct_meanA_minA_maxA_meanRe_minRe_maxRe_maxAlp...
= predictcpcurve_sm(fileCtang, fileCnorm, sol, tsr, uRe);
tableCtang = load(fileCtang);
tableCnorm = load(fileCnorm);
global a_ACCURACY a_MIN a_MAX axisAlpha axisRe;
a_ACCURACY = 0.01; a_MIN = -0.49; a_MAX = 0.99;
axisAlpha = [0 : 0.1 : 0.1*(length(tableCtang(1,:))-1)];
axisRe = [0 : 10e3 : 10e3*(length(tableCtang(:,1))-1)];
delta_theta = 5; %all num. that can divide 90
theta = [delta_theta/2 : delta_theta : 180]';
%For every tsr
for itsr = 1 : 1 : length(tsr)
%For every theta position (streamtube)
for itheta = 1 : 1 : length(theta)
%Find induction factor and properties for this theta tube
%under current tsr speed
%NOTE: Crosswind force cancelled out due to symmetry up&downwind
[a, Ur, alpha, Re, CtubeThru, Ctang, Cnorm]...
= getTubeProps(tableCtang, tableCnorm, sol, uRe(itsr),...
tsr(itsr), theta(itheta));
%Insert every theta tube record in a constant tsr table
theta_a_Ur_alpha_Re_CtubeThru_Ctang_Cnorm(itheta,:) ...
= [theta(itheta) a Ur alpha Re/1e6 CtubeThru Ctang Cnorm];
end
%Print to screen some info of progress for some tsr
if rem(tsr(itsr), 1) == 0
PRINT_tsr... %Redundant variables
= tsr(itsr)
PRINT_theta_a_Ur_alpha_Re_Thr... %Matlab clips off long name
= theta_a_Ur_alpha_Re_CtubeThru_Ctang_Cnorm(2:3:end, 1:6)
%%skip some angles and columns
end
%Tag if any reverse far wake, no convergence,
%and mean induction factor in all tubes in the swept area
vectA = theta_a_Ur_alpha_Re_CtubeThru_Ctang_Cnorm(:, 2);
meanA = mean(vectA);
maxA = max(vectA);
minA = min(vectA);
%Tag Reynolds number for mean, max, and min values
vectRe = theta_a_Ur_alpha_Re_CtubeThru_Ctang_Cnorm(:, 5);
meanRe = mean(vectRe);
maxRe = max(vectRe);
minRe = min(vectRe);
%Tag angle of attack also
vectAlp = theta_a_Ur_alpha_Re_CtubeThru_Ctang_Cnorm(:, 4);
maxAlp = max(vectAlp);
%Find torque and power coefficient for this tsr
vecUr = theta_a_Ur_alpha_Re_CtubeThru_Ctang_Cnorm(:, 3);
vecCtang = theta_a_Ur_alpha_Re_CtubeThru_Ctang_Cnorm(:, 7);
n_theta = length(theta);
Ct = sol * sum(vecUr.^2 .* vecCtang) / n_theta;
Cp = Ct * tsr(itsr);
%Insert record for every tsr for this input solidity
tsr_Cp_Ct_meanA_minA_maxA_meanRe_minRe_maxRe_maxAlp(itsr,:)...
= [tsr(itsr) Cp Ct meanA minA maxA meanRe minRe maxRe maxAlp];
end
PRINT_output_table_format... %Redundant var for print to screen info
= 'tsr_Cp_Ct_meanA_minA_maxA_meanRe_minRe_maxRe_maxAlp'
%=====================================================================
%Function to find corrects properties for a tube at certain theta
%Trial and error method used first to find correct induction factor
%in the specified interval at certain accuracy
%Used the slow but sure converge bisection method
function [a, Ur, alpha, Re, CtubeThru, Ctang, Cnorm]...
= getTubeProps(tableCtang, tableCnorm, sol, uRe, tsr, theta);
global a_ACCURACY a_MIN a_MAX;
a_left = a_MIN;
a_right = a_MAX;
[Ur, alpha, Re, CtubeThru, Ctang, Cnorm, f_left] = ...
getTrialTubeProps(tableCtang, tableCnorm, sol, uRe, tsr,...
theta, a_left);
[Ur, alpha, Re, CtubeThru, Ctang, Cnorm, f_right] = ...
getTrialTubeProps(tableCtang, tableCnorm, sol, uRe, tsr,...
theta, a_right);
%Check if sure root exist
if f_left*f_right <= 0
%While interval not accurate and no lucky boundary
while (a_right-a_left)>a_ACCURACY & (f_left*f_right)~= 0
%Define new interval
a_center = (a_left + a_right) / 2;
[Ur, alpha, Re, CtubeThru, Ctang, Cnorm, f_center]...
= getTrialTubeProps(tableCtang, tableCnorm, sol, uRe, tsr,...
theta, a_center);
if f_left * f_center <= 0
a_right = a_center;
f_right = f_center;
else % > 0
a_left = a_center;
f_left = f_center;
end
end
%Find root value from final interval
if f_left == 0 %lucky boundary
a = a_left;
elseif f_right == 0
a = a_right;
else %(a_right-a_left)<accuracy
a = (a_left + a_right) / 2;
end
%else no sure root
else
%Give out of range root value
a = a_MAX*2;
end
%Return properties from correct induction factor for this tube
[Ur, alpha, Re, CtubeThru, Ctang, Cnorm, CthruDiff]...
= getTrialTubeProps(tableCtang, tableCnorm, sol, uRe, tsr, theta, a);
%=====================================================================
%Function to calculate a trial(might incorrect) tube properties
%for a given guessed induction factor
function [Ur, alpha, Re, CtubeThru, Ctang, Cnorm, CthruDiff]...
= getTrialTubeProps(tableCtang, tableCnorm, sol, uRe, tsr, theta, a);
[Ur, alpha, Re, Ctang, Cnorm] = ...
getFoilProps(tableCtang, tableCnorm, uRe, tsr, theta, a);
CtubeThru1 = getFoilThrust(sol, theta, Ur, Ctang, Cnorm);
CtubeThru2 = getWindThrust(a);
CtubeThru = CtubeThru1; %simply pick one
CthruDiff = CtubeThru1 - CtubeThru2;
%=====================================================================
%Get thrust coefficient acting ON the airfoil in TWICE entry into tube,
%for tangent and normal coefficients, normalised relative velocity,
%and airfoil area to swept area ratio => solidity + theta
function CtubeThru = getFoilThrust(sol, theta, Ur, Ctang, Cnorm);
CtubeThru = sol*(Ur)^2*2/pi*(-Cnorm - Ctang/tan(theta/180*pi));
%Error at theta = 0/180 degree as no actual swept area exist here
%=====================================================================
%Get thrust coefficient BY wind momentum or Glauert empirical formula
%for a given induction factor
%A crude straight line approximation for Glauert formula is used
%between 0.4 < a < 1.0, 0.96 < CtubeThru < 2.0
function CtubeThru = getWindThrust(a);
if a < 0.4
CtubeThru = 4 * a * (1 - a);
else
CtubeThru = 26 / 15 * a + 4 / 15;
end
%=====================================================================
%Function to calculate airfoil properties at certain theta position,
%tsr rotation, and induction factor
function [Ur, alpha, Re, Ctang, Cnorm]...
= getFoilProps(tableCtang, tableCnorm, uRe, tsr, theta, a);
[Ur, Un, Ut, alpha] = getVelocityDiagram(theta, tsr, a);
Re = Ur*uRe;
[Ctang, Cnorm]...
= getAirfoilCoefficients(tableCtang, tableCnorm, Re, alpha);
%=====================================================================
%Function to calculate normalised relative, normal, tangent velocity,
%and angle of attack for a given theta, tsr, and induction factor
%No a=1 and tsr=0 as Ur zero magnitude
%so alpha undefined (all 360degree??)
function [Ur, Un, Ut, alpha] = getVelocityDiagram(theta, tsr, a);
Ur = sqrt(((1-a)*sin(theta/180*pi))^2 ...
+ ((1-a)*cos(theta/180*pi) + tsr)^2);
Un = (1-a)*sin(theta/180*pi);
Ut = ((1-a)*cos(theta/180*pi) + tsr);
alpha = atan2(Un, Ut)/pi*180; %return between 0 to pi and 0 to -pi
%use special inverse tangent (4 quadrants inv. tangent) bc.
%normal inv tangent always lose info when --,-+,++(all conditions)
%complain at a=1, tsr=0 because no actual direction & zero mag. Ur
%=====================================================================
function [Ctang, Cnorm]...
= getAirfoilCoefficients(tableCtang, tableCnorm, Re, alpha);
%Get unsigned lift and drag coefficients
global axisRe axisAlpha;
Ctang = interp2(axisAlpha, axisRe, tableCtang, alpha, Re, '*linear');
Cnorm = interp2(axisAlpha, axisRe, tableCnorm, alpha, Re, '*linear');
%Assigns sign to coefficients (SIGN CONVECTION set earlier)
%CAREFUL: org foil data => both lift & drag sign relative to wind dir
if alpha < 0 %never happen in upwind multiple tube :( STUPID
Cnorm = -Cnorm;
end
A sample run of “predictcpcurve_sm.m” in Matlab
» tsr = [1.5:0.25:8]'
tsr =
1.5000
1.7500
2.0000
2.2500
2.5000
2.7500
3.0000
3.2500
3.5000
3.7500
4.0000
4.2500
4.5000
4.7500
5.0000
5.2500
5.5000
5.7500
6.0000
6.2500
6.5000
6.7500
7.0000
7.2500
7.5000
7.7500
8.0000
» ure = ones(size(tsr)) * 53e3
ure =
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
53000
» cp = predictcpcurve_sm('Ctang.A10.NACA12.table',...
'Cnorm.A10.NACA12.table', 0.1, tsr, ure)
PRINT_tsr =
2
PRINT_theta_a_Ur_alpha_Re_Thr =
7.5000 0.0332 2.9612 2.4424 0.1569 0.1225
22.5000 0.0679 2.8833 7.1064 0.1528 0.2548
37.5000 0.0910 2.7768 11.4946 0.1472 0.3387
52.5000 0.0910 2.6532 15.7713 0.1406 0.3355
67.5000 0.0795 2.5013 19.8776 0.1326 0.3007
82.5000 0.0737 2.3112 23.4137 0.1225 0.2667
97.5000 0.0621 2.0952 26.3468 0.1110 0.2290
112.5000 0.0505 1.8569 28.1894 0.0984 0.1868
127.5000 0.0390 1.6073 28.3170 0.0852 0.1436
142.5000 0.0274 1.3636 25.7332 0.0723 0.1051
157.5000 0.0216 1.1583 18.8587 0.0614 0.0811
172.5000 0.0159 1.0323 7.1482 0.0547 0.0683
PRINT_tsr =
3
PRINT_theta_a_Ur_alpha_Re_Thr =
7.5000 0.0448 3.9490 1.8093 0.2093 0.1810
22.5000 0.0968 3.8500 5.1508 0.2040 0.3566
37.5000 0.1546 3.7066 7.9810 0.1964 0.5219
52.5000 0.1777 3.5608 10.5563 0.1887 0.5898
67.5000 0.1893 3.3939 12.7494 0.1799 0.6069
82.5000 0.1546 3.2213 15.0815 0.1707 0.5298
97.5000 0.1315 3.0123 16.6099 0.1597 0.4648
112.5000 0.1141 2.7840 17.0960 0.1476 0.4048
127.5000 0.0968 2.5528 16.3017 0.1353 0.3535
142.5000 0.0910 2.3451 13.6485 0.1243 0.3260
157.5000 0.0795 2.1782 9.3073 0.1154 0.2971
172.5000 0.0274 2.0397 3.5683 0.1081 0.0977
PRINT_tsr =
4
PRINT_theta_a_Ur_alpha_Re_Thr =
7.5000 0.0621 4.9314 1.4225 0.2614 0.2433
22.5000 0.1315 4.8139 3.9590 0.2551 0.4490
37.5000 0.2009 4.6595 5.9931 0.2470 0.6420
52.5000 0.2587 4.4900 7.5267 0.2380 0.7675
67.5000 0.2991 4.3170 8.6263 0.2288 0.8410
82.5000 0.3049 4.1484 9.5624 0.2199 0.8435
97.5000 0.2876 3.9703 10.2475 0.2104 0.8135
112.5000 0.2587 3.7789 10.4422 0.2003 0.7662
127.5000 0.2240 3.5809 9.8999 0.1898 0.7021
142.5000 0.1835 3.3889 8.4339 0.1796 0.6002
157.5000 0.1084 3.1945 6.1317 0.1693 0.3914
172.5000 0.0274 3.0384 2.3946 0.1610 0.1097
PRINT_tsr =
5
PRINT_theta_a_Ur_alpha_Re_Thr =
7.5000 0.0852 5.9081 1.1580 0.3131 0.3109
22.5000 0.1604 5.7846 3.1841 0.3066 0.5464
37.5000 0.2471 5.6160 4.6812 0.2977 0.7436
52.5000 0.3165 5.4432 5.7175 0.2885 0.8689
67.5000 0.3685 5.2740 6.3511 0.2795 0.9252
82.5000 0.3916 5.1151 6.7719 0.2711 0.9492
97.5000 0.3859 4.9574 7.0551 0.2627 0.9509
112.5000 0.3570 4.7909 7.1233 0.2539 0.9186
127.5000 0.3049 4.6100 6.8702 0.2443 0.8402
142.5000 0.2240 4.4097 6.1498 0.2337 0.7007
157.5000 0.1373 4.2159 4.4915 0.2234 0.4773
172.5000 0.0332 4.0434 1.7884 0.2143 0.1267
PRINT_tsr =
6
PRINT_theta_a_Ur_alpha_Re_Thr =
7.5000 0.1084 6.8850 0.9686 0.3649 0.3851
22.5000 0.2009 6.7452 2.5986 0.3575 0.6383
37.5000 0.2934 6.5747 3.7515 0.3485 0.8346
52.5000 0.3801 6.3963 4.4099 0.3390 0.9370
67.5000 0.4263 6.2421 4.8707 0.3308 0.9990
82.5000 0.4437 6.0976 5.1899 0.3232 1.0302
97.5000 0.4379 5.9528 5.3719 0.3155 1.0313
112.5000 0.4148 5.8013 5.3478 0.3075 0.9905
127.5000 0.3570 5.6317 5.1974 0.2985 0.9239
142.5000 0.2702 5.4392 4.6849 0.2883 0.7929
157.5000 0.1662 5.2394 3.4916 0.2777 0.5502
172.5000 0.0332 5.0431 1.4339 0.2673 0.1373
PRINT_tsr =
7
PRINT_theta_a_Ur_alpha_Re_Thr =
7.5000 0.1315 7.8619 0.8262 0.4167 0.4656
22.5000 0.2355 7.7118 2.1740 0.4087 0.7231
37.5000 0.3454 7.5299 3.0337 0.3991 0.9078
52.5000 0.4263 7.3633 3.5437 0.3903 1.0111
67.5000 0.4668 7.2209 3.9118 0.3827 1.0808
82.5000 0.4841 7.0858 4.1391 0.3755 1.1120
97.5000 0.4841 6.9515 4.2193 0.3684 1.1029
112.5000 0.4610 6.8120 4.1921 0.3610 1.0645
127.5000 0.4148 6.6599 3.9976 0.3530 0.9819
142.5000 0.3165 6.4711 3.6867 0.3430 0.8654
157.5000 0.1893 6.2587 2.8413 0.3317 0.6116
172.5000 0.0332 6.0428 1.1966 0.3203 0.1396
PRINT_tsr =
8
PRINT_theta_a_Ur_alpha_Re_Thr =
7.5000 0.1662 8.8274 0.7064 0.4679 0.5484
22.5000 0.2702 8.6787 1.8440 0.4600 0.7917
37.5000 0.4032 8.4813 2.4551 0.4495 0.9606
52.5000 0.4668 8.3353 2.9090 0.4418 1.0790
67.5000 0.5073 8.2012 3.1820 0.4347 1.1427
82.5000 0.5246 8.0758 3.3458 0.4280 1.1707
97.5000 0.5188 7.9515 3.4396 0.4214 1.1744
112.5000 0.4957 7.8209 3.4153 0.4145 1.1362
127.5000 0.4495 7.6773 3.2614 0.4069 1.0519
142.5000 0.3570 7.5001 2.9919 0.3975 0.9222
157.5000 0.2066 7.2734 2.3923 0.3855 0.6645
172.5000 0.0332 7.0426 1.0267 0.3733 0.1338
PRINT_output_table_format =
tsr_Cp_Ct_meanA_minA_maxA_meanRe_minRe_maxRe_maxAlp
cp =
Columns 1 through 7
1.5000 0.0263 0.0176 0.0358 0.0043 0.0679 0.0874
1.7500 0.0383 0.0219 0.0446 0.0043 0.0852 0.0991
2.0000 0.0551 0.0275 0.0544 0.0043 0.1026 0.1113
2.2500 0.0762 0.0339 0.0664 0.0043 0.1199 0.1236
2.5000 0.1006 0.0403 0.0799 0.0043 0.1430 0.1362
2.7500 0.1325 0.0482 0.0954 -0.0015 0.1604 0.1488
3.0000 0.1685 0.0562 0.1138 -0.0015 0.1893 0.1616
3.2500 0.2100 0.0646 0.1339 -0.0015 0.2066 0.1744
3.5000 0.2646 0.0756 0.1561 -0.0073 0.2355 0.1873
3.7500 0.3247 0.0866 0.1766 -0.0073 0.2702 0.2003
4.0000 0.3429 0.0857 0.1967 -0.0130 0.3107 0.2133
4.2500 0.3538 0.0833 0.2155 -0.0130 0.3396 0.2264
4.5000 0.3540 0.0787 0.2285 -0.0188 0.3570 0.2395
4.7500 0.3478 0.0732 0.2404 -0.0188 0.3743 0.2525
5.0000 0.3364 0.0673 0.2514 -0.0246 0.3916 0.2656
5.2500 0.3198 0.0609 0.2635 -0.0304 0.4090 0.2787
5.5000 0.3043 0.0553 0.2742 -0.0362 0.4205 0.2918
5.7500 0.2952 0.0513 0.2850 -0.0420 0.4321 0.3049
6.0000 0.2879 0.0480 0.2954 -0.0420 0.4437 0.3180
6.2500 0.2763 0.0442 0.3057 -0.0477 0.4552 0.3312
6.5000 0.2601 0.0400 0.3158 -0.0535 0.4668 0.3443
6.7500 0.2407 0.0357 0.3255 -0.0593 0.4784 0.3574
7.0000 0.2164 0.0309 0.3346 -0.0709 0.4899 0.3705
7.2500 0.1871 0.0258 0.3446 -0.0766 0.4957 0.3837
7.5000 0.1548 0.0206 0.3536 -0.0824 0.5073 0.3968
7.7500 0.1180 0.0152 0.3621 -0.0882 0.5130 0.4099
8.0000 0.0767 0.0096 0.3708 -0.0940 0.5246 0.4230
Columns 8 through 10
0.0269 0.1313 40.6860
0.0401 0.1443 33.6779
0.0533 0.1572 28.6477
0.0666 0.1702 24.7672
0.0798 0.1831 21.7764
0.0928 0.1961 19.2697
0.1060 0.2093 17.0960
0.1192 0.2222 15.1699
0.1322 0.2352 13.3574
0.1454 0.2481 11.7295
0.1584 0.2614 10.4422
0.1716 0.2743 9.2365
0.1846 0.2872 8.4055
0.1978 0.3002 7.7365
0.2108 0.3131 7.1367
0.2237 0.3261 6.5960
0.2366 0.3390 6.1509
0.2496 0.3520 5.7435
0.2628 0.3649 5.3892
0.2758 0.3778 5.0637
0.2887 0.3908 4.7637
0.3017 0.4037 4.4864
0.3143 0.4167 4.2291
0.3273 0.4293 4.0120
0.3402 0.4423 3.8111
0.3531 0.4552 3.6018
0.3661 0.4679 3.4462
» save('So010.A10.Naca12.Re53e3.curve', 'cp', '-ascii')
» plot(cp(:,1), cp(:,2))
»
Plotted graph of the predicted Cp curve
|
