Tuulivoiman teknisen potentiaalin määrittäminen

2. ohjaus kurssilla FYSS481 Tuulienergia, kesä 2015

Contents

Vaihe 0. Valmistelut

  1. Luo työasemalle hakemisto c:\MyTemp\FYSS481\Ohjaus2
  2. Kopioi tehtävänannossa mainitut aineistot eli i) mittausten tulokset ja ii) tuuliturbiinien tekniset tiedot luomaasi hakemistoon
  3. Käynnistä MATLAB ja luo Editor-ikkunassa komennolla Ctrl-n uusi script-tiedosto, johon kokoamme ohjaukset aikana käytettävät koodin pätkät.

Vaihe 1. Tuuliturbiinin ominaisuudet

Aluksi määritellään valittujen tuuliturbiinien ominaisuudet:

ja lisäksi

Tehtävä 1.1. Hae valittujen tuuliturbiinien ominaisuudet annetuista teknisistä tiedoista ja tallenna ne alla oleviin muuttujiin.

% Tuuliturbiinien ominaisuudet: Bornay 6000
v0 = 3.5;  % Cut-in speed
v1 = 12.0; % Rated speed
v2 = 14.0; % Cut-off speed
Pn = 6.0;  % Rated or Nominal power
D  = 4.0;  % Rotor diameter
h1 = 9.0;  % Hub height 1
h2 = 12.0; % Hub height 2

Vaihe 2. Tuuliturbiinin tehokäyrä

Justus et al. (1976) ovat esittäneet, että tuuliturbiinin tehokäyrää $P_{wec}(v)$ voidaan kuvata funktiolla

$$P_{wec}(v) = 0,\ \mbox{kun}\ v \le v_0$$

$$P_{wec}(v) = A + Bv + C v^2,\ \mbox{kun}\ v_0 < v \le v_1$$

$$P_{wec}(v) = P_r,\ \mbox{kun}\ v_1 < v \le v_2$$

$$P_{wec}(v) = 0,\ \mbox{kun}\ v > v_2$$

missä $v$ on tuulen nopeus roottorin navan korkeudella, $P_r$ tuuliturbiinin nimellisteho ja kertoimet $A$, $B$ ja $C$ saadaan ratkaisemalla ehtojen (1)-(3) muodostama lineaarinen yhtälöryhmä:

$$(1) \quad A + Bv_0 + C v_0^2 = 0$$

$$(2) \quad A + Bv_1 + C v_1^2 = P_r$$

$$(3) \quad A + Bv_c + C v_c^2 = P_r (v_c/v_1)^3$$

missä $v_c = (v_0 + v_1)/2$.

Tehtävä 2.1. Määrittele matriisiyhtälö $A x = y$ ja ratkaise vektorin $x = (A, B, C)^{T}$ alkioiden arvot.

A = zeros(3,3);
y = zeros(3,1);

vc = 0.5*(v0 + v1);
A = [1 v0 v0^2; 1 v1 v1^2; 1 vc vc^2]
y = [0 Pn Pn*(vc/v1)^3]'
x = A\y

A = x(1), B = x(2), C = x(3)

A =

    1.0000    3.5000   12.2500
    1.0000   12.0000  144.0000
    1.0000    7.7500   60.0625


y =

         0
    6.0000
    1.6163


x =

    0.7470
   -0.4815
    0.0766


A =

    0.7470


B =

   -0.4815


C =

    0.0766

Tehtävä 2.2. Määrittele tuuliturbiinin tehokäyrä käyttäen yllä annettua funktiota ja sille saatuja kertoimia.

u = 0:0.1:20;
Pwec = zeros(size(u));
Pwec(v0 < u & u <= v1) = A + B*u(v0 < u & u <= v1) + ...
    C*u(v0 < u & u <= v1).^2;
Pwec(v1 < u & u <= v2) = Pn;

Tehtävä 2.3 Piirrä tehokäyrän kuvaaja ja vertaa sitä teknisissä tiedoissa annettuu käyrään. Vastaako määritetty funktio riittävän hyvin todellista tehokäyrää? Miltä osin funktio ei vastaa todellista tehokäyrää?

figure(1)
plot(u,Pwec)
axis([0 20 0 10])
grid
xlabel('Tuulen nopeus napakorkeudella')
ylabel('Teho, kW')

Vaihe 3. Tuulimittausten tulosten käsittely

Sovitetaan Weibull jakauman tuulen nopeus havaintoihin kussakin tuulen suunnassa. Käytettävä aineisto on tallennettu otos07.mat tiedostoon ja ensimmäinen tehtäväsi on ladata aineisto MATLABiin ja tallentaan taulukot eri nimille.

load otos07;
v = WindSpeed;
d = WindDir;

Tehtävä 3.1 Määritä tuulen suunnan suhteelliset frekvenssit kun tuulen suunnan luokat ovat $0,1,2,\ldots,16$.

Id = 0:16;         % Luokat
fd = hist(d,Id);   % Frekvenssit
pd = fd/sum(fd)    % Suhteelliset frekvenssit

pd =

  Columns 1 through 7

    0.0449    0.0557    0.0540    0.0335    0.0480    0.0229    0.0202

  Columns 8 through 14

    0.0298    0.0439    0.0979    0.1510    0.1214    0.0801    0.0422

  Columns 15 through 17

    0.0243    0.0397    0.0906

Tehtävä 3.2 Piirrä histogrammi tuulen suunnan suhteellisista frekvensseistä.

figure(2)
bar(Id,pd)
xlabel('Tuulen suunta')
ylabel('Suhteellinen frekvenssi')

Tehtävä 3.3 Muunna tuulen suunnan suhteellinen frekvenssijakauma ajalliseksi frekvenssijakaumaksi kertomalla suhteelliset frekvenssit vuotuisella tuntimäärällä. Piirrä histogrammi tuulen suunnan ajallisesta frekvenssijakaumasta.

Td = pd*8760

figure(3)
bar(Id,Td)
xlabel('Tuulen suunta')
ylabel('Ajallinen frekvenssi, h')

Td =

   1.0e+03 *

  Columns 1 through 7

    0.3931    0.4883    0.4730    0.2932    0.4204    0.2002    0.1769

  Columns 8 through 14

    0.2608    0.3849    0.8577    1.3226    1.0631    0.7013    0.3699

  Columns 15 through 17

    0.2131    0.3479    0.7936

Tehtävä 3.4 Määritä tuulen nopeuden suhteelliset frekvenssit koko havaintoaineistolle, kun osavälien leveys on 0.25 m/s.

Iv = 0:0.25:20;   % Luokat
k  = length(Iv);  % Luokkien lukumäärä
fv = hist(v,Iv);  % Frekvenssit
pv = fv/sum(fv)   % Suhteelliset frekvenssit

pv =

  Columns 1 through 7

    0.1111         0    0.1243         0    0.0964    0.0815         0

  Columns 8 through 14

    0.0874         0    0.0941         0    0.0872    0.0742         0

  Columns 15 through 21

    0.0607         0    0.0484         0    0.0369         0    0.0272

  Columns 22 through 28

         0    0.0195    0.0138         0    0.0097         0    0.0073

  Columns 29 through 35

         0    0.0053    0.0038         0    0.0030         0    0.0023

  Columns 36 through 42

         0    0.0018         0    0.0013    0.0009         0    0.0007

  Columns 43 through 49

         0    0.0004         0    0.0003    0.0002         0    0.0001

  Columns 50 through 56

         0    0.0001         0    0.0001         0    0.0001         0

  Columns 57 through 63

    0.0000    0.0000         0    0.0000         0    0.0000         0

  Columns 64 through 70

         0         0         0         0         0         0         0

  Columns 71 through 77

         0         0         0         0         0         0         0

  Columns 78 through 81

         0         0         0         0

Tehtävä 3.5 Piirrä histogrammi tuulen nopeuden suhteellisesta frekvenssijakaumasta.

figure(4)
bar(Iv,pv)
xlabel('Tuulen nopeus')
ylabel('Suhteellinen frekvenssi')

Vaihe 4. Weibull jakauman sovittaminen mittausten tuloksiin

Sovitetaan Weibull jakauma erikseen kultakin tuulen suunnalta tehtyjen mittausten tuloksiin. Suoritetaan sovitukset käyttäen MATLABin Weibull funktioita:

Tehtävä 4.1 Alla oleva koodi laskee Weibull jakauman parametrit erikseen kullekin tuulen suunnan sektorille. Hae Googlella hakusanoilla "mathworks matlab wblfit" funktion dokumentaatio ja selvitä mitä alla olevat koodinpätkät tekevät?

Huom. Sovituksessa otetaan huomioon anemometrin mittausalue, joka alkaa 1 m/s eli tätä pienemmät mittaustulokset sensuroidaan pois havainnoista.

c = zeros(1,16); % Taulukko skaalauskertoimille
cci = zeros(2,16);
k = zeros(1,16); % Taulukko muotokertoimille
kci = zeros(2,16);

figure(5)
u = 0:0.1:16;
for id = 1:16
    idx = find(v~=0.0 & d==id);
    [parmhat,parmci] = wblfit(v(idx),0.05,v(idx)<=1);
    c(1,id) = parmhat(1); cci(:,id) = parmci(:,1);
    k(1,id) = parmhat(2); kci(:,id) = parmci(:,2);
    %
    subplot(4,4,id)
    plot(u,wblpdf(u,parmhat(1),parmhat(2))*0.1)
    axis([0 16 0 0.08])
    title(int2str(id))
end

c_cci = [c' cci'], k_kci = [k' kci']

c_cci =

    3.9649    3.9384    3.9916
    3.7146    3.6926    3.7368
    3.2812    3.2533    3.3093
    2.6099    2.5925    2.6274
    2.1251    2.1036    2.1469
    2.0717    2.0494    2.0942
    2.2012    2.1802    2.2224
    2.4497    2.4307    2.4688
    2.9701    2.9552    2.9852
    3.6671    3.6523    3.6821
    4.0817    4.0637    4.0997
    4.4470    4.4223    4.4718
    3.3415    3.3143    3.3688
    3.0993    3.0648    3.1343
    3.3216    3.2951    3.3483
    3.8518    3.8329    3.8708


k_kci =

    2.1394    2.1189    2.1601
    2.4419    2.4177    2.4663
    2.3331    2.3022    2.3644
    2.6793    2.6487    2.7103
    3.3538    3.2858    3.4231
    3.5229    3.4439    3.6038
    3.1351    3.0792    3.1921
    2.8284    2.7900    2.8674
    2.4748    2.4547    2.4951
    2.2761    2.2621    2.2901
    2.2267    2.2121    2.2413
    2.1446    2.1271    2.1621
    2.2320    2.2056    2.2586
    2.2674    2.2307    2.3047
    2.2627    2.2359    2.2900
    2.2699    2.2527    2.2871

Tehtävä 4.2 Laske tuulen keskimääräinen nopeus kullekin sektorille ja vertaa saatuja arvoja MATLABin mean funktion antamiin arvoihin.

ka = zeros(1,16);
for id = 1:16
    idx = find(v~=0.0 & d==id);
    ka(1,id) = mean( v(idx) );
end
mu = c.*gamma(1 + 1./k);
mu_vs_ka = [mu' ka']

mu_vs_ka =

    3.5114    3.0590
    3.2941    2.8999
    2.9074    2.3122
    2.3203    1.7223
    1.9078    1.1048
    1.8646    1.0568
    1.9696    1.1383
    2.1822    1.3803
    2.6346    1.9418
    3.2484    2.6515
    3.6150    3.1556
    3.9383    3.4825
    2.9595    2.2937
    2.7454    1.9475
    2.9422    2.3613
    3.4119    3.0059

Vaihe 5. Weibull parametrien ekstrapolointi turbiinin napakorkeudelle

Käytetään Justuksen ja Mikhailin (1976) ehdottamaan menetelmään tuulen nopeuden arvojen ekstrapolointiin turbiinin napakorkeudelle. Ekstrapoloinninssa käytettävät yhtälöt ovat

$$k_2/k_1 = [1 - 0.0881 \ln(z_1/10)]/[1 - 0.0881 \ln(z_2/10)]$$

$$c_2/c_1 = (z_2 / z_1)^n$$

missä karheuseksponentti $n$ lasketaan kaavalla

$$n = (0.37 - 0.0881 \ln c_1)/[1 - 0.0881 \ln (z_1/10)]$$

Tehtävä 5.1 Määritä Weibull parametrit käyttäen edellä annettuja yhtälöitä. Anemometrin korkeutta merkitään yhtälöissä $z_1$:llä ja sille voimme käyttää arvoa 4 m.

z0 = 4; % Anemometrin korkeus

n = (0.37 - 0.0881*log(c))./(1 - 0.0881*log(z0/10)); % Karheuseksponentti

c1 = c.*(h1/z0).^n;
k1 = k.*(1 - 0.0881*log(z0/10))./(1 - 0.0881*log(h1/10));

c2 = c.*(h2/z0).^n;
k2 = k.*(1 - 0.0881*log(z0/10))./(1 - 0.0881*log(h2/10));


printmat([n' c1' k1' c2' k2'],'Weibull parameters by sector',...
    strcat(int2str(1:16)),'n c_1 k_1 c_2 k_2')

 
Weibull parameters by sector = 
                         n          c_1          k_1          c_2          k_2
            1      0.23007      4.77811      2.29085      5.10507      2.34986
            2      0.23539      4.49589      2.61472      4.81088      2.68207
            3      0.24550      4.00399      2.49824      4.29700      2.56259
            4      0.26416      3.23337      2.86898      3.48866      2.94288
            5      0.28091      2.66882      3.59117      2.89345      3.68367
            6      0.28299      2.60606      3.77232      2.82709      3.86949
            7      0.27804      2.75789      3.35706      2.98756      3.44353
            8      0.26932      3.04763      3.02862      3.29314      3.10663
            9      0.25362      3.64834      2.64998      3.92448      2.71824
           10      0.23644      4.44218      2.43719      4.75485      2.49996
           11      0.22771      4.90944      2.38428      5.24181      2.44569
           12      0.22072      5.31865      2.29636      5.66732      2.35551
           13      0.24402      4.07262      2.38995      4.36879      2.45152
           14      0.25015      3.79636      2.42787      4.07963      2.49041
           15      0.24450      4.05000      2.42291      4.34514      2.48532
           16      0.23243      4.65078      2.43053      4.97239      2.49314

KYSYMYS: Vastaavatko menetelmän antamat karheuseksponentin arvot ympäröivään maastoa? Eri karheuseksponentin arvoja vastaavat maastokuvaukset on annettu 1. luennon kalvossa 16.

Vaihe 6. Tuulen teho- ja energiatiheys

Lasketaan ensin tuulen teho- ja energiatiheys kullekin tuulen suunnalle seuraavilla yhtälöillä

$$P_{wind,i} = \frac{1}{2} \rho \int_0^{\infty} v^3 f_{W,i}(v) dv$$

$$E_{wind,i} = P_{wind,i} \times \tau_i$$

missä $f_{W,i}(v)$ on tuulen suunnalle $i$ määritetty tuulen nopeuden Weibull jakauma ja $\tau$ kertoo kuinka kauan vuoden aikana on tuullut kyseisestä suunnasta.

Edelleen vuotuinen tuulen energiatiheys saadaan kaavalla

$$E_{wind,tot} = \sum_{i=1}^{16} E_{wind,i} = \sum_{i=1}^{16} P_{wind,i} \times \tau_i$$

Tehtävä 6.1 Määritellään ensin funktio, joka laskee tuulen nopeuden kolmannen potenssin painotettuna tuulen nopeutta vastaavalla Weibull funktion arvolla:

fun = @(x,c,k) (x.^3).*wblpdf(x,c,k);

Tehtävä 6.2 Määritellään ilman tiheys:

rho = 1.225;                    % Air density, kg/m3

Tehtävä 6.3 Lasketaan tuulen teho- ja energiatiheys kullekin tuulen suunnalle ja kahdelle valitulle napakorkeudelle:

% 1. napakorkeus
for id = 1:16
    % Numeerinen integrointi sektoreittain
    E_x3(id) = integral(@(x)fun(x,c1(id),k1(id)),0,16);
end
P1_wind = 0.5*rho*E_x3;            % Wind power density by sector, W/m2
E1_wind = 1e-3*P1_wind.*Td(2:end); % Wind energy density by sector, kWh/m2

% 2. napakorkeus
for id = 1:16
    % Numeerinen integrointi sektoreittain
    E_x3(id) = integral(@(x)fun(x,c2(id),k2(id)),0,16);
end
P2_wind = 0.5*rho*E_x3;            % Wind power density by sector, W/m2
E2_wind = 1e-3*P2_wind.*Td(2:end); % Wind energy density by sector, kWh/m2

Tehtävä 6.4 Lasketaan tuulen vuotuinen energiatiheys:

E1_tot_wind = sum(E1_wind)         % Total wind energy density, kWh/m2
E2_tot_wind = sum(E2_wind)         % Total wind energy density, kWh/m2

E1_tot_wind =

  468.3022


E2_tot_wind =

  563.8065

Tehtävä 6.5 Tulostetaan lasketut arvot taulukkona:

printmat([P1_wind' P2_wind' E1_wind' E2_wind'],'Wind power density',...
    strcat(int2str(1:16)),'P_wind_1 P_wind_2 E_wind_1 E_wind_2')

 
Wind power density = 
                  P_wind_1     P_wind_2     E_wind_1     E_wind_2
            1     78.40397     93.76765     38.28549     45.78773
            2     59.64306     72.02214     28.21128     34.06660
            3     43.34009     52.70923     12.70864     15.45595
            4     21.12250     26.22399      8.87887     11.02329
            5     10.95895     13.87704      2.19373      2.77786
            6     10.08334     12.80252      1.78338      2.26431
            7     12.32906     15.55224      3.21579      4.05649
            8     17.26905     21.56762      6.64600      8.30031
            9     31.62440     38.81130     27.12539     33.28985
           10     60.17674     72.54765     79.59204     95.95427
           11     82.50708     98.63667     87.71697    104.86505
           12    107.92252    128.03104     75.68773     89.79015
           13     47.01866     57.01308     17.38994     21.08639
           14     37.66136     45.93743      8.02497      9.78846
           15     45.79112     55.57781     15.93234     19.33747
           16     69.18935     83.11641     54.90962     65.96232

Vaihe 7. Tuuliturbiinin energiantuotanto

Kun tunnemme tuulen nopeuden todennäköisyysjakauman voimme laskea tuuliturbiinin keskimääräisen tehon yhtälöllä

$$P_{ave} = \int_0^{\infty} P_{wec}(v) f(v) dv$$

missä $P_{wec}(v)$ on tuuliturbiinin tehokäyrä tuulen nopeuden funktiona turbiinin napakorkeudella ja $f(v)$ on tuulen nopeuden todennäköisyysjakauma.

Tehtävä 7.1 Määritellään ensin funktio, joka laskee tuuliturbiinin tehon painotettuna annettua tuulen nopeutta vastaavalla Weibull jakauman arvolla:

fun = @(x,v0,v1,v2,Pn,c,k) WECpower(x,v0,v1,v2,Pn).*wblpdf(x,c,k);

Tehtävä 7.2 Lasketaan tuuliturbiinin keskimääräinen teho kullekin tuulen suunnalle:

% 1. napakorkeus
for id = 1:16
    % Numeerinen integrointi sektoreittain
    P_ave(id) = integral(@(x)fun(x,v0,v1,v2,Pn,c1(id),k1(id)),0,16);
end
E1_wec = P_ave.*Td(2:end); % Energiantuotanto sektoreittain, kWh

% 2. napakorkeus
for id = 1:16
    % Numeerinen integrointi sektoreittain
    P_ave(id) = integral(@(x)fun(x,v0,v1,v2,Pn,c2(id),k2(id)),0,16);
end
E2_wec = P_ave.*Td(2:end); % Energiantuotanto sektoreittain, kWh

Tehtävä 7.3 Lasketaan turbiinin kapasiteettikerroin sektoreittain prosentteina:

CF1 = 100.0*E1_wec./(Pn*Td(2:end));
CF2 = 100.0*E2_wec./(Pn*Td(2:end));

Tehtävä 7.4 Lasketaan turbiinin vuotuinen energiantuotanto:

E1_tot_wec = sum(E1_wec)         % Total wind energy density, kWh/m2
E2_tot_wec = sum(E2_wec)         % Total wind energy density, kWh/m2

E1_tot_wec =

   1.8062e+03


E2_tot_wec =

   2.3185e+03

Tehtävä 7.5 Lasketaan turbiinin vuotuinen kapasiteettikerroin prosentteina:

CF1_ann = 100.0*E1_tot_wec/(Pn*8760)
CF2_ann = 100.0*E2_tot_wec/(Pn*8760)

CF1_ann =

    3.4365


CF2_ann =

    4.4112

Tehtävä 7.6 Tulostetaan lasketut arvot taulukkona:

A = pi*(D/2)^2 % Roottorin pinta-ala, m2

printmat([A*E1_wind' E1_wec' A*E2_wind' E2_wec'],'WEC energy',...
    strcat(int2str(1:16)),'E1_wind E2_wec E2_wind E2_wec')

A =

   12.5664

 
WEC energy = 
                   E1_wind       E2_wec      E2_wind       E2_wec
            1    481.10961    166.11858    575.38565    209.10102
            2    354.51336    105.35737    428.09355    138.00571
            3    159.70145     40.29672    194.22524     54.45512
            4    111.57517     12.48614    138.52272     19.84077
            5     27.56717      0.40754     34.90759      1.01362
            6     22.41068      0.19154     28.45419      0.54762
            7     40.41084      1.10753     50.97541      2.37281
            8     83.51615      6.25855    104.30476     10.79166
            9    340.86770     65.12393    418.33262     93.10524
           10   1000.18301    305.58824   1205.79691    396.72362
           11   1102.28393    384.17488   1317.77309    483.33894
           12    951.12006    358.96995   1128.33626    440.16474
           13    218.52840     59.21004    264.97939     78.72464
           14    100.84479     23.53859    123.00538     32.23828
           15    200.21166     53.07304    243.00184     70.91928
           16    690.01458    224.29681    828.90694    287.19565

Vaihe 8. Tallenna tulokset Excel-tiedostoon

Lopuksi avaa uusi Excel-tiedosto, johon luodaan taulukko ohjauksen aikana kerätyille tiedoille kummastakin tuuliturbiinista. Kun olet saanut laskettua halutut tiedot Bornay 6000 tuuliturbiinille, aloita alusta ja ota tällä kertaa tarkasteluun Lely Aircon 10 -tuuliturbiini.

Vaihe 9. Lopuksi

Vertaillaan Excel-tiedostoon kerättyjä tietoja valituista tuuliturbiineista. Mitä suosittelisit näiden tietojen perusteella? Mitä mieltä olet tuuliturbiinien soveltuvuudesta energiantuotantoon tällä paikalla?


Published with MATLAB® R2013a