Æw. 4-2012 - Przep³yw w stopniu turbiny

Ćw. 4 Modelowanie przepływu w stopniu

turbiny

W  rzeczywistości  stopień  turbiny  jest  układem
trój-wymiarowym.  Składa  się  on  z  nieruchomych
łopatek  kierowniczych  (stator)  i  ruchomego
wirnika (rotora).

Podobnie  jak  w  przypadku  spręŜarki  osiowej
dokonujemy  uproszczenia:  analizujemy  przepływ
dla  jednej  pary  łopatek  na  pewnym  promieniu.
Przepływ

trójwymiarowy

zastępujemy

więc

przepływem dwu-wymiarowym.
Obliczenia

przeprowadzimy

dla

segmentu

obejmującego jedną łopatkę nieruchomą i jedną
ruchomą na powierzchni cylindrycznej którą
łopatki przecinają. Nasz stopień turbiny składa się
zatem z 2 elementów.

Tworzenie geometrii – GAMBIT
Do wykonania siatki wykorzystamy gotową
geometrię, pokazaną na rysunku, którą
zaimportujemy do Gambita (geometria.dbs).

-2;3.86

-2;-3.86

6.46;-1.76

6.46;-9.47

6.46;-2.30

6.46;-10.01

6.98;-5.47

17;-4.81

17; 2.90

1. Wykorzystując istniejące krawędzie tworzymy 4
powierzchnie: stator, lop-stator, rotor, lop-rotor

lop-stator

lop-rotor

rotor

stator

2. Od powierzchni stator odejmujemy (Subtract
Real Faces) powierzchnię lop-stator a od
powierzchni rotor – lop-rotor

rotor

stator

3. Dzielimy krawędzie spływu obu łopatek, aby
moŜna było im nadać osobne nazwy w części
górnej i dolnej – opcja Split Edge with Point

4. Linkujemy odpowiadające sobie krawędzie
statora i rotora ze względu na późniejsze
zastosowanie warunku brzegowego typu
PERIODIC (jest to konieczne przed siatkowaniem
krawędzi) –Mesh > Edge > Link Edge Meshes z
włączoną opcją Periodic.

Uwaga: śeby operacja

linkowania była poprawna, obie krawędzie muszą
mieć ten sam kierunek! (jeśli tak nie jest, to
koerunek krawędzi moŜna zmienić naciskając
Shift i środkowy przycisk myszy).

rotor

stator

Link

5. Siatkujemy krawędzie jak pokazano na
rysunku. Uwaga: liczba węzłów na wszystkich
krawędziach powinna być parzysta, z uwagi na
rodzaj siatki (Quad/Pave)

100 succ ratio 1

110 succ ratio 1

80 succ ratio 1

6. Siatkujemy obie powierzchnie powierzchnie
(Quad/Pave) Liczba elementów (Face): stator – ok.
5470, rotor – ok. 5130

7. Nadajemy warunki brzegowe

wlot

pressure_inlet

wylot

pressure_outlet

rotor-d

wall

rotor-g

wall

stator-d

wall

stator-g

wall

inter-stator

interface

inter-rotor

interface

periodic-rotor

periodic

periodic-stator

periodic

8. Nadajemy warunki na Continuum (obie
powierzchnie jako typu Fluid ale z osobnymi
nazwami „stator” i „rotor”)

9. Zapamiętujemy geometrię, siatkę i warunki
brzegowe a następnie eksportujemy siatkę dwu-
wymiarową.

Obliczenia - Fluent

Uruchomić Fluenta w wersji: 2D, Double
Precision, Serial

Wczytanie, sprawdzenie i przeskalowanie siatki:
wymiary siatki w cm

A. Ustawienia ogólne

General:
Density Based, Absolute, Steady, Planar
Ustawić jednostki ciśnienia (Units): bar (10

Pa)

Włączyć równanie energii
Models > Energy > Edit: Energy Equation: On
Ustawić model lepkości płynu
Models > Viscous > Edit:
model turbulencji Spalart-Allmaras (1 eqn): On
Wybrać rodzaj płynu
Materials > Fluid > Create/Edit
Fliud Fluent Materials: air, Density: ideal-gas >
Change/Create > Close
Ustawić warunki odniesienia
Cell Zone Conditions > Operating Conditions:

Operating Pressure = 0 bar

Warunki brzegowe:
Boundary Conditions > Edit
wlot: pressure_inlet: Gauge Total Pressure = 5 bar,
Supersonic/Initial Gauge Pressure = 4.9 bar, Total
Temperature = 900 K
wylot: pressure_outlet: Gauge Pressure = 3.5 bar,
Backflow Total Temperature = 900 K

Ustalenie interfejsów:
Mesh Interfaces > Create/Edit:
utworzyć interface z krawędzi inter-stator oraz inter-
rotor z włączoną opcją (Interface Options) Periodic
Repeats

B. Obliczenia dla pozornie ruchomego
wirnika

Cel: znalezienie takiej prędkości ruchu rotora, przy
której prędkość na wylocie z łopatek rotora będzie
pozioma (taka prędkość powinna być na wlocie do
łopatek kierownicy następnego stopnia).

Cell Zone Conditions > Rotor > Edit
w zakładce Motion i w polu Motion Type ustawiamy
opcję Moving Reference Frame, poniŜej ustawiamy
wartość prędkości ruchu łopatki (Translational
Velocity Speed) Y = -170 m/s

Inicjalizacja
Inicjalizację przeprowadzamy z warunków wlotu.
Iteracje
Wykonujemy ok. 400 iteracji (do osiągnięcia
zbieŜności 10

-3

potrzeba ok. 1000 iteracji) i

sprawdzamy wektory prędkości na wlocie, wlocie do
rotora (interface-rotor) oraz wylocie (scale= 20,
skip=4).

•

PoniewaŜ prędkość na wylocie do następnego
stopnia nie jest pozioma, oznacza to, Ŝe musimy
zwiększyć prędkość obrotu rotora. Zmieniamy ją na
190 m/s.

•

iterujemy następne 200 kroków (juŜ bez ponownej
inicjalizacji!) i ponownie sprawdzamy wektory
prędkości na wylocie. W razie potrzeby ponownie
zmieniamy ją na inną aŜ do osiągnięcia dobrego
wyniku. Wynik dla prędkości Y=-250 m/s
pokazano poniŜej

Analizujemy wyniki obliczeń:

•

pole ciśnień

•

pole prędkości (liczby Macha)

•

wektory prędkości (scale = 10, skip = 20)

•

prędkości (liczby Macha) na wlocie i wylocie

•

ciśnienia na wlocie i wylocie

•

rozkład ciśnień na łopatce rotora (w kierunku X

oraz Y)

Do obliczenia wartości współczynników siły nośnej
i oporu potrzebne jest ustalenie wielkości
referencyjnych - moŜliwe są róŜne warianty:

•

wariant 1 (z wlotu do kanału statora)

•

wariant 2 (z wlotu na łopatki rotora)

W przypadku analizowania łopatki rotora przyjmiemy
jako referencyjne dane z wariantu 2, tzn.
Reference Values >
Area (m2) = 0.083 (iloczyn cięciwy łopatki i
wymiaru w głąb=1 m)
Density (kg/m3) = 1.68 (przyjmiemy średnią wartość
obliczoną dla interface-rotor)
Depth (m) = 1 (wymiar w głąb)
Enthalpy (J/kg) – nie musimy obliczać wartości,
poniewaŜ nie jest potrzebna do wyznaczenia
współczynników
Length (m) = 0.083 (równy cięciwie łopatki rotora)
Pressure (bar) = 4.15 (przyjmiemy średnią wartość
obliczoną dla interface-rotor)
Temperature (K) = 858 (przyjmiemy średnią wartość
obliczoną dla interface-rotor)
Velocity (m/s) = 310 (przyjmiemy średnią wartość
obliczoną dla interface-rotor)
Pozostałe wielkości parametrów (lepkość i współ. C

)

pozostawimy bez zmiany.

Rozkłady współczynnika ciśnienia na powierzchni
łopatek statora i rotora (Pamiętajmy, Ŝe wartości są
poprawne tylko dla rotora z uwagi na przyjęte
parametry referencyjne)

Wartości sił i współczynników: siły nośnej oraz
oporu dla rotora

Reports > Forces >
W kierunki poziomym

W kierunku pionowym

C. Przepływ nieustalony

Zmieniamy ustawienia wirnika
Cell Zone Conditions > Rotor > Edit
w zakładce Motion i w polu Motion Type zmieniamy
ustawienie z Moving Reference Frame na Moving
Mesh. Prędkości ruchu łopatki (Translational
Velocity Speed) Y = -250 m/s pozostawiamy bez
zmiany

Ustawienie sceny wizualizacji:
Zwiększamy liczbę widocznych na ekranie kanałów
do 3
Display > Views > Periodic Repeats > Define
Periodic Type: Translational
Translation: X=0, Y= - 0.0771, Z=0
Number of Repeats: 3

Zmiana liczby widocznych okien graficznych
Z paska narzędzi naciskamy ostatnią ikonę (Arrange
the graphics window layout) i wybieramy opcję z 4
oknami.
Następnie wyłączamy wyświetlanie rezydułów
Monitors > Residuals > Off
Ustawienie animacji:
Zrobimy 2 filmy: pole ciśnień i pole liczb Macha
Calculation Activities > Solution Animations >
Create/Edit
W polu Animation Sequences ustawiamy wartość =2.

Film 1 – Nazwę sequence-1 zmieniamy na cisnienie,
w polu Every pozostawiamy wartość 1, w polu When
zmieniamy Iteration na Time Step. Naciskamy
Define i ustawiamy dalej.
W oknie Window pozostawiamy 1 i naciskamy Set.
W polu Display Type wybieramy Contours,

naciskamy Edit i przechodzimy do nowego okna. W
polu Contours of wybieramy Pressure i Static
Pressure. Naciskamy Display. W oknie nr 1 ukazuje
się obraz pola ciśnienia. Jeśli chcemy zmienić
wielkość wyświetlanego obrazka naciskamy klawisz
Fit To Window. Wychodzimy z tego okna naciskając
Close oraz z okna Animation sequences naciskając
Ok. Przechodzimy do definiowania drugiego filmu.

Film 2 – Nazwę sequence-2 zmieniamy na Mach. W
polu Every pozostawiamy wartość 1, w polu When
zmieniamy Iteration na Time Step. Naciskamy
Define i ustawiamy dalej.

W oknie Window zmieniamy numer okna na 2 i
naciskamy Set W polu Display Type wybieramy
Contours, naciskamy Edit i przechodzimy do
nowego okna. W polu Contours of wybieramy
Velocity i Mach Number. Naciskamy Display. W
oknie nr 2 ukazuje się obraz pola liczb Macha. Jeśli
chcemy zmienić wielkość wyświetlanego obrazka
naciskamy klawisz Fit To Window. Wychodzimy z
tego okna naciskając Close oraz z okna Animation
sequences naciskając Ok.
Wychodzimy z okna Solution Animation przez Ok.

Monitorowanie przebiegu współczynnika siły
nośnej i oporu dla łopatek rotora
Siła nośna:
Monitors > Lift > Edit

W oknie Options pozostawiamy wciśniętą tylko opcję
Plot, w oknie Window ustawiamy numer okna na 3,
uaktywniamy opcję Write pozostawiając nazwę pliku
Cl-history. W oknie Wall Zones uaktywniamy
rotor-g i rotor-d. W oknie Force Vector ustawiamy
X=0 i Y=1 (interesuje nas siła zgodnie z kierunkiem
ruchu łopatek). Potwierdzamy ustawienia (Ok.).

Siła oporu:
Monitors > Drag > Edit

W oknie Options pozostawiamy wciśniętą tylko opcję
Plot, w oknie Window ustawiamy numer okna na 4,
uaktywniamy opcję Write pozostawiając nazwę pliku
Cd-history. W oknie Wall Zones uaktywniamy
rotor-g i rotor-d. W oknie Force Vector ustawiamy
X=1 i Y=0 (interesuje nas siła w kierunku
prostopadłym do kierunku ruchu łopatek).
Potwierdzamy ustawienia (Ok.).

Ustalenie kroku czasowego
Przed przystąpieniem do iteracji musimy ustalić
wielkość kroku czasowego. Wynika on z
następującego rozumowania:

Najmniejszy wymiar liniowy siatki w Gambicie
wynosi ∆x = 0.1, co we Fluencie (po przeskalowaniu)
daje
∆

x = 0.001 m = 1 mm

Krok czasowy wynika ze wzoru

∆

t = ∆x / (V + a)

Maksymalną prędkość przyjmiemy jako V = 420
m/s, a prędkość dźwięku (dla temperatury 900K)
przyjmiemy a = 600 m/s
∆

t = ∆x / (V + a) = 10-3/1020 = 0.98 x 10-6 s ~ 1e

-6

PoniewaŜ obliczenia przebiegałyby wtedy bardzo
wolno przyjmiemy krok czasowy 10 razy większy
tzn. ∆t = 10-5 s.

Dla prędkości przesuwu (unoszenia) rotora równej
250 m/s i odległości między łopatkami rotora równej
7.7 cm ( po przeskalowaniu) otrzymamy czas
odpowiadający przejściu tej odległości jako

T = 7.7 x 10-2 / 250 = 3.08 x 10-4 s
Uwzględniając przyjęty krok czasowy otrzymamy
ilość kroków odpowiadającą jednemu przejściu jako

TS = T / ∆t = 3.0 x 10-4 / 10-5 s = 30
Zatem w trakcie jednego cyklu otrzymamy 30
zdjęć podczas animacji

Iteracje:
Najpierw wyłączamy zapisywanie animacji i
wykonujemy ok. 150 kroków czasowych do ustalenia
się amplitud wartości współczynnika siły nośnej i
oporu (w skali makro)

Teraz włączamy zapisywanie animacji i wykonujemy
kolejne 50 kroków czasowych.

Wyniki (dla pół prędkości)

0 kroków

5 kroków

10 kroków

15 kroków

20 kroków

25 kroków

30 kroków