Windenergieanlagen - aerodynamische Auslegung und
numerische Berechnung mittels CFD
Bachelorthesis
eingereicht von
Name: Tom Hammelstein
Matrikelnummer: 492891
Anschrift: Rotdornweg 7
Abgabetermin: 28.02.2013
Betreuer
Prof. Dr. Frank Kameier
Maschinenbau und Verfahrenstechnik
Josef-Gockeln-Str. 9
40474 Düsseldorf
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Erklärung Hiermit versichere ich, Tom Hammelstein, die vorliegende Bachelor-Thesis selbständig
verfasst und keine weiteren als die angegebenen Hilfsmittel und Quellen benutzt zu haben.
Dies ist die von der Fachhochschule Düsseldorf zu bewertende Version.
Ort, Datum ________________________ Unterschrift ____________________
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Inhaltsverzeichnis 1. Einleitung ......................................................................................................................................... 7
1.1 Vorwort ................................................................................................................................... 7
1.2 Zielsetzung ............................................................................................................................... 7
2. Windturbine – Theoretische Grundlagen ........................................................................................ 9
2.1 Physik der Windenergie........................................................................................................... 9
2.2 Theorie von Betz .................................................................................................................... 10
2.3 Froude-Rankinesches Theorem ............................................................................................. 12
2.4 Eulersche Strömungsmaschinenhauptgleichung .................................................................. 14
2.5 Blattelementmethode ........................................................................................................... 16
3. Strömungssimulation - Computational Fluid Dynamics ................................................................ 24
3.1 Kontinuitätsgleichung ............................................................................................................ 24
3.2 Impulserhaltung (Navier-Stokes-Gleichung) ......................................................................... 26
4. Modellierung ................................................................................................................................. 28
4.1 NACA-Profile .......................................................................................................................... 28
4.2 Geometrieaufbereitung in Inventor ...................................................................................... 29
4.3 Randbedingungen ................................................................................................................. 32
4.4 Simulationsraum ................................................................................................................... 32
4.5 Netzgenerierung .................................................................................................................... 34
5. 2D-Profilumströmung mit dem Programm XFOIL ......................................................................... 37
5.1 XFOIL-Ergebnisse ................................................................................................................... 37
6. Ergebnisse...................................................................................................................................... 39
6.1 Residuenstudie ...................................................................................................................... 40
6.2 Stromröhrenaufweitung nach Betz / Froude-Rankine .......................................................... 43
6.3 Vergleich der Drehmomente ................................................................................................. 49
6.5 Drallbestimmung nach EULER ............................................................................................... 57
6.6 Kräfte am Rotor ..................................................................................................................... 62
6.6 Anlagenkennlinien der Modelle ............................................................................................ 66
6.7 Das Q-Kriterium ..................................................................................................................... 68
7. Zusammenfassung ......................................................................................................................... 70
Literaturverzeichnis ............................................................................................................................... 71
Anhang .................................................................................................................................................. 72
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Nomenklatur
Abkürzungen
2D Zweidimensional
3D Dreidimensional
NACA National Advisory Committee for Aeronautics
SST Shear-Stress-Transport-Modell
XFOIL Programm eines 2D-Panelverfahrens
WEA Windenergieanlage
num Numerisch
analy Analytisch
Lateinische Zeichen
A [m2] Stützstelle, Querschnittsfläche
b [m] Profilbreite
B Stützstelle
ca Auftriebsbeiwert
caxial Geschwindigkeitskomponente in Rotorachsenrichtung
cm Momentenbeiwert
cP, cP Betz Rotorleistungsbeiwert, Rotorleistungsbeiwert nach BETZ
cp Druckbeiwert
cw Widerstandsbeiwert
c1,c2,c3 [m/s] Windgeschw. weit vor, am und weit hinter dem Rotor
c2u [m/s] Geschwindigkeitskomponente in Umfangsrichtung
d [m] Profildicke
dr [m] Breite eines Kreisringsegments
E [kg m2/s2] Kinetische Energie
[N] Feldkräfte
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
, [N] Auftriebskraft, Widerstandskraft
, [N] Tangentialkraft, Normalkraft
, [N] Umfangskraft, Schubkraft
i Laufindex
l [m] Profiltiefe, Länge der Profilsehne
M [Nm] Rotordrehmoment
[kg/s] Massenstrom
N Kreissegmente
n [1/min] Drehzahl
ncrit Critical Amplification Ratio
P [W] Leistung, Abgegebene Generatorleistung
pb [bar] Barometrischer Druck
p-2, p+2 [N/m2] Statischer Druck kurz vor / kurz hinter der Rotorebene
R [m] Radius Rotorspitze, Radius Kreiszylinder
r [m] Radius
rN [m] Nasenradius
Re Reynolds-Zahl
t [s] Zeit
Tu Turbulenzgrad
u [m/s] Umfangsgeschwindigkeit
w [m/s] Relativgeschwindigkeit
x, y, z [m] Ortskoordinate
xc [m] Abszisse der Skelettlinie / Wölbungslinie
xd [m] Rücklage der größten Dicke
xf [m] Rücklage der größtenWölbung
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
xo, xu [m] Abszisse der Profiloberseite / Profilunterseite
yc [m] Ordinate der Skelettlinie / Wölbungslinie
yo, yu [m] Ordinate der Profiloberseite / Profilunterseite
Y spezifische Stutzenarbeit
Z Blattanzahl
Griechische Zeichen
αA [°] Aerodynamischer Anstellwinkel
α1, a [°] Anströmwinkel weit vor und am Rotor
β [°] Abströmwinkel
ε Gleitzahl
η [N s/m2] Dynamische Viskosität
λ Schnelllaufzahl
ν [m2/s] Kinematische Viskosität
Kreiszahl
ρ [kg/m3] Dichte
Sonstige
∇ Nabla-Operator
Δ Laplace Operator
Partielle Ableitung nach der Zeit
` Schwankungsgröße
¯ Mittelwert; Parametrisierung
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
1. Einleitung
1.1 Vorwort
Die Motivation eine Thesis über regenerative Energien zu schreiben, ist ebenfalls tief mit
meinen ideologischen Werten, einen Lebensraum für die nachkommenden Generationen zu
gewährleisten, gekoppelt. Die eingehende Zerstörung des Klimas durch Verwendung
konventioneller, fossiler Energieträger und die abzusehende Ressourcenknappheit und
damit verbunden zukünftig steigender Preise, muss zwangläufig weltweit zum Umdenken in
der Energiepolitik führen.
In Deutschland wurde durch die Einführung des Erneuerbare Energie Gesetz und Subvention
der Einspeisevergütung von Erneuerbaren Energien ein großer Schritt in Richtung
Klimaschutz und Umweltschutz umgesetzt.
Es wurden in Deutschland im Jahr 2010 rund 17 % der gesamten Energiebereitstellung aus
erneuerbaren Energien bewerkstelligt, was einer Verdreifachung seit 1998 entspricht, wobei
die Windkraft, obwohl es ein schwaches Windjahr war, mit 6 % den höchsten Anteil daran
hatte [1].
1.2 Zielsetzung
Ausgangspunkt der Thesis war das Interesse eine gesamte WEA unter Verwendung von
Ansys CFX dreidimensional numerisch zu untersuchen und anschließend qualitativ und
quantitativ mit der Theorie zu vergleichen. Dafür wurde zuerst die Anwendung der Software
erlernt. Ebenso musste ein Auftriebsläufer in Inventor konstruiert und dieser in Ansys CFX
übertragen werden. Für die numerische Berechnung ist es notwendig, das Modell zu
vernetzen (Flächen die über Knotenpunkte verbunden sind) und sinnvolle Randbedingungen
auf das zu berechnende Strömungsfeld aufzubringen.
Ziele dieser Arbeit waren grundsätzlich, ein geeignetes Modell zu schaffen, welches die
Strömungsverhältnisse einer WEA realitätsnah abzubilden, um Theorien von Betz, Froude-
Rankine, Euler sowie die Blattelementemethode anwenden zu können.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Bei geeigneter Wahl des strömungsphysikalischen Modells ist es möglich, Aufschlüsse über
Leistungsausbeute, Kräfteverteilung, Strömungsverluste, Akustik und Topologie im Nachlauf
der WEA, zu erhalten.
Es wird bewusst ein großer Luv-Schnellläufer, also Dreiflügler mit hoher Drehzahl simuliert,
da sich diese Technik wegen geringer Schallemissionen und ruhigem Laufverhalten und in
großtechnischen Anlagen durchgesetzt hat. Außerdem liefern sie, durch ihre niedrigen
Drehmomenten und der hohen Drehzahl, optimale Bedingungen, um Generatoren zu
betreiben und sind wirtschaftlicher als z.B. Vertikalläufer [3]. Außerdem geht der Trend in
Industrieländern zu immer größeren Anlagen, da die Infrastruktur des Stromnetzes
flächendeckend vorhanden ist und der gewonnene Strom einer dezentralen Klein-WEA in
Konkurrenz mit dem aus Großwindenergieanlagen und Großkraftwerken steht. Ebenso
müssen weiterhin Normen hinsichtlich der Betriebssicherheit von Kleinwindenergieanlagen
geschaffen werden.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
2. Windturbine – Theoretische Grundlagen
2.1 Physik der Windenergie
Die Energie, die dem Wind entzogen werden kann, wird als kinetische Energie der Luftmasse
verstanden,
2
12
1cmE
•
= ,
(2.1)
die in einer bestimmten Zeit die Fläche durchströmt.
1cAdt
dxAm ρρ ==
•
.
(2.2)
Da dieser Luftmassenstrom selbst noch der Geschwindigkeit proportional ist, ergibt sich für
die Leistung (Energie pro Zeiteinheit):
3
1
2
12
1
2
1cAcmEPWind ρ===
•
.
(2.3)
Dies ist die theoretisch maximale Leistung, die dem Wind entzogen werden kann, wenn man
davon ausgeht, dass in der Rotorebene eine mittlere axiale Geschwindigkeit von An- und
Abströmung vorliegt. Die entnommene Leistung wächst also mit der Geschwindigkeit in der
dritten Potenz, wobei eine Verdoppelung der Geschwindigkeit eine Verachtfachung der
Leistung bedeutet. Da die Windgeschwindigkeit, durch die Höhenlage des Rotors im
Luftraum, nur bedingt beeinflusst werden kann, muss folglich die durchströmte Fläche
größer werden, um die Energieausbeute weiter zu steigern.
Durch den Rotor wird die kinetische Energie des Windes durch Abbremsen und Umlenkung
(Drall) in mechanische Energie umgewandelt. Stellt sich die Frage, in wie weit der Wind
abgebremst werden kann, ohne die nachkommenden Luftmassen zu beeinflussen, um ein
Optimum an Ausbeute zu erzielen? Ferner stellt sich die Frage, wie groß der Drall stromab
der Windturbine im Sinne der Eulerschen Strömungsmaschinenhauptgleichung ist.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
2.2 Theorie von Betz
Pioniere auf diesem Gebiet waren Betz [2] und Rankine [3], die in ihren Theorien idealisierte
Windkonverter betrachteten, um die maximale Leistungsentnahme zu bestimmen. Dabei
stellten sie fest, dass dies der Fall sei, wenn die ursprüngliche Windgeschwindigkeit hinter
dem Rotor auf ein Drittel abgebremst wird
,
wobei dann in der Rotorebene eine Geschwindigkeit von
23
vorherrscht.
Zusammengefasst beträgt dann die maximale theoretische verlustfreie Leistungsentnahme,
BetzBetz cpcAP3
12
1ρ=
,
(2.4)
wobei der Leistungsbeiwert 16/27 ~ 0,59 beträgt, siehe Abbildung 2.1.
Abbildung 2.1 - Betzkoeffizient cp in Abhängigkeit vom
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Geschwindigkeitsverhältnis vor und nach dem Rotor
Betz nimmt bei seiner Theorie eine homogene Windströmung an, die sich hinter der
Rotorebene verzögert und sich daher aus Kontinuitätsgründen aufweiten muss, siehe
Abbildung 2.2.
Abbildung 2.2 - Windgeschwindigkeitsreduktion mit einhergehender Stromröhrenaufweitung und
Druckverlauf vor und nach dem Rotor einer idealisierten WEA [3]
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
332211 AcAcAc ρρρ == .
(2.5)
Die Dichte ist, wegen der geringen Druckänderungen zwischen Ein- und Austritt der
Rotorebene, nahezu konstant.
Die dem Wind entnommene Leistung ist demnach:
)(2
1 2
3
2
1 ccmPE Wind −==•
.
(2.6)
Wie oben schon erwähnt, muss also durch Verzögerung der Luftmassen Leistung
entnommen werden. Dabei wird bei zu starker Verzögerung der Durchsatz zu gering (c3 = 0)
und bei keiner Verzögerung ist die Leistungsentnahme zu gering (c1 = c3).
Zwischen diesen beiden Auslegungen muss es eine optimale Leistungsentnahme geben, bei
dem der Durchsatz
2cAm ρ=•
(2.7)
in der Radebene bekannt ist. c2 ist die mittlere Geschwindigkeit zwischen Rotorein- und
austritt.
Gleichung 2.8 in Gl. 2.7 eingesetzt ergibt:
)(2
1 2
3
2
12 cccAPE Wind −== ρ .
(2.8)
2.3 Froude-Rankinesches Theorem
Aufbauend auf der Arbeit von Rankine bestimmt Froude die Windgeschwindigkeit c2 in
Abhängigkeit von c1 und c3. Dazu wird die Schubkraft F durch den Impulssatz ausgedrückt.
( − ) ( − ) . (2.9)
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Ebenso kann die Energie durch die Bernoulli-Gleichung, jeweils links und rechts von der
Rotorebene, aufgestellt werden.
² + #
$ + %& '(² + #'(
$ + %&) + * +,-.
- ,
(2.10)
/(² + #/(
$ + %&) .² + #.
$ + %& + * +,-.
-( . (2.11)
Die Einheit ist N m / kg ( Energie pro Masse), welche durch Multiplikation von ρ die Form in
Drücken in N/m² = Pa erhält.
Die Strömung sei stationär, inkompressibel und reibungsfrei. Vor und nach der Turbine
wirken keine äußeren Kräfte auf das Fluid ein, also können die Gleichungen 2.11 und 2.12 zu
$ ² + 0 $
)² + 0) und
(2.12)
$ 1² + 01 $
² + 0 . (2.13)
vereinfacht werde.
Die statischen Drücke weit vor und weit hinter der Rotorebene sind gleich dem
barometrischen Umgebungsdruck: 02 0 0. Durch Subtraktion der Gleichung (2.13)
von (2.12) ergibt sich somit:
$ ( − ) 0) − 01 . (2.14)
Die Schubkraft F kann auch über die Differenz des statischen Drucks ausgedrückt werden
(0) − 01) . (2.15)
Setzt man nun Gleichung 2.15 in 2.14 ein und das Ergebnis wiederum in 2.9, erhält man die
Geschwindigkeit c2 gemäß dem Rankine - Froude Theorem
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
12(1 + 3)(1 − 3) 2(1 − 3) → 2 (1+.)
2 . (2.16)
Durch Einsetzen von Gleichung 2.16 in Gl. 2.8 erhält man, die aus dem Wind enthaltene
maximale Leistung:
−
+=
•2
1
3
1
3
1 1*12
1³
2
1
c
c
c
ccAE ρ
.
(2.17)
Der Ausdruck in den eckigen Klammern wird durch den Betzkoeffizient cp ersetzt und besitzt
ein Maximum in einem bestimmten Geschwindigkeitsverhältnis.
2.4 Eulersche Strömungsmaschinenhauptgleichung
In Turbinen wird die innere Energie des strömenden Mediums in mechanische
Wellenleistung umgesetzt. Die Wellenleistung setzt sich dabei aus dem Drehmoment und
der Winkelgeschwindigkeit zusammen [6].
5 67 , (2.18)
wobei die Winkelgeschwindigkeit mit
7 28 , (2.19)
definiert ist.
Das Drehmoment bzw. der Drehimpuls ist mit
6 9 (2.20)
zu bestimmen.
Daraus folgt
6 (9: − 9:) . (2.21)
Mit der Winkelgeschwindigkeit
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
5 (79: − 79:) , (2.22)
oder anders ausgedrückt
5 (;: − ;:) . (2.23)
Die spezifische Stutzenarbeit Y P/m muss postiv sein, dann gilt für Turbinen
@ ;: − ;: (2.24)
als Eulersche Strömungsmaschinenhauptgleichung für eine Turbine.
Da die Umfangsgeschwindigkeiten vor und nach dem Rotor gleich sind u2=u1 und wegen der
drallfreien Anströmung : 0, muss die Abströmgeschwindigkeit cC negativ sein und sich
somit entgegen der Rotation ausbreiten (Abb. 2.3)
@ ;: − ;: . (2.25)
Der Drall der Abströmung ist also immer entgegen der Laufrichtung des Rotors gerichtet.
Geschwindigkeitsdreiecke stellen einen Zusammenhang zwischen dem Absolutsystem und
dem Relativsystem des Rotors dar, wobei die Windgeschwindigkeit c im Absolutsystem
gleich der Anströmgeschwindigkeit w im Relativsystem am Rotor mit dessen
Umfangsgeschwindigkeit u ist.
F + ; . (2.26)
0
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 2.3 - Geschwindigkeitsdreieck vor und hinter dem Rotor
Die Geschwindigkeitsdreiecke sind das grundlegende Konzept für die Beschreibung von
Strömungsbedingungen in der Blade-Element-Theory (Blattelementmethode). In ihr werden
der Wind und die Rotation der Flügel als vektorielle Größen betrachtet.
2.5 Blattelementmethode
Im Gegensatz zu Widerstandsläufer sind die modernen Auftriebsläufer von der Ausbeute her
wesentlich effizienter und erreichen mit cp = 0,50 annähernd den idealen Betzschen
Leistungsbeiwert. Dies hängt zum einen damit zusammen, dass auf die Tragflügelprofile
durch Anströmung des Körpers nicht nur eine Widerstandskraft W in Richtung der
Anströmung sondern auch eine senkrecht zu ihr gerichtete Auftriebskraft A resultiert [3].
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
GHIIIIJ Auftriebskraft GKIIIIJ Tangentialkraft
GLIIIIIJ Widerstandskraft αBau Bauwinkel
GMIIIIJ Schubkraft α Anströmwinkel
GNIIIIJ Umfangskraft αA Aerodynamischer Anstellwinkel
Abbildung 2.4 - Auftriebs- und Widerstandskraft in Abhängigkeit vom Anstellwinkel [7]
F O( + ;). (2.27)
Ist die Anströmgeschwindigkeit bestimmt, lässt sich auch die Reynolds-Zahl wie folgt
ermitteln:
QR F ∗ TU . (2.28)
Die Reynolds-Zahl ist eine wichtige dimensionslose Ähnlichkeitsgröße der Fluidmechanik und
charakterisiert die Strömungsform in Abhängigkeit von der Reibung des Fluid. Die kritische
Reynolds-Zahl ist eine Kenngröße, die den Umschlag von laminarer zur turbulenten
Strömung in der Grenzschicht angibt. Bei den großen Auftriebsläufern liegen die kritischen
Reynolds-Zahlen im Bereich von 0,5 105 bis 5 105. Bedingt durch die höhere Reynolds-Zahl
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
erfolgt hier der Übergang von der laminaren zur turbulenten Grenzschicht entweder ohne
Ablösung oder innerhalb einer laminaren Ablöseblase. Die Grenzschicht sollte in jedem Fall
größtenteils anliegend sein [9].
Praktisch sämtliche technisch eingesetzten Profile werden kritisch betrieben. Die
Flügelprofile sind Gegenstand der Forschung und ihre charakteristischen Kurven geben
dimensionslosen Auftriebsbeiwert cV und Widerstandsbeiwerte cW über dem Anstellwinkel
αA an (Abb. 2.5).
Abbildung 2.5 - Charakteristik eines Tragflügels in Abhängigkeit vom Anstellwinkel [3]
Wird das Profil vom Flügel variiert, so entstehen unterschiedliche Kurven, die im Windkanal
empirisch oder aber durch Simulationen näherungsweise ermittelt werden können. Welches
Profil nun optimal ist, hängt von dem Einsatzgebiet des Flügels ab. In der Regel werden die
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
herrschenden Druckbedingungen am Tragflügel und damit auch der Auftrieb über gemittelte
Werte berechnet. Der sogenannte drag/lift-Coefficient oder auch Gleitzahl genannt, stellt
das reziproke Verhältnis aus Widerstandsbeiwert und Auftriebsbeiwert in Abhängigkeit vom
aerodynamischen Anstellwinkel dar.
X (Y)(Y)
(2.29)
Die dimensionslosen Auftriebsbeiwerte cV und Widerstandsbeiwerte cW für die analytische
Auswertung werden in dieser Arbeit mit XFOIL bestimmt, siehe Kapitel 5.
Die Auftriebskraft an einem bestimmten Flügelschnitt wird aus der Tragflügeltheorie
bestimmt:
IIIIJ 2FT+9Z(Y), (2.30)
IIIIIJ 2FT+9\(Y), (2.31)
IIIJ IIIIJ],(Y) ↔ IIIIIJ,_8(Y). (2.32)
Nicht nur die Auftriebskraft beeinflusst die Schubkraft, auch die Widerstandskraft muss man
einfließen lassen. Da beide Kräfte in die gleiche Richtung weisen, werden diese dann addiert.
IIIJ IIIIJ],(Y) + IIIIIJ,_8(Y). (2.33)
Die Umfangskraft wird ähnlich ermittelt,
:IIIJ \IIIIJ],(Y) ↔ IIIIJ,_8(Y), (2.34)
mit dem Unterschied, dass die Widerstandskraft entgegen der Umfangskraft gerichtet ist.
:IIIJ IIIIJ,_8(Y) −\IIIIJ],(Y). (2.35)
6(9) :IIIJ9 . (2.36)
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Die Kräfte und Momente am gesamten Rotor ergeben sich dann aus der Summe der
einzelnen Kräfte aller Flügelschnitte, siehe Abb. 2.6:
`+(9)a
, (2.37)
b`+(9)a
, (2.38)
6 b`+(9)a
9, (2.39)
5 67 . (2.40)
Abbildung 2.6 - Luftkräfte am Flügelschnitt [3]
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 2.7 zeigt die gesamten Geschwindigkeitsdreiecke an einem Tragflügelsegment.
Das Geschwindigkeitsdreieck (rote Linien), bestehend aus c und u, aus denen sich w und
α bestimmen lassen, liegt weit vor dem Rotor. Das blau-farbige Dreieck ist folglich hinter
dem Rotor und wird unter der Annahme einer axialen Geschwindigkeitsverzögerung von 2/3
gebildet. Die grüne gestrichelte Linie bildet das fiktive Strömungsdreieck in der Rotorebene,
wobei sich die relative Anströmgeschwindigkeit aus
F F],(Y − Y). (2.41)
bestimmen lässt.
Aus den gegebenen Geschwindigkeitsvektoren unter Bildung eines rechtwinkligen Dreiecks
lässt sich nun der Anströmwinkel vor dem Rotor, aus
de8Y : , (2.42)
berechnen. Aus den theoretischen Anströmwinkeln vor dem Rotor resultiert der
Anströmwinkel nach Betz in der Rotorebene.
Y Y . (2.43)
Durch die in der Rotorebene entstehende Verzögerung verringert sich die Geschwindigkeit
c1,axial , folglich wird der resultierende Anströmwinkel α ebenfalls kleiner. Betz ging davon
aus, dass die Strömung von der Geschwindigkeit c1 weit vor dem Rad auf c3,axial hinter der
Rotorebene verzögert wird, ohne eine axiale Richtungsänderung einzubeziehen. Schmitz
jedoch berücksichtigt auch die Drallkomponente Δu in Umfangsrichtung in der
Umfangsebene und stromab des Rotors. Diese entsteht durch die Umlenkung der Strömung
aufgrund der Geometrie und ist in der Strömungsmaschinentheorie für den Arbeitsumsatz
verantwortlich. Nach Schmitz sind aber nur 50% davon nutzbar, so dass die Drallkomponente
Δu zu fC (Nomenklatur Schmitz) oder zu cC/2 (Nomenklatur Strömungsmaschinenbauer)
wird. Der Drall einer WEA ist immer entgegen der Umfangsgeschwindigkeit u gerichtet, siehe
Abbildung 2.7.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 2.7 – Geschwindigkeitsdreiecke von Schmitz und Betz [3]
Mit Änderung der Windgeschwindigkeit und der daraus resultierenden Nenndrehzahl bilden
sich neue Geschwindigkeitsdreiecke. Abbildung 2.8 veranschaulicht die Strömungs-
verhältnisse vor dem Rotor an den einzelnen Profilen im Rotorblattquerschnitt in einem
Betriebspunkt, die für die Auslegung einer Anlage relevant sind.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 2.8 - Geschwindigkeitsdreiecke vor dem Rotor im jeweiligen Profilschnitt
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
3. Strömungssimulation - Computational Fluid Dynamics
Durch die wachsenden Speicherkapazitäten und Rechnerleistungen ist es möglich,
zunehmend komplexere Strömungsprobleme unter vertretbarem Aufwand, akzeptabler Zeit
und ausreichender Genauigkeit und Sicherheit zu lösen. Die CFD wird dadurch eine immer
wichtigere Ergänzung der experimentellen Strömungsmechanik. Numerische Verfahren sind
Methoden und Algorithmen zum näherungsweisen Lösen von Modellgleichungen mittels
einfacher arithmetischer Operationen an diskreten Stellen des Lösungsgebietes, das in
Gitterstrukturen, in sog. finite Elemente aufgeteilt wird. Die mathematischen Näherungen
(Approximationen) sind durch physikalische Erhaltungsbedingungen geprägt und dienen zur
Simulation. In der Strömungslehre werden dazu vier Grundgleichungen verwendet [9]:
• die Bilanzgleichung für die Masse (Kontinuitätsgleichung): die zeitliche
Änderung der Masse eines materiellen Volumens ist null.
• die Bilanzgleichung für den Impuls (Navier-Stokes-Gleichung): die zeitliche
Änderung des Impulses eines materiellen Volumens ist gleich der am Volumen
angreifenden äußeren Kraft.
• die Bilanzgleichung für den Drehimpuls (Drehimpulssatz): die zeitliche
Änderung des Drehimpulses eines materiellen Volumens ist gleich dem am
Volumen angreifenden Drehmoment.
• die Bilanzgleichung für die Energie (Energiesatz): die zeitliche Änderung der
inneren und der kinetischen Energie eines materiellen Volumens ist gleich der
von außen angreifenden Leistung (unter Vernachlässigung von Strahlung und
innerer Dissipation). Die Energiebilanz ist relevant für kompressible
Strömungen, die im Rahmen der vorliegenden Arbeit nicht betrachtet werden
müssen.
3.1 Kontinuitätsgleichung
Nach dem Massenerhaltungssatz müssen in jeder Stromröhre folgende Bedingungen
herrschen:
Kontinuitätsgleichung - Massenbilanz
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
(3.1)
(3.2)
In diesem Fall handelt es sich um eine eindimensionale kompressible Kontinuitätsgleichung,
da die Bedingungen aus Abbildung 3.1 gelten.
Abbildung 3.1 - Einteilung von Strömungsmaschinen [8]
Wenn die Dichte g]8,d ist, lässt die Gleichung sich zu
(3.3)
vereinfachen.
Die Kontinuitätsgleichung für mehrdimensionale inkompressible und instationäre
Strömungen gestaltet sich schwieriger.
Dazu betrachten wir einen Raum, indem das einströmende Volumen hijk gleich groß dem
ausströmenden Volumen h:- ist, also
hijk − h:- 0. (3.4)
Dabei ist das einströmende Volumen in X-Achsenrichtung
l +d+m+& , (3.5)
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
gleich dem ausströmenden Volumen
l n l +o+d+m+&. (3.6)
, wobei
n l . (3.7)
die Änderung der Geschwindigkeit in x-Achsenrichtung ist.
Die Differenz der beiden lautet dann
n l +o+m+&+d 0. (3.8)
Dies lässt sich nun auf die y und z-Achse übertragen und lautet in differentieller Form:
pnpq +
prps + pt
pu 0 oder in vektorieller Schreibweise +_U ] .
3.2 Impulserhaltung (Navier-Stokes-Gleichung)
Impuls wird durch drei Wege transportiert: Durch Konvektion, Diffusion und zusätzlich über
Druck. Die Diffusion- und der Druckterm sind in den Spannungen enthalten. Die
Schubspannungen haben nur diffusive Anteile (Reibung). Die Normalspannungen haben
einen Druckanteil und einen Diffusionsanteil (Reibungs- bzw. Viskositätsanteil).
Die diffusiven Flüsse sind die Spannungen mit negativen Vorzeichen. Das negative
Vorzeichen ist selbstverständlich, weil der diffusive Impulsfluss in die Richtung stattfindet, in
der der Impuls (Geschwindigkeit, Druck) abnimmt, d.h. die Gradienten von der
Geschwindigkeit bzw. des Druckes negativ sind.
Der Impuls ist eine vektorielle Größe, d.h. er hat eine Richtung. Deshalb müssen
Impulsbilanzen in allen Koordinatenrichtungen erstellt werden. Auf die Herleitung der
Gleichung wird an dieser Stelle verzichtet. Die NS-Gleichung in vektorieller Form lautet:
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
vv − w0 + xy. (3.9)
η Dynamische Viskosität ∇ Nabla-Operator
Feldkräfte Δ Laplace-Operstor
Dabei beschreibt der erste Term die Impulsänderung in drei Koordinatenrichtungen, da der
Geschwindigkeitsvektor drei Geschwindigkeitskomponenten l, z, in Abhängigkeit der
Zeit beschreibt. Er beinhaltet eine konvektive Beschleunigung %9e+w, die eine
Geschwindigkeitsänderung von Ort zu Ort beschreibt und eine lokale Beschleunigung , die
zeitabhängig ist. Der zweite Term drückt die dreidimensionalen Feldkräfte. Die letzten
beiden Terme beschreiben die Oberflächenkräfte/Spannungen, die sich aufgrund von Druck-
und Zähigkeitskräften einstellen. Ist der Reibungsterm xy vernachlässigbar, geht die NS in
die Euler-Gleichung über.
Wird die NS durch die Dichte dividiert, lautet die Gleichung:
| + | l~ j −
$ # l| + | l~². (3.10)
Dabei wird die kinematische Viskosität ν aus dem Verhältnis zwischen dynamischer
Viskosität η und der Dichte ρ gebildet. Setzt man die Indizes _ o;8+ o, m, & ein, so
erhält man die NS in x-Koordinatenrichtung. Analog lässt sich die Gleichung ebenfalls für die
y, z-Richtung aufstellen. Die Schwerkraft wird dabei häufig in technischen Strömungen
vernachlässigt.Zur Berechnung dreidimensionaler inkompressibler Strömungen wird die NS mit der
Kontinuitätsgleichung kombiniert, wobei die NS aus drei Gleichungen (für die x, y, z-
Koordinate) besteht. Es stehen also vier Gleichungen zur Verfügung, mit denen die vier
Unbekannten (Geschwindigkeit in x, y, z-Richtung und der Druck) näherungsweise berechnet
werden können, wenn die Randbedingungen am Ein- und Austritt bekannt sind. Bei
transienten Simulationen sind außer den Randbedingungen noch Anfangs- bzw.
Startbedingungen von Nöten
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
4. Modellierung
4.1 NACA-Profile
Für die Beschreibung der Profilgeometrie werden in dieser Arbeit, die in Abbildung 4.1
dargestellten Bezeichnungen nach Erich Hau [7] verwendet. Die Profilkoordinaten der Kontur
in Tabellenform wurden aus Profilgenerator von Jens Trapp und Robert Zores entnommen
[5]. Die NACA Profile sind mit mehrstelligen Kennziffern beschrieben, die Hinweise auf die
Profilgeometrie geben.
l Profiltiefe, Länge der Profilsehne yc Ordinate der Skelettlinie
x0 Abszisse der Profiloberseite (Saugseite) yo Ordinate der Profiloberseite
xc Abszisse der Skelettlinie yu Ordinate der Profilunterseite
xu Abszisse der Profilunterseite (Druckseite) d größte Profildicke
xd Rücklage der größten Dicke d f Wölbung
xf Rücklage der Wölbung f m Wölbungsmaß m =f/l
θ Steigungswinkel der Skelettlinie t Dickenverhältnis t=d/l
rN Nasenradius
Abbildung 4.1 – Profilgeometrie mit Berechnungsparametern
Die NACA-Profile mit vier Ziffern (4-digit) lassen sich wie folgt beschreiben:
1. Ziffer: Maximale Wölbung der Profiltiefe [%]
2. Ziffer: Wölbungsrücklage in Zehntel der Tiefe
3./4. Ziffer: Maximale Dicke in Prozent der Tiefe
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Bspw. besitzt das Profil im Blattspitzenbereich des Rotors NACA4415 die Eigenschaften:
- 4% Wölbung bei 40% der Tiefe mit einer
- maximalen Dicke von 12% der Tiefe,
- wobei die Dickenrücklage für alle vierziffrigen Profile bei 30% liegt.
Eine Übersicht der verwendeten Profile für die Geometrie des Rotors.
Profil αBau [°] Länge [m] R [m]
NACA 4415 0 0,72 27
NACA 4420 2 1,51 17
NACA 4425 4 1,97 12
NACA 4450 8 2,92 5
Tabelle 4.1 – Übersicht der verwendeten Profile und deren Eigenschaften
4.2 Geometrieaufbereitung in Inventor
Die Modellbildung des Flügelprofils erfolgte mit Autodesk Inventor 2011. Bei der
Konstruktion muss auf einige Schritte in der Modellbildung eingegangen werden, die enorme
Wichtigkeit zur Durchführung der Simulation in Ansys CFD beitragen. Die verwendeten
NACA-Profile 4415, 4420, 4425 und 4450 wurden als 2D-Koordinatenpunkte aus den
Profilkatalogen von Eppler entnommen und in eine Excel-Tabelle übertragen. Diese wird in
Inventor in eine 2D-Skizze importiert und mit Splines innerhalb der Skizze in Inventor
miteinander verbunden, um die Profilkontur zu generieren. Die Approximation an die
eigentliche Profilkontur ist dabei von großer Bedeutung, denn sie muss krümmungsstetig
sein, um Druck- und Geschwindigkeitssprünge entlang der Profilkontur in der numerischen
Berechnung zu vermeiden. Des Weiteren wurde darauf geachtet, dass die Profilhinterkante
nicht in einem zu spitzen Winkel zuläuft. Ansys CFD produziert Fehler bei der automatischen
Erstellung der Netze, wenn ein Winkel innerhalb einer Geometrie zu spitz zuläuft. Um den
Fehler zu vermeiden, wurde der letzte Punkt an der Abströmkante entfernt und eine Gerade
zwischen Profilunter- und Profiloberseite gezogen (Abbildung 4.2).
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 4.2 - Abgeschnittene Profilhinterkante des NACA Profils
Die Profilnasen der verwendeten Profile liegen dabei radial in einer Flucht, siehe Abbildung
4.3. Mit der Funktion Erhebung entstehen nun die Flächen des Rotors. Dabei dürfen die
Oberflächen nicht zu grob gestaltet sein. Mit der automatischen Erstellung werden die
Koordinatenpunkte beliebig radial zwischen den einzelnen Profilen verbunden, was zu einer
unebenen Rotoroberfläche führt.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 4.3 - Anordnung der einzelnen Profile im Rotorblatt
In Inventor kann die Erhebung manuell bearbeitet werden, was eine genauere Verbindung
der Koordinatenpunkte innerhalb der Skizzen zulässt, siehe Abbildung 4.4.
Abbildung 4.4 - Radiale Verknüpfung der Profile mit Splines
Ebenso sollte man möglichst viele Splines in radialer Richtung einbauen, um eine homogene
Oberflächen zu schaffen. Durch die Tatsache, dass lediglich ein 120° Ausschnitt
(Drittelschnitt des gesamten Rotors) simuliert wurde, um die Rechenzeit und –aufwand zu
verkürzen, wurde die Nabe mit Flansch für die Rotorblätter stark vereinfacht und es wurde
auf eine Pitchregelung verzichtet, genauso wie auf eine Gondel (Maschinenraum für
Getriebe / Generator).
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
4.3 Randbedingungen
In der Simulation ist das Shear Stress Transport-Modell verwendet worden. Dabei
kombiniert das SST-Modell die Turbulenzmodelle k-ω-Modell und k-ε-Modell. Das k-ω-
Modell wird in wandnähe und das k-ε-Modell im Bereich der Hauptströmung verwendet
[11]. Auf die Transportgleichungen des SST-Modells wird hier nicht eingegangen und sei auf
das CFX-Handbuch verwiesen.
4.4 Simulationsraum
Der Simulationsraum (Abb. 4.5) ist ein 120° Kreisausschnitt über eine definierte Länge vor
und hinter dem Rotor. Dieser Raum wurde in drei Teile unterteilt. Der erste Teil wird in den
Simulationen nur für stationäre Strömungen verwendet und bildet den Raum vor der WEA,
indem die drallfreie Anströmung erfolgt. Dabei werden dem Inlet (durchströmte Fläche am
Eintritt des Raumes) drei Windgeschwindigkeiten (5m/s, 10m/s und 12m/s als
Betriebspunktvariation) aufgeprägt. Der zweite Teil ist der Rotationsraum, indem das zu
untersuchende Modell je nach Betriebspunkt mit unterschiedlichen Drehzahlen rotiert. Der
Rotor wird als Wand (Wall – umströmtes Volumen) definiert, der um die z-Achse rotiert.
Abbildung 4.5 - Schematischer Aufbau des Simulationsraumes
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Die Rotoroberfläche ist als hydraulisch glatt definiert. Das Rotorblatt ist dabei in y-
Achsenrichtung ausgerichtet und dreht in x-Achsen Richtung mit der Winkelgeschwindigkeit
Ω, siehe Abbildung 4.6
Abbildung 4.6 - Schematische Darstellung der WEA im Zusammenhang mit dem Kartesischen
Koordinatensystem
Der Abstand von Blattflügelspitze bis zum äußeren Rand des Simulationsraumes soll
mindestens das zweieinhalb-fache des Rotorradius betragen und der Nachlauf in etwa das
fünffache [10]. Der Nachlauf der WEA bildet den dritten Teil. Dabei sind die stationären
Räume 100 Meter lang, der instationäre Berechnungsraum wurde verkleinert auf 75 Meter
mit einem Radius von 100m, um Elemente einzusparen. Der Rotor hat dabei einen Abstand
von 114 Metern zum Inlet (Eintritt) bzw. 161 Meter bis zum Outlet (Austritt).
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
4.5 Netzgenerierung
Die Netzgenerierung ermöglicht das Modell in viele verschiedene einzelne Elemente zu
unterteilen (Finite Elemente). Man unterscheidet die Elemente in ihrer Geometrie: bei der
Verwendung von Hexaedern, also Würfel oder Quader, spricht man von strukturierten
Netzen; bei Tetraedern (Pyramiden) hingegen wird das Netz als unstrukturiert bezeichnet.
Die Größe, Anzahl sowie die Form dieser einzelnen Elemente haben großen Einfluss auf die
Genauigkeit der Ergebnisse. Es muss ein Kompromiss zwischen ausreichender Simulations-
genauigkeit und Rechenaufwand gewählt werden. Um wandnahe Effekte wie die
Grenzschichtströmung abbilden zu können, müssen umströmte Körper mit einer
Prismenschicht, einer Verfeinerung des Netzes hin zu Wänden, belegt werden. Es wurde
zwar versucht, ein fein strukturiertes Hexadernetz mit Prismenschicht in ICEM zu erstellen,
jedoch ist durch die Verwindung der Profile und gleichzeitige Verjüngung des Rotors eine
strukturierte Vernetzung nicht ohne erheblichen Zeitaufwand möglich gewesen. Lediglich die
stationären Bereiche bestehen aus Hexaedern und wurden mittels Verfeinerungsfunktion
orthogonal angeordnet. Im rotierenden Bereich wurde auf den Vernetzungsalgorithmus von
Ansys zurückgegriffen. Der Rotationsraum wurde mit einem flexiblen Tetraedernetz
umgesetzt und wird mit abnehmendem Abstand zum Profil immer feiner, wobei die
Elementgröße auf der Fläche des Rotors nur noch 80/40mm (Modell I/ Modell II) beträgt.
Modell I (λ=7,92)
Domain Knoten Elemente Volumen [m³]
Stat. Bereich I 155907 147750 20943,95
Stat. Bereich II 155958 147950 20943,95
Rot. Bereich 353795 2012212 15707,96
Gesamt 665660 2307912 57595,86
Tabelle 4.2 - Netzstatistik von Modell I
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Modell II (λ=7,07)
Domain Knoten Elemente Volumen [m³]
Stat. Bereich I 270351 262500 20943,95
Stat. Bereich II 270351 147950 20943,95
Rot. Bereich 621107 3570675 15707,96
Gesamt 1161809 4095675 57595,86
Tabelle 4.3 - Netzstatistik von Modell II
Beim Übergang vom Hexaeder- zum Tetraedernetz in z-Achsenrichtung liegen die
Knotenpunkte nicht übereinander und es muss interpoliert werden (Abbildung 4.8), weshalb
Ansys dort spezielle Schnittstellen, das Generalised Grid Interface, verlangt. Das Generalised
Grid Interface (GGI) erlaubt das beliebige Zusammenfügen von beliebigen Teilgittern. Weder
die Knotenverteilung noch die Gitterstruktur müssen an der gemeinsamen Schnittfläche
übereinstimmen. Durch das Generalised Grid Interface wird die Netzgenerierung einfacher
und schneller. Das gesamte Rechengebiet kann in einfachere Teilgebiete zerlegt werden, die
separat mit der jeweils optimalen Topologie vernetzt werden können [10].
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 4.7 - Schematische Darstellung des Simulationsraumes, Modell I
Die Frozen Rotor Methode ist ebenfalls ein stationäres Modell. Im Gegensatz zum Stage
Interface findet hier aber keine Umfangsmittelung der Strömungsgrößen statt, sondern nur
ein Wechsel in das jeweilige Bezugssystem. Das Frozen Rotor Interface kommt dann zum
Einsatz, wenn die Variation der Strömung in Umfangsrichtung sehr stark ist[10].
Abbildung 4.8 – Verwendung des GGI-Interface bei verschiedener Netzstruktur [Ansys]
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
5. 2D-Profilumströmung mit dem Programm XFOIL
Das Programm XFOIL wurde am Massachusetts Institute of Technology (MIT) von Mark Drela
[4] in Fortran entwickelt und steht unter der GNU General Public License (GPL) zur freien
Verfügung. Es verwendet ein 2D-Panelverfahren zur Berechnung von Auftriebsbeiwert ,
Widerstandsbeiwert und Momentenbeiwert eines gegebenen Profils bei
subsonischen Strömungsgeschwindigkeiten (Ma<1). XFOIL verwendet standardmäßig den
Wert 8aj = 9, was laut XFOIL USER GUIDE dem Turbulenzgrad eines durchschnittlichen
Windkanals entspricht. Größere Werte von 8aj bedeuten weniger Turbulenz, kleinere
dagegen einen höheren Turbulenzgrad. Bei der Simulation mit XFOIL kann über den als
critical amplification ratio bezeichneten Parameter 8aj Einfluss auf den Turbulenzgrad
genommen werden.
5.1 XFOIL-Ergebnisse
Um die Beiwerte zu ermitteln, muss die Anströmgeschwindigkeit F in der Rotorebene aus
dem Kapitel 2.7 in jedem Profilschnitt bekannt sein. Ebenso muss die jeweilige Reynolds-Zahl
mit der Länge des Profils bei gegebenen Stoffdaten des Mediums ermittelt werden.
Der aerodynamische Anstellwinkel ist die Differenz aus dem Anströmwinkel in der
Rotorebene und dem Bauwinkel des Profils im Modell.
Y Y(9) − YZ:(9). (2.18)
Die Beiwerte werden benötigt, um die vorherrschenden Kräfte aus der Tragflügeltheorie am
Rotorblatt herzuleiten, siehe Abbildung 5.1.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 5.1 – Druckverlaufplot des NACA 4415 Profils aus XFOIL
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
6. Ergebnisse
Bei den 3D-Simulationen werden zunächst die Residuen untersucht, um zu schauen, ob die
Lösung ausreichend konvergiert ist. Wenn dies der Fall ist, wird die Simulation auf eine
mögliche Verletzung der Kontinuitätsgleichung überprüft und der Einfluss des
Profilkoeffizienten auf das Rotordrehmoment und den axialen Geschwindigkeiten im
Fernbereich untersucht. Nach diesen Parameterstudien werden die Strömungsverhältnisse
in der Meridianansicht mittels Stromlinienverläufen entlang der Blattoberfläche visualisiert.
Ebenso wird der Drall, der sich entgegen der Umfangsgeschwindigkeit ausbreitet, quantitativ
erfasst, um eine Leistungsabschätzung nach der Eulerschen Strömungshauptmaschinen-
gleichung durchführen zu können. Anschließend werden die analytischen Auftriebs- und
Widerstandskräfte einzelner Teilflächen des Rotors in Schub- und Umfangskraft überführt
und den vorherrschenden Kräften im Modell gegenübergestellt.
Abschließend werden die theoretischen Leistungsberechnungen von Betz, Schmitz und Euler
mit den ermittelten mechanischen Leistungen aus der Simulation verglichen.
Die Modelle sind wie schon erwähnt nicht transient simuliert worden und deshalb sind
Strömungsgrößen nur ortsabhängig und nicht zeitabhängig. Man könnte es mit einer
Momentaufnahme der Fluidteilchen in einem Strömungsfeld erklären, indem die Strömung
durch Stromlinien, eine Raumkurve, welche die Geschwindigkeitsvektoren berührt, darstellt.
Die transiente Berechnung, gekennzeichnet durch die zeitliche Ableitung, brach schon nach
wenigen Iterationen ab. Der Solver gab den Hinweis, keine Überlappung (permit no
intersection) zwischen den Schnittstellen (Interfaces) der rotierenden Domain und der
Stationären zu zulassen. Leider führte auch diese Option zu keinem Ergebnis und die
Simulation brach erneut ab.
Hier nochmal eine Übersicht der Modelle:
Modellname c1 [m/s] n [1/min] λ Netzelemente
Modell I 5/10/12 12,5/25/30 7,92 ~2,3 106
Modell II 5/10/12 14/28/33,6 7,07 ~4,1 106
Tabelle 6.1 – Übersicht der Modellparameter
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Die verwendeten Stoffdaten von Luft, die für die Simulation verwendet wurden, sind im
Anhang aufgeführt. Für die Auswertung wurden die Geschwindigkeiten bezogen auf das
kartesische Koordinatensystem mit dem Turbomaschinentool in Zylinderkoordinaten
(o, m, & − , 9, &) übertragen.
Um die Fülle an Daten der numerischen Berechnung besser auswerten zu können, werden
Stützstellen verwendet, siehe Abbildung 6.1.
Abbildung 6.1 – Einführung von Stützstellen zur Auswertung der Ergebnisse
Die Orte A und B liegen dabei in den stationären Teilen des Raums, wobei A mit 100 Metern
weit vor dem Rotor liegt und B 150 Meter dahinter. Die Orte 1 und 2 liegen in der
rotierenden Domain nahe der Rotorebene, 3 Meter davor und 2 Meter dahinter.
6.1 Residuenstudie
Die Residuen bilden die Veränderung der Ergebnisse eines Iterationsschritts zum Vorherigen
ab. Dabei stammen die Residuen U-Mom, V-Mom und W-Mom aus der Impulsgleichung in x,
y und z-Richtung, während P-Mass der Kontinuitätsgleichung entspricht.
Graphisch werden diese dann über die Anzahl der Iterationsschritte aufgetragen. Die
iterative Rechenoperation kann durch Erreichen eines bestimmten Residuum
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
(Konvergenzkriterium) oder durch die vorgegebene Anzahl an Iterationsschritten beendet
werden. Bei dem Modell I wurde ein Konvergenzkriterium von 1e-6 eingestellt und maximal
tausend Iterationsschritte durchgeführt. Modell II wurden nur noch 800 Iterationen
durchgeführt, da die Residuen einen relativ konstanten Verlauf hatten.
Abbildung 6.2 - Residuen von Modell I - =10m/s, λ=7,92
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.3 - Residuen vom Modell II - =10m/s, λ=7,02
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
6.2 Stromröhrenaufweitung nach Betz / Froude-Rankine
Um die Verletzung der Massenerhaltung zu überprüfen, wurde der Volumenstrom am Inlet
(Kreisfläche am Eintritt), mit denen am Entrainment (Mantelfläche) und am Opening
(Kreisfläche am Austritt) verglichen, siehe Abbildung 6.4. Nach dem Kontinuitätsprinzip muss
ijk − ( :-,ika + :-,#k) 0 (6.1)
ergeben.
Der prozentuale Fehler lag dabei unter 0,01 % und ist vernachlässigbar.
Abbildung 6.4 – Schematischer Aufbau zur Prüfung der Massenerhaltung
Eine andere Möglichkeit zur Überprüfung der Massenerhaltung kann über die Stromröhre
geschehen. Dabei werden an den zwei Orten (A, B) die arithmetischen Mittel der
Geschwindigkeiten ZljZ über den Querschnitt der Stromröhre gebildet
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
ZljZ 1`ZljZ,j
j.
(6.2)
Mit der mittleren Geschwindigkeit und der Querschnittsfläche an den Orten A und B lassen
sich wieder die Massenströme bestimmen, wobei
ZljZ, − ZljZ, 0, (6.3)
sein muss.
Die Prüfung der Massenerhaltung durch die Stromröhre ergab einen prozentualen Fehler
von etwa 1% und ist ebenfalls vernachlässigbar.
Wie in Kapitel 2 erwähnt, wird die Windgeschwindigkeit c1 auf c3 abgebremst. Umgekehrt
proportional zur Strömungsgeschwindigkeit nimmt die Querschnittsfläche der Stromröhre
von AA nach AB zu (Abbildung 6.5).
Abbildung 6.5 – Aufweitung der Stromröhre, die schwarze Linie stellt die Position des Rotors dar,
λ=7,92
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Da die Konstruktion des Rotors nicht auf der Betzschen Optimalauslegung beruht, wird die
erwartete axiale Geschwindigkeitsverzögerung weit hinter dem Rotor von zwei Drittel nicht
erreicht.
Tabelle 6.2 – Parameter der Stromröhrenaufweitung, Modell I
Abbildung 6.6 - Stromlinienplot der Absolutgeschwindigkeit– Aufweitung der Stromröhre, Modell
=10m/s, λ=7,92
In Abbildung 6.7 ist die Strömungsverzögerung in Richtung der Rotorachse quantitativ
erfasst worden. Weit vor dem Rotor ist die Geschwindigkeit noch auf dem Niveau der
c1 [m/s] [m²] c3 [m/s] [m²] [/, /, 5 2290 3,1 3696 22902 22513
10 2290 6,1 3696 11451 11288
12 2290 7,3 3696 27483 27018
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
eingestellten Strömungsgeschwindigkeit, nimmt aber infolge der Stauung zum Rotor hin ab.
Der Geschwindigkeitsanstieg unmittelbar hinter dem Rotor entsteht, wenn die Strömung die
Stauung durch den Rotor überwunden hat und drallbehaftet fortgetragen wird.
Abbildung 6.7 – Änderung der axialen Absolutgeschwindigkeit in z-Koordinatenrichtung
Vergleicht man die Aufweitung der Stromröhre mit der axialen Strömungsgeschwindigkeit,
wird die Proportionalität durch die Massenerhaltung deutlich. An dieser Stelle sei erwähnt,
dass die axialen Geschwindigkeiten aus Abbildung 6.7 und 6.8 immer auf die gleiche
Querschnittfläche bezogen worden sind.
Das Modell II mit der Schnelllaufzahl 7,07 weist, entgegen der Erwartungen, einen ganz
anderen axialen Geschwindigkeitsverlauf auf (Abbildung 6.8). Anscheinend hat bei diesem
Modell die Übergabe der Geschwindigkeitskomponenten von der rotierenden Domain in die
Stationäre nicht funktioniert, wobei das gleiche Interface (GGI) wie in Modell I verwendet
wurde. Der rotierende Bereich ist davon jedoch nicht betroffen.
2
3
4
5
6
7
8
9
10
11
12
0 25 50 75 100 125 150 175 200 225 250 275
c_axial [m/s]
Z [m]
Axiale Absolutgeschwindigkeit - λ=7,92
c1=10 m/s c1=5 m/s c1=12 m/s Rotor
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.8 – Änderung der axialen Absolutgeschwindigkeit in z-Koordinatenrichtung
Abbildung 6.9 bildet die axiale Geschwindigkeitsverteilung über den Radius ab. Die axiale
Geschwindigkeit wurde dabei auf den Ringflächen, die aus zwei konzentrischen Kreisen mit
dem Abstand dr = 1 m besteht, gemittelt. Weit vor dem Rotor (A) ist die Geschwindigkeit
konstant. Am Rotoreintritt (1) wird die Strömung durch die Nabe gestaut und erfährt eine
Verzögerung. Im Bereich des Rotors von Profil NACA4450 bis NACA4415 ist die Schwankung
jedoch wieder relativ gering. Am Rotoraustritt (2) entsteht durch die Kante an der
Nabenrückseite eine Ablösung mit einer Totzone, in der eine Rückströmung stattfindet.
Durch die konvexe Geometrie der Nabe am Rotoreintritt erfährt die Strömung eine
Beschleunigung, die ihr Maximum kurz vor dem NACA4450 erreicht. Danach stellt sich
erneut ein relativ konstanter Wert über den gesamten Rotor ein. Weit hinter dem Rotor (B)
nimmt die Geschwindigkeit mit zunehmendem Radius zu.
2
3
4
5
6
7
8
9
10
11
12
0 25 50 75 100 125 150 175 200 225 250 275
c_axial [m/s]
Z [m]
Axiale Absolutgeschwindigkeit - λ=7,07
c1=10 m/s c1=5 m/s c1=12 m/s Rotor
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.9 – Änderung der axialen Absolutgeschwindigkeit über den Radius an den vier
Stützstellen
-6
-4
-2
0
2
4
6
8
10
12
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
c_axial [m/s]
r [m]
Axiale Absolutgeschwindigkeit, c=10m/s, λ=7,92
A B 1 2 Nabe
NACA4450 NACA4425 NACA4420 NACA4415
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
6.3 Vergleich der Drehmomente
Das Modell I wurde mit drei verschiedenen Windgeschwindigkeiten und Drehzahlen und
unterschiedlichen Blattneigungswinkeln simuliert. Dabei ist das höchste Drehmoment bei
einem Blattneigungswinkel von 0°, ausgehend von der Lage des Profils NACA 4415 im
Blattspitzenbereich, ermittelt worden. Die übrigen simulierten Neigungswinkel wirken sich
dagegen negativ auf die Aerodynamik des Rotors bei gegebener Nenndrehzahl aus, siehe
Abbildung 6.10. Der negative Anstellwinkel bedeutet, dass die Profilnase in Windrichtung
geneigt ist, positive Anstellwinkel entsprechend umgekehrt.
Abbildung 6.10 - Drehmoment in Abhängigkeit vom Anstellwinkel
Bei dem Modell II, mit fast doppelt so vielen Elementen, ist das ermittelte Drehmoment
wesentlich höher und zeigt gute Übereinstimmungen mit den analytischen Werten aus der
Blattelemententheorie, siehe Tabelle 6.11.
0
50000
100000
150000
200000
250000
300000
350000
-10 -8 -6 -4 -2 0 2 4 6 8 10
Dre
hm
om
en
t [
Nm
]
Anstellwinkel [°]
Drehmoment bei verschiedenen Anstellwinkeln
Modell I - 5 m/s
Modell I - 10 m/s
Modell I - 12 m/s
Modell II - 5 m/s
Modell II - 10 m/s
Modell II - 12 m/s
Poly. (Modell I - 5 m/s)
Poly. (Modell I - 10 m/s)
Poly. (Modell I - 12 m/s)
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Modell [m/s] HH [Nm] N [Nm]
Modell I ,
5 49184 37717
10 207738 151875
12 301666 225588
Modell II ,
5 42178 52080
10 210945 203143
12 283914 302400
Tabelle 6.3 – Vergleich - Analytisches und numerisches Drehmoment
Anscheinend hat die schlechte Diskretisierung in der rotierenden Domain zur Verfälschung
der numerischen Werte bei Modell I geführt (Abb. 6.12).
Abbildung 6.11 – Beispiel eines Vergleichs zwischen guter und schlechter Diskretisierung der
Knotenpunkte
0 50 100 150 200 250 300-5
0
5
10
15
20
25
x[m]
c[m
/s]
schlechte Diskretissierung
verbesserte Diskretisierung
Spline "schlechte Diskr."
Spline "verbesserte Diskr."
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Beim Vergleich der Plots (Abbildung 6.13 und Abbildung 6.15) im selben Rotorschnitt wird
deutlich, dass die grobe Netzstruktur massiven Einfluss auf die Rundung der Profilnase hat,
denn bei Modell I ist diese deutlich eckiger. Durch die Unstetigkeit der Krümmung können
Druck- und Geschwindigkeitssprünge entlang des Profils entstehen. Ebenso berücksichtigt
die Tragflügeltheorie keine radialen Strömungsvorgänge, die natürlich real sowie in der
Simulation vorhanden sind.
Die hohen Unterschiede der Drehmomente zwischen den beiden Modellen können aber
nicht mit der Schnelllaufzahl zusammenhängen, denn das erste Modell wurde mit
unterschiedlichen Anstellwinkeln simuliert.
Abbildung 6.12 – Stromlinienplot der Relativgeschwindigkeit, Modell II - NACA 4420 mit r = 17 m,
HN = 2°
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.13 – Vektorplot der Relativgeschwindigkeit, Modell II - NACA 4420 mit r = 17 m,
HN = 2°
Die Abbildungen 6.14 – 6.21 zeigen die Relativgeschwindigkeiten in Form von Stromlinien
und Vektoren auf einer Ebene im Profil NACA 4420 bei unterschiedlichen Anstellwinkeln von
Modell I. Dabei wird ersichtlich, dass sich die maximale Relativgeschwindigkeit, mit
abnehmenden Anstellwinkeln, von der Profilnase saugseitig in Richtung der Hinterkante
verschiebt. Die Auftriebskraft ist dadurch größtenteils nicht mehr, wie gewünscht in
Umfangsrichtung vorhanden, sondern in weist in Richtung der Rotationsachse. Zudem
verringert sich die Strömungsgeschwindigkeit auf der Saugseite, durch den abnehmenden
Anstellwinkel.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.14 – Stromlinienplot der Relativgeschwindigkeit, Modell I - NACA 4420 mit r = 17 m,
HN = 2°
Abbildung 6.15 – Vektorplot der Relativgeschwindigkeit, Modell I - NACA 4420 mit r = 17 m,
HN = 2°
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.16 – Stromlinienplot der Relativgeschwindigkeit, Modell I - NACA 4420 mit r = 17 m,
HN = 4°
Abbildung 6.17 – Vektorplot der Relativgeschwindigkeit, Modell I - NACA 4420 mit r = 17 m,
HN = 4°
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.18 – Stromlinienplot der Relativgeschwindigkeit, Modell I - NACA 4420 mit r = 17 m,
HN = 6°
Abbildung 6.19 – Vektorplot der Relativgeschwindigkeit, Modell I - NACA 4420 mit r = 17 m,
HN = 6°
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.20 – Stromlinienplot der Relativgeschwindigkeit, Modell I - NACA 4420 mit r = 17 m,
HN = 10°
Abbildung 6.21 – Vektorplot der Relativgeschwindigkeit, Modell I - NACA 4420 mit r = 17 m, HN
= 10°
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
6.5 Drallbestimmung nach EULER
Wie in Kapitel 2.4 erwähnt, lässt sich die Leistung einer axialen Turbine auch durch die
Eulersche Strömungsmaschinenhauptgleichung bestimmen. Die Anwendung der Gleichung
setzt jedoch eine schaufelkongruente Strömung voraus. Die absolute Geschwindigkeits-
komponente in Umfangsrichtung, die für den Drall verantwortlich ist, wurde über
Ringflächen (dr=1 Meter) gemittelt. Weit vor dem Rotor (A) ist die Strömung drallfrei. Am
Rotoreintritt (1) schwankt die Geschwindigkeit um den Nullpunkt, ist aber vernachlässigbar
klein und wird ebenfalls als drallfrei angenommen. Der Drall am Rotoraustritt (2) ist nun
entscheidend für die Leistungsbestimmung nach EULER.
Abbildung 6.22 – Änderung der Absolutgeschwindigkeit in Umfangsrichtung
Schaut man sich den Kurvenverlauf des Dralls direkt hinter dem Rotor an, lässt dieser sich ab
einem Radius von 6 Metern am ehesten mit einem Potentialwirbel vergleichen, bei dem die
-1
0
1
2
3
4
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
c2u [m/s]
r [m]
Absolutgeschwindigkeit in Umfangsrichtung , c=10m/s, λ=7,92
A B 1 2 Nabe
NACA4450 NACA4425 NACA4420 NACA4415
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Strömungsgeschwindigkeit in Umfangsrichtung mit zunehmendem Radius abnimmt, siehe
Abbildung 6.23/6.24.
Abbildung 6.23 – Potentialwirbel nach Hameln-Oseen [11]
Abbildung 6.24 – Stromlinienplot der Absolutgeschwindigkeit zur Darstellung des Dralls
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Um eine mögliche Verletzung der Betz-Theorie zu überprüfen, wurde die spezifische Arbeit
Y, unter Verwendung der Annahme, dass die axiale Geschwindigkeitskomponente sich auf
ein Drittel der Ausgangsgeschwindigkeit hinter dem Rotor verzögert, bestimmt. Für die
vorliegende Anlage ergeben sich gemäß dieses Vorgehens folgende Werte:
n 28 1/min
n_sec 0,467 1/s
rho 1,185 kg/m^3
D 54 m
c_Wind 10 m/s
P_Wind 1,357 MW
P_Betz 0,801 MW
m_pkt_ein 27139 Kg/s
m_pkt_aus 9046 Kg/s
m_pkt_mittel 18093 kg/s
Y_Betz,
P=m_pkt_mittel*Y 44,25
m²/s²
Tabelle 6.4 – Anlagenparameter zur Ermittlung von Y_Betz
Der Massenstrom kann aus dem arithmetischen Mittel zwischen Rotoreintritt und –austritt
gebildet werden.
ijk 9. (6.4)
:- 9 13 . (6.5)
18` j
k
j
(6.6)
Die spezifische Arbeit kann nun durch Umformung der Euler Gleichung berechnet werden.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
@ 5 .
(6.7)
Um den Drall über den Radius bestimmen zu können, muss die spezifische Arbeit über die
Fläche gemittelt werden.
@ 1@(9) +. (6.8)
@ 1 @(9)
¡¡+9+.
(6.9)
@ 2 @(9)
¡¡+9.
(6.10)
Dividiert man die Summe der mittlere spezifische Arbeit Y¢£¤u durch die gesamte
Rotorfläche erhält man wieder Y¢£¤u .
Unter der Annahme, dass die mittlere spezifische Arbeit über den Radius konstant bleibt,
kann man den Drall c2u bestimmen.
:(9) @; . (6.11)
In der Nähe der Nabe ergeben sich unrealistisch hohe Werte für die Drallkomponente c2u,
die größer sind als die Umfangsgeschwindigkeit des Rotors, so daß der Verlauf erst ab einem
Radius von etwa 3 Metern interpretierbar ist.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.25 – Analytische Drallbestimmung mittels Betz Theorie, Annahme spezifische Arbeit
ist konstant
Vergleicht man den Kurvenverlauf des analytisch ermittelten Dralls (Abb. 6.25) mit dem Drall
aus der numerischen Berechnung (Abbildung 6.22, an Stelle 2), erkennt man, dass diese sich
ab einem Radius > 6 Meter ähneln. Jedoch ist die ermittelte mittlere spezifische Arbeit aus
der numerischen Berechnung höher. Dies könnte mit einer abgelösten Profilströmung in der
Simulation zusammenhängen, die eine höhere Umlenkung ergibt, als die Theorie von Betz
zulässt. Es kann somit davon ausgegangen werden, dass die Strömung in der Simulation
nicht schaufelkongruent und die Eulersche Strömungsmaschinenhauptgleichung nicht
anwendbar ist.
Um Ablösungen generell abbilden zu können, muss das Profil von einer Prismenschicht
umgeben sein und eine hohe Diskretisierung vorhanden sein. Um die Fluktuation der
Ablösung zu visualisieren, bedarf es einer transienten Berechnung der Modelle.
Betriebspunkt ¥¦§¨ ¥N¦©ª«¬
c1 = 10 m/s, λ = 7,92 44,3 71,3
Tabelle 6.5 – Vergleich der spezifischen Arbeit nach Euler
y = 15,091x-1
0
2
4
6
8
10
12
14
16
18
20
0 5 10 15 20 25 30
c2u [m/s]
r [m]
c_2u
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
6.6 Kräfte am Rotor
Im Kapitel 2.5 werden die Auftriebs- und Widerstandskräfte am Rotor in Schub- und
Umfangskraft überführt. Um einen Vergleich mit den numerischen Werte herstellen zu
können, muss das Rotorblatt im Modell erst unterteilt und an den Teilflächen die Schub- und
Umfangskraft ermittelt werden (Abbildung 6.26).
Abbildung 6.26 – Teilflächen zur Ermittlung der Schub- und Umfangskräfte der Modelle
Für die analytische Bestimmung der Kräfte am Rotorblatt wurden die Beiwerte aus der 2D -
Simulation mit XFOIL (Kapitel 5) in die Blattelemententheorie miteinbezogen und die Länge
sowie die Beiwerte des Profils linear über den Radius interpoliert.
Ebenso wurden Verluste miteinbezogen, um möglichst realistische Werte zu erhalten. Die
Verluste, die miteinfließen, sind die Profilverluste
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
xa®¯j 1 − 329Q°± ,
(6.12)
die durch den Profilwiderstand resultieren. Ebenfalls wird der Tip-Verlust für eine
Auslegungsschnelllaufzahl ° 2 mit
xj# 1 − 1,84b° ,
(6.13)
bestimmt. Durch die Umströmung der Blattspitze von der Druckseite zur Saugseite, nimmt
der Auftrieb an dieser Stelle ab (Abb. 6.27) und es bildet sich ein aufweitender Wirbel, der
mit der Strömung davon getragen wird.
Abbildung 6.27 – Schematische Darstellung des Tip –Verlusts [Gasch]
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.28 –Vektorplot der Relativgeschwindigkeit an der Blattspitze
Abbildung 6.29 – Tip-Verlust an der Blattspitze, infolge des Druckunterschieds
Abbildungen 6.30 und 6.31 stellen die numerischen und die analytisch aus der
Tragflügeltheorie ermittelte Schubkräfte gegenüber. Das Modell II stimmt nahezu mit den
ermittelten Werten aus der Tragflügeltheorie überein. Bei Modell I hingegen könnte wieder
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
die schlechte Diskretisierung ein Grund für die hohe Abweichung sein. Im Naben- und
Blattspitzenbereich sind aber generell bei beiden Modellen Abweichungen zu verzeichnen.
Abbildung 6.30 – Vergleich der Schubkräfte ab einem Radius > 5 Meter
Abbildung 6.31 – Vergleich der Schubkräfte ab einem Radius > 5 Meter
0
500
1000
1500
2000
2500
3000
3500
0,00 0,10 0,20 0,30 0,40 0,50 0,60 0,70 0,80 0,90 1,00
N
r/R
Schubkraft (λ=7,92/c_wind=10m/s)
Fs analytisch Fs numerisch Fs analytisch - verlustbehaftet
0
500
1000
1500
2000
2500
3000
0,00 0,10 0,20 0,30 0,40 0,50 0,60 0,70 0,80 0,90 1,00
N
r/R
Schubkraft (λ=7,07/c_wind=10m/s)
Fs analytisch Fs numerisch Fs analytisch - verlustbehaftet
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Die Gesamtschubkräfte am Rotor in tabellarischer Form:
Modell [m/s] GM,HH [N] GM,N [N]
Modell I ,
5 31158 30025
10 130174 120233
12 190006 170644
Modell II ,
5 28536 27207
10 118849 113814
12 172151 164279
Tabelle 6.6 – Vergleich der Gesamtschubkräfte am Rotor
6.6 Anlagenkennlinien der Modelle
Abbildung 6. stellt die Leistung der numerisch berechneten Modelle mit der Leistung gemäß
der Theorie von Betz gegenüber. Obwohl bei Modell I bessere Gleitzahlen an den jeweiligen
Profilschnitten mit XFOIL ermittelt wurden, weist Modell II eine bessere Leistungsausbeute
auf. Dies bestärkt wieder die Annahme, dass die Netzstruktur erheblichen Einfluss auf die
Ergebnisse hat.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.32 – Vergleich der Leistung nach der Theorie von Betz mit den Modellen
Aus den ermittelten numerischen Leistungen, können nun die Leistungsbeiwerte der
Modelle bestimmt werden.
c1 Leistung
[MW]
Leistung Wind
[MW]
Leistungsbeiwert
cp_Modell
Modell I
5 [m/s] 0,055 0,170 0,326
10 [m/s] 0,445 1,356 0,328
12 [m/s] 0,779 2,343 0,333
Modell II
5 [m/s] 0,068 0,170 0,400
10 [m/s] 0,532 1,356 0,392
12 [m/s] 0,950 2,343 0,406
Tabelle 6.7 – Vergleich von Leistung und Leistungsbeiwerte
0100200300400500600700800900
100011001200130014001500
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
Leistung [kW]
c_Wind [m/s]
Anlagenkennlinien der Modelle
P_Betz P_Mech Modell I P_Mech Modell II
Poly. (P_Mech Modell I) Poly. (P_Mech Modell II)
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
6.7 Das Q-Kriterium
Das Q-Kriterium (Wirbelidentifikation) beschreibt die Stärke des Wirbels, indem die
Scherrate S und die WirbelstärkeΩ in folgender Beziehung stehen:
µ 12(Ω − ) 0
(6.14)
Überwiegt die Wirbelstärke der Scherung, dann ist Q > 0. Das Q-Kriterium wird verwendet,
um mögliche Ablösungen darstellen zu können.
Abbildung 6.33 – Q-Kriterium von 0,0014 druckseitig
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Abbildung 6.34 – Q-Kriterium von 0,0014 saugseitig
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
7. Zusammenfassung
In der vorliegenden Arbeit wurde sowohl eine Windenergieanlage vom Durchmesser 54 m
aerodynamisch dimensioniert, konstruiert als auch unter aerodynamischen Gesichtspunkten
simuliert. Festigkeitsbetrachtungen erfolgten keine.
Qualitativ wurden Strömungsgrößen ermittelt, die mit der Theorie von Betz
(Anströmgeschwindigkeit wird abgebremst) und Rankine-Froude (Stromröhre weitet sich
auf, da in der Rotorebene eine mittlere Geschwindigkeit aus An- und Abströmung gilt)
übereinstimmen. Quantitativ können die Berechnungsergebnisse allerdings nur grob
bewertet werden, da es sich bei der Windturbine um eine fiktive Maschine handelt, zu der
keine Messdaten vorliegen.
Die Maschine wurde zwar als rotationssymmetrischer Drittelschnitt ökonomisch mit
periodischen Randbedingungen in Umfangsrichtung simuliert, dennoch war die
Netzauflösung mit Netzschrittweiten von bis zu 2,5 m erheblich zu grob, um Detailanalysen
zu ermöglichen. Weder Strömungsablösung an einzelnen Profilschnitten noch
aerodynamische Verluste können aufgrund des erheblich zu groben Netzes zuverlässig
bewertet werden. Dennoch zeigen die Ergebnisse ein prinzipiell richtiges Vorgehen, das
ohne weiteres mit einem erheblich höheren Rechenaufwand zu realistischen Ergebnissen
führen würde. Vermutlich muss das Gitter um einen Faktor 10 bis 100 verfeinert werden, so
dass mit Rechenzeiten von einigen Tagen eventuell sogar mehreren Wochen, statt wie im
Rahmen der vorliegenden Studie von einigen Stunden zu rechnen ist.
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Literaturverzeichnis
[1] Bundesverband der Energie- und Wasserwirtschaft e.V. (2010). 17% Anteil der
Erneubaren Energien an der Deckung des Strombedarfs in 2010. Berlin.
[2] Betz, A. (1926). Windenergie und ihre Ausnutzung durch Windmühlen. Göttingen.
[3] Gasch, R. / Twele, J. (2010). Windkraftanlagen. Wiesbaden: Teubner.
[4] Drela, M. (2000). XFOIL. Abgerufen am 17. Januar 2013 von
http://web.mit.edu/drela/Public/web/xfoil/
[5] Jens Trapp / Robert Zores. (kein Datum). NACA 4 Digits Series. Abgerufen am 9.
September 2012 von http://www.ppart.de/aerodynamics/profiles/NACA4.html
[6] Siekmann, H. E. (2000). Strömungslehre. Heidelberg: Springer.
[7] Hau, E. (2008). Windkraftanlagen - Grundlagen, Technik, Einsatz, Wirtschaftlichkeit.
Heidelberg: Springer.
[8] Kameier, F. (1999/2010). Vorlesung Strömungslehre.
[9] Sigloch, H. (2008). Technische Fluidmechanik. Berlin, Heidelberg, New York: Springer.
[10] ANSYS. Handbuch / Hilfe / Support
[11] Schade, H. / Kunz, E. (2007). Strömungslehre. Berlin: de Gruyter
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Anhang
Stoffdaten der Luft
Temperatur 25°
Dichte 1,185 kg m-³
Kin. Viskosität 1,5452 10-5 m2 s--1
Umgebungsdruck 1 Bar
Ermittlung der Auftriebs- und Widerstandsbeiwerte mit XFOIL
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
Plots
0
500
1000
1500
2000
2500
3000
3500
4000
0,00 0,10 0,20 0,30 0,40 0,50 0,60 0,70 0,80 0,90 1,00
N
r/R
Schubkraft (λ=7,07/c1=12m/s)
Fs analytisch Fs numerisch Fs analytisch - verlustbehaftet
0
100
200
300
400
500
600
700
800
0,00 0,10 0,20 0,30 0,40 0,50 0,60 0,70 0,80 0,90 1,00
N
r/R
Schubkraft (λ=7,07/c_wind=5m/s)
Fs analytisch Fs numerisch Fs analytisch - verlustbehaftet
Windenergieanlagen – aerodynamische Auslegung und numerische Berechnung mittels Ansys CFX
0
100
200
300
400
500
600
700
800
0,00 0,10 0,20 0,30 0,40 0,50 0,60 0,70 0,80 0,90 1,00
N
r/R
Schubkraft (λ=7,07/c_wind=5m/s)
Fs analytisch Fs numerisch Fs analytisch - verlustbehaftet
0
500
1000
1500
2000
2500
3000
3500
4000
0,00 0,10 0,20 0,30 0,40 0,50 0,60 0,70 0,80 0,90 1,00
N
r/R
Schubkraft (λ=7,07/c1=12m/s)
Fs analytisch Fs numerisch Fs analytisch - verlustbehaftet