Numeryczna symulacja zjawisk przepływowych w turbinach cieplnych

ZESZYTY NAUKOWE

POLITECHNIKI ŚLĄSKIEJ

? .% 5 Ojoo

WŁODZIMIERZ WRÓBLEWSKI

NUMERYCZNA SYMULACJA ZJAWISK PRZEPŁYWOWYCH W TURBINACH CIEPLNYCH

ENERGETYKA

z. 132

GLIWICE

2000

POLITECHNIKA SLĄSKA ZESZYTY NAUKOM

Nr 1450

00

WŁODZIMIERZ WRÓBLEWSKI

NUMERYCZNA SYMULACJA ZJAWISK PRZEPŁYWOWYCH W TURBINACH CIEPLNYCH

Gliwice ²⁰⁰⁰

Prof. dr hab. inż. Ryszard. Qryboś Prof. dr hab. inż. Jerzy Krzyżanowski

/

f

k o l e g i u m r e d a k c y j n e

REDAKTOR NACZELNY — Prof. dr hab. Zygmunt Kleszczewski REDAKTOR DZIAŁU — Dr hab. inż. Andrzej Witkowski.

Profesor Politechniki Śląskiej SEKRETARZ REDAKCJI — Mgr Elżbieta Lesko

REDAKCJA Mgr Anna Błażkiewicz

REDAKCJA TECHNICZNA Alicja Nowacka

Wydano za zgodą Rektora Politechniki Śląskiej

PL ISSN 0372-9796

Wydawnictwo Politechniki Śląskiej ul. Akademicka 5

44-100 Gliwice tel./fax (0-32) 237-13-81

Dział Sprzedaży i Reklamy (0-32) 237-18-48

w w w . w v d a w n i c t w o . D o l s l . g l i w i c e . p l

wydawnictwo@polsl.gliwice.pl

N akl. 110 + 83 A rk . w yd. 14 A rk . d ru k . 13,375 P a p ie r offset. Id. III 7 0 x 100, 80 g O d d a n o do d r u k u 01.02.2000 r. P o d p is, do d ru k u 01.02.2000 r. D r u k ukończ, w lu ty m 2000 r.

F oto k o p ie, d r u k i o p ra w ę w ykonał „ R O L E K ” , G liw ice, ul. K azim ierz a W ielkiego 4

Spis treści

Strona

W a żn iejsze ozn aczen ia ... 9

R ozdział 1. C el i zakres pracy 11

R ozdział 2. M atem atyczny opis zasad zachow an ia 15

2.1. O gólne rów nania zachow ania 15

2.2. R ów nania zachow ania E ulera 17

2.3. R ów nania zachow ania dla przepływ u lepkiego 19

2.4. N orm ow anie w ielkości obliczeniow ych 20

2.5. K rzyw oliniow y układ w spółrzędnych 21

2.5.1. T ransform acja rów nań E ulera 23

2.5.2. T ransform acja rów nań N aviera-Stokesa 24

2.6. W arunki brzegow e 26

2.7. Z w iązki zam ykające dla przepływ u turbulentnego 32

2.8. R ów nanie stanu gazu 34

R ozdział 3. O b szar ob liczen iow y i sposoby jego d yskretyzacji ... 36

3.1. W prow adzenie 36

3.2. T ransform acja obszaru oraz podział m etod generacji

siatek 37

3.3. P odstaw ow e typy siatek regularnych 39

3.4. C echy regularnych siatek obliczeniow ych 41

3.5. R ów nania do generacji siatki regularnej 42

3.6. M etoda generacji siatki regularnej 44

3.6.1. B rzeg w ew nętrzny - rozkład punktów na profilu 45 3.6.2. R ozkład punktów na brzegu zew nętrznym ... 46 3.6.3. R ozw iązanie rów nań P oissona ... 46 3.6.4. G eneracja algebraicznej siatki prostej ... 48

3.7. G eneracja siatki trójw ym iarow ej 49

R ozdział 4. N u m eryczn e rozw iązania rów nań zachow an ia 51

4.1. B ilansow anie członów konw ekcyjnych 51

4.2. D ekom pozycja różnic strum ieni. M etody G odunow a 52 4.2.1. M etoda G odunow a z dokładnym rozw iązaniem

zagadnienia R iem anna 53

4.2.2. O kreślenie lokalnego zagadnienia R iem anna 54

4.2.3. P rzybliżenie akustyczne 59

4.2.4. M etoda O shera 60

4.2.5. M etoda R oe 61

4 Spis treści

4.3. S chem aty z rozdzieleniem w ektora strum ieni ... 65

4.4. P o rów nanie w ybranych schem atów num erycznych ... 67

4.5. P o łąc ze n ie schem atów FV S i FD S ... 74

4.6. S chem aty typu ‘pod p rą d ’ w yższego rzędu dokładności ... 75

4.7. Jaw na m etoda całkow ania w zględem czasu ... 82

4.7.1. O kreślenie kroku czasow ego ... 83

4.7.2. S chem at E ulera (Laxa) ... 84

4.7.3. S chem at M ac C orm acka ... 84

4.7.4. S chem at dw ukrokow y dła m etod ‘pod p rą d ’ w yższego rzęd u ... 85

4.7.5. S chem at R ungego-K utty ... 85

4.8. P rzejście do zagadnienia w ielow ym iarow ego ... 87

4.9. B ilansow anie czło n ó w lepkich ... 90

4.10. N um eryczna realizacja w arunków brzegow ych ... 92

4.10.1. W arunki brzegow e n a w locie ... 95

4.10.2. W arunki brzegow e na w ylocie ... 97

4.10.3. W arunki na ściance ... 100

4.10.4. W arunek okresow ości ... 101

4.10.5. P ołączenie obszarów w obliczeniach stopnia ... 103

4.10.6. Im plem entacja w arunków n ieo dbijających ... 108

4.11. D o b ó r w arunków początkow ych ... 113

R ozdział 5. O b lic z e n ia s tr u k tu r y p rz e p ły w u w k a n a ła c h tu r b in c iep ln y c h ... 118

5.1. O bliczenia przepływ u w turbinie pow ietrznej ... 119

5.2. O bliczenia ostatniego stopnia turbiny parow ej ... 129

5.3. O bliczenia przepływ u gazu lepkiego p rzez kanał łopatkow y ... 147

5.4. A naliza w arunków przepływ u na w locie do części niskoprężnej turbiny ... 153

R ozdział 6 . M o d e lo w a n ie p rz e p ły w u p a r y w o d n e j w k a n a ła c h ło p a tk o w y c h ... 163

6.1. M otyw acja ... 163

6.2. Z ależności dla czynnika rzeczyw istego ... 164

6.3. P orów nanie m odeli czynnika roboczego ... 167

6.4. M odelow anie przepływ u p ary w odnej z k o n d en sacją ... 172

6.4.1. T w orzenie zarodków kondensacji ... 173

6.4.2. O kreślenie udziału fazy ciekłej ... 174

6.4.3. M etoda obliczeń przepływ u pary w odnej z k o n d en sac ją ... 176

6.5. S truktura przepływ u p a iy w odnej z k o n d en sacją ... 178

6.5.1. P rzepływ p rzez płaski kanał łopatkow y ... 179

6.5.2. T rójw ym iarow a struktura przepływ u ... 191

Rozdział 7. K on cepcja m odułow ej analizy przepływ u w układach

łop atk ow ych turbin cieplnych 195

7.1. Wprowadzenie 195

7.2. Algorytm analizy przepływu w grupie stopni 196

Rozdział 8. P odsum ow an ie 200

L iteratu ra 203

S treszczen ie 213

S um m ary 214

Spis treści_________________________________ 5

Contents

Page

N o tatio n 9

Chapter 1. A im and scope o f the thesis 11

Chapter 2. M ath em atical d escrip tion o f the con servation law s ... 15

2.1. Basic conservation equations ... 15

2.2. Euler conservation equations ... 17

2.3. Conservation equations for the viscous flow 19 2.4. Non-dimensional quantities 2 0 2.5. Curvilinear coordinate system 21 2.5.1. Transformation o f the Euler equations ... 23

2.5.2. Transformation o f the Navier-Stokes equations ... 24

2.6. Boundary conditions 26 2.7. Closure relations for the turbulent flow 32 2 .8..Gas equation o f state 3 4 Chapter 3. C alculation domain and its discretisation methods ... 36

3.1. Introduction 3 6 3.2. Domain transformation and division o f generation methods ... 3 7 3.3. Basic types o f regular grids 3 9 3.4. Properties o f regular numerical grids 41 3.5. Equations for generation o f regular grid 42 3.6. Generation method o f regular grid ... 4 4 3.6.1. Inner boundary. Points distribution on the profile ... 4 5 3.6.2. Outer boundary 4 6 3.6.3. Poisson equations solution 46 3.6.4. Generation o f algebraic grid 4 8 3 .7 .3 D grid generation 4 9 Chapter 4. N um erical solution o f conservation equations 51 4.1. Calculation o f convective terms 51 4.2. Flux difference splitting. Godunov-type methods 52 4.2.1. Godunov method with exact Riemann solver 5 3 4.2.2. Description o f the local Riemann problem 5 4 4.2.3. Acoustic approximation 5 9 4.2.4. Osher method 60 4.2.5. Roe method 61 4.3. Flux Vector Splitting schemes 65 4.4. Comparison o f chosen numerical schemes 67 4.5. C om bination o f FVS and FD S schem es ... 74

4.6. U pw ind schem es o f high o rder accuracy ... 75

4.7. E xplicit m ethod o f integration in tim e ... 82

4.7.1. T im e step calculation ... 83

4.7.2. E uler (Lax) schem e ... 84

4.7.3. M ac C orm ack schem e ... 84

4.7.4. T w o-steps schem e for upw ind m ethods o f high order accuracy ... 85

4.7.5. R unge-K utta schem e ... 85

4.8. Step to 3D problem s ... 87

4.9. C alculation o f viscous term s ... 90

4.10. N um erical im plem entation o f boundary conditions ... 92

4.10.1. Inlet boundary conditions ... 95

4.10.2. O utlet boundary conditions ... 97

4.10.3. W all conditions 100 4.10.4. P eriodic conditions 101 4.10.5. D om ain connection for stage calculations 103 4.10.6. Im plem ent o f nonreflecting boundary conditions ... 108

4.11. S election o f initial conditions 113 C hapter 5. C alcu lation o f the flow stru cture in therm al turbine channels 118 5.1. C alculation o f the flow in the air turbine 119 5.2. C alculation o f last stage o f th e steam turbine ... 129

5.3. C alculation o f the viscous flow in the blade cascade 147 5.4. F low conditions analysis at the inflow to the LP part o f the turbine ... 153

C hapter 6 . M o d e llin g o f th e ste a m flow in b la d e ch a n n e ls 163 6 .1. M otivation ... 163

6.2. Real gas equations ... 164

6.3. C om parison o f gas m odels 167 6.4. M odelling o f the steam flow w ith condensation ... 172

6.4.1. N ucleation 173 6.4.2. W etness fraction description 174 6.4.3. C alculation m ethod o f steam flow w ith condensation ... 176

6.5. Steam flow structure w ith condensation 178 6.5.1. F low in the blade cascade ... 179

6.5.2. 3D structure o f the flow 191 C hapter 7. C oncep tion o f the m odular analysis o f the flow in th e rm a l tu r b in e b la d e c h a n n els 195 7.1. Introduction 195 7.2. A lgorithm o f analysis o f the flow in the group o f stages ... 196

8 Contents

Chapter 8. S u m m ary 200

R eferen ces 203

P olish su m m ary 213

E n glish su m m ary 214

Ważniejsze oznaczenia

a - prędkość dźw ięku

A,B, - m acierz

CFL - liczba C ouranta-Friedrichsa-L evy’ego

cp - pojem ność cieplna w łaściw a przy stałym ciśnieniu c, - pojem ność cieplna w łaściw a przy stałej objętości

e - p E

E - energia całkow ita, w ektor strum ienia (kierunek x) F - w ektor strum ienia (kierunek y)

G - w ektor strum ienia (kierunek z)

h - entalpia w łaściw a

H - uogólniony w ektor strum ieni, w ektor źródłow y J - jak o b ian transform acji

J ' - prędkość tw orzenia zarodzi kondensacji k - stała B oltzm anna

L - m acierz lewych w ektorów w łasnych L - ciepło utajone, lewy w ektor własny

I - droga m ieszania

Is - droga sw obodna

M - liczba M acha

n - w ektor norm alny

P

-

P r - liczba Prandtla

Q

-

-

R - stała gazow a

r - prom ień kropli, prom ień

r - prom ień krytyczny

Re - liczba R eynoldsa s - funkcja ograniczająca

S - przesycenie

t ^{- czas}

T - tem peratura

U, V, w - składow e prędkości w układzie bezw zględnym

u , v , w - składow e kontraw ariantne w układzie bezw zględnym U ' . V . W - składow e kontraw ariantne w układzie względnym W - w ektor zm iennych charakterystycznych

X.y.z - kierunki w k artezjańskiego układu w spółrzędnych

y - nierów now agow y stopień w ilgotności z - w spółczynnik ściśliw ości

a - w spółczynnik

P - w spółczynnik

A T - przechłodzenie

- kierunki transform ow anego układu w spółrzędnych

K - w ykładnik izentropy

X - w artość własna

A - m acierz w artości w łasnych

A - lepkość dynam iczna

P - gęstość

a - naprężenie pow ierzchniow e

z - naprężenie

Q * p rędkość kątowa, w irow ość

eksy

A - w ielkości w krzyw oliniow ym układzie w spółrzędnych + - odległość bezw ym iarow a od ścianki

0 - param etry całkow ite

c - w ielkość w układzie bezw zględnym

g - gaz

g r - graniczny

ij.k - indeksy siatki

l - lam inam y

max - m aksim um

min - m inim um

n - dyskretny poziom czasow y

s - nasycenie

t - turbulentny

x ,y,z,ę, ij,ę - p o chodna w kierunku x,y,z,Ł,,

//,

....( 0 -

w - w oda, w ielkość w układzie w zględnym , ścianka

w ew - w ew nętrzny

zew - zew nętrzny

Rozdział 1

Cel i zakres pracy

W pracy przedstawiono zagadnienia związane z zastosowaniem metod typu ‘pod prąd” do obliczeń przepływu w układach łopatkowych turbin cieplnych. Adaptacja do złożonych geom etrycznie układów w iąże się z pokonaniem licznych trudności, które w przypadku prostych geometrycznie obszarów przepływu nie występują.

Szereg złożonych zjawisk, jakie towarzyszą przepływom w turbinach, szczególnie w ostatnich stopniach turbin parowych, stanowi wyzwanie zarówno dla badań eksperymentalnych, jak i dla numerycznego modelowania. Wydaje się jednak, że badania eksperymentalne ze względu na skalę rozpatrywanego obiektu, jakim jest turbina, będą się koncentrowały na wycinkowym badaniu zjawisk i będą stanowiły w ażny element weryfikacji badań numerycznych. W modelowaniu numerycznym pozostaje nadzieja na poznanie struktury przepływu w całej maszynie wirnikowej z m ożliw ością szybkiej korekty geometrii kanałów oraz poszukiwania największych źródeł strat. Droga do osiągnięcia tego poziomu symulacji jest jednak jeszcze daleka. Podstawowe trudności, jakie spotyka się w przypadku modelowania num erycznego zjawisk przepływowych w turbinach cieplnych, to, po pierwsze, złożona geometria, stanowiąca trudność dla dyskretyzacji oraz dla określania warunków brzegowych. Składają się na nią wysokie łopatki, wąskie szczeliny m iędzyw ieńcow e, złożone kształty kanałów międzyłopatkowych oraz niekorzystne proporcje wymiarowe do konstrukcji siatek w przypadku łopatek długich. Po drugie, trudność stanowi obecność złożonych zjawisk przepływowych, powiązana z występowaniem przepływu w szerokim zakresie prędkości; od poddźwiękowych poprzez transoniczne do naddźwiękowych. Po trzecie, utrudnienie stanowi konieczność uwzględniania rzeczywistych własności czynnika roboczego oraz w ystępowanie dodatkowych zjawisk termodynamicznych, jak kondensacja pary wodnej w przestrzeni międzyłopatkowej. Prace skoncentrowane wokół tych problem ów stanowią podstawowy zakres tematyczny niniejszego opracowania.

1 W literaturze obcojęzycznej ten schemat określa się terminem ‘upwind’.

12 Rozdział 1

G łów ne cele, jakie starano się osiągnąć, podejmując problematykę zagadnień przepływow ych w turbinach cieplnych, można sformułować następująco:

• opracowanie własnych metod numerycznych do analizy procesu ekspansji gazu w kanałach przepływowych turbin cieplnych ,

• pokonanie trudności merytorycznych i obliczeniow ych w e wprowadzaniu w łasności gazu rzeczyw istego do analizy przepływowej,

• uw zględnienie w analizie wpływu zjawiska kondensacji, towarzyszącej ekspansji pary wodnej w ostatnich stopniach,

• określenie m ożliw ości badania struktury przepływu w w ielow ieńcow ych układach łop atk ow ych .

Prezentowane programy i algorytmy obliczeniow e opracowywane były przez szereg lat. W yniki uzyskiwane w ramach tych prac zawiera szereg publikacji zam ieszczonych w wykazie literatury. W niniejszym opracowaniu podano również rezultaty, które nie były publikowane. D otyczą one głów nie m odelowania przepływu w stopniu turbiny powietrznej, analizy przepływu na w locie do części niskoprężnej turbiny parowej oraz procesu m odelowania zjawiska kondensacji w kanałach łopatkowych w konfiguracji dwu- i trójwymiarowej. Ponieważ w Instytucie M aszyn i Urządzeń Energetycznych Politechniki Śląskiej wykorzystywano do m odelowania przepływu przez ostatnie stopnie turbin parowych metodę krzywizny linii prądu, starano się stworzyć efektywne narzędzia uzupełniające i rozszerzające aktualne m ożliw ości analizy i syntezy.

W rozdziale drugim pracy przedstawiono podstawowy opis matematyczny rozpatrywanego zadania mechaniki płynów, który uwzględnia specyfikę maszyn wirnikowych poprzez wprowadzenie bezw zględnego i w zględnego układu odniesienia odpowiednio dla kierownicy i wirnika. Podano w nim podstawowe zależności, opisujące przepływ, sformułowane w krzywoliniowym układzie współrzędnych, stanowiącym punkt wyjścia do konstrukcji rozwiązania i budowy algorytm ów numerycznych. Ważnym uzupełnieniem omawianych równań jest opis konstrukcji warunków brzegowych, zw iązków zamykających dla przepływu turbulentnego oraz równanie stanu.

R ozdział trzeci p ośw ięcony jest zagadnieniom związanym z dyskretyzacją obszaru przepływow ego. Konstrukcja i dobór siatki obliczeniow ej jest ważnym elementem procesu rozwiązywania numerycznego zadań dynamiki gazów.

W rozdziale przedstawiono różne warianty dyskretyzacji obszaru obliczeniow ego kanałów łopatkowych oraz om ówiono ich podstawowe cechy, decydujące o przydatności do rozwiązywania zagadnień związanych z m odelowaniem przepływu w m aszynach przepływowych.

W rozdziale czwartym opisano numeryczne sposoby rozwiązania równań zachowania. S zczególną uwagę zwrócono w nim na metody typu ‘pod prąd’, które są przydatne przy analizie przepływów transonicznych i naddźwiękowych. Om ówiono

Cel i zakres pracy 13

kilka wariantów tych metod i przedstawiono ich porównanie z punktu widzenia efektywności obliczeń i dokładności otrzymywanych rozwiązań. Zbadano efektywność schem atów w yższego rzędu dokładności w obliczeniach złożonych zjawisk falowych. Przedstawiono konstrukcję schematów numerycznych dla zagadnień wielowym iarowych oraz do obliczeń przepływu z zastosowaniem modelu gazu lepkiego. Szczególną uwagę zwrócono na sformułowanie warunków brzegow ych dla obliczeń palisadowych. Uwzględniono w tym przypadku szereg wariantów warunków brzegowych, występujących przy rozpatrywaniu zarówno izolow anych profilów, jak i układów łopatkowych stopnia osiowej maszyny wirnikowej. Ważnym elementem procesu numerycznego rozwiązywania zadania początkowo-brzegowego, decydującym o jego czasochłonności, jest dobór warunków początkowych. Zaproponowano różne warianty tych warunków i prześledzono ich w pływ na proces iteracyjny.

R ozwiązania pola przepływu dla szeregu rzeczywistych konfiguracji geom etrycznych turbin cieplnych przedstawiono w rozdziale piątym. Analizowano przepływ przez stopień turbiny powietrznej oraz stopień części niskoprężnej turbiny parowej dużej m ocy. Przeprowadzono również obliczenia różnych wariantów przepływu przez w ieniec łopatkowy, wykorzystując model gazu lepkiego.

Uzupełnieniem tych obliczeń jest przeprowadzenie numerycznej analizy warunków przepływu na w locie do części niskoprężnej turbiny średniej mocy. Przeprowadzone w tym rozdziale obliczenia uwzględniały przyjęcie założenia, że czynnik roboczy jest gazem doskonałym.

R ozdział szósty poświęcono modelowaniu przepływu pary wodnej w kanałach łopatkowych turbin parowych. W rozdziale tym czynnik roboczy, jakim jest para wodna, traktowany jest jako gaz rzeczywisty. Podano zależności do określania zw iązków m iędzy parametrami stanu dla pary wodnej, przydatne do obliczeń numerycznych. Wykorzystano w tym przypadku wirialne równanie stanu. Zwrócono uwagę na te aspekty, w których konieczne jest uwzględnienie w obliczeniach modelu gazu rzeczyw istego dla poprawnego określenia przepływu pary wodnej. Omówiono ważny element modelowania pracy w ieńców łopatkowych ostatnich stopni, związany z pojawieniem się procesu kondensacji pary wodnej. Uwagę skupiono na m odelowaniu procesu kondensacji spontanicznej. Przedstawiono sposób m odelowania przepływu czynnika dwufazowego oraz podano wyniki obliczeń przepływu z kondensacją w kanale turbinowym z wykorzystaniem różnych wariantów modeli wzrostu i tworzenia się kropli. Poddano analizie różne aspekty numerycznego modelowania procesu kondensacji, biorąc pod uwagę sposób dyskretyzacji obszaru przepływowego oraz korekty modeli fizycznych. Określenia struktury przepływu dwufazowego pary wodnej dokonano również dla przestrzennego kanału m iędzyłopatkowego kierownicy ostatniego stopnia turbiny dużej mocy.

W rozdziale siódmym przedstawiono koncepcję modułowej analizy przepływu w układach łopatkowych turbin cieplnych, natomiast rozdział ósm y zawiera podsumowanie przeprowadzonych badań.

W szystkim tym, którzy wykazywali zainteresowanie i nie szczędzili cennych uwag podczas wykonywania niniejszej pracy, a w szczególności Prof. Tadeuszowi Chm ielniakowi, a także recenzentom pracy Prof. Jerzemu Krzyżanowskiemu, Prof. Ryszardowi Grybosiowi i Dr. hab. Janowi Górskiemu, składam serdeczne podziękowania.

Rozdział 2

Matematyczny opis zasad zachowania

2.1. Ogólne równania zachowania

Przepływ substancji ciągłej jest opisywany przez równania Naviera-Stokesa, wyrażające zasadę zachowania masy, pędu i energii1. W przypadku gdy pominiemy

X

Rys.2.1. Obszar bilansowy Fig.2.1. Finite volume

oddziaływanie sił zewnętrznych oraz założymy, że w przepływie nie występują źródła ciepła, możem y zapisać całkową postać równań Naviera-Stokesa dla objętości V, opisującą niestacjonarny przepływ płynu ściśliw ego, lepkiego i przewodzącego ciep ło w następującej postaci:

Ń o - n d A = 0 (2.1)

V 0 1 A

1 Termin równania Naviera-Stokesa będzie używany do wszystkich równań zachowania, tzn. masy, pędu i energii, pomimo że ściśle odnosi się on jedynie do równania zachowania pędu.

16 Rozdział 2

gdzie

p

r

V.

Q = p v ^II

T^° p v v +

f

< e ,

v

y

W p ow yższych zależnościach Q o znacza w ektor zm iennych zachow aw czych, H0 rep rezen tu je uogólniony w ektor strum ieni przepływ ających p rzez pow ierzchnię A, a « j e s t w ektorem norm alnym do pow ierzchni A skierow anym na zew nątrz.

W p o w yższych zw iązkach e oznacza energię całk o w itą p rzy p a d ają cą na jed n o stk ę o bjęto ści, f rep rezen tu je tensor naprężeń, składający się z członu ciśnieniow ego o raz czło n u lepkiego:

f = p l + t (2 .2 )

nato m iast q je s t w ektorem strum ienia ciepła, który określa praw o Fouriera:

q = - X v T (2.3)

gdzie X je s t w spółczynnikiem przew odzenia ciepła.

B io rąc p o d uw agę to, że funkcja Q je s t ciągła w obszarze przepływ u, stosując tw ierdzenie G aussa-O strogradzkiego i przyrów nując w yrażenie podcałkow e m ożem y rów n an ie (2 . 1) zapisać w postaci różniczkow ej:

SQ -t

- J f + v - / ? o = 0 (2.4)

T a form a rów nań je s t fo rm ą zachow aw czą lub diw ergentną. R ów nania N aviera-S to k esa o p isu ją m ieszane zagadnienie fizyczne, konw ekcyjno-dyfuzyjne, a co za tym idzie, m atem atycznie p o sia d ają hiperboliczno-paraboliczny charakter.

W celu zam knięcia układu rów nań zachow ania m usim y zastosow ać dodatkow e zw iązki. Pierw szym z nich je s t rów nanie stanu, a drugim m odel turbulencji, który u zpełnia rów nania N aviera-Stokesa w p rzypadku przepływ u turbulentnego.

Z punktu w idzenia tw orzonych m odeli przepływ u w ygodne je s t rozdzielenie uogólnionego w ektora strum ieni H0 na część n i e l e p k ą / / i le p k ą H v:

H0= H + H V (2.5)

F orm y zarów no całkow a, ja k i różniczkow a rów nań zachow ania stanow ią punkt w yjścia do tw orzenia schem atów num erycznych. Z astosow any w niniejszej pracy schem at num eryczny bierze za p unkt w yjścia różn iczk o w ą postać rów nań zachow ania.

Matematyczny opis zasad zachowania 17

2.2. Równania zachowania Eulera

W przypadku gdy w rów naniach N aviera-Stokesa (2.1) lub (2.4) pom iniem y te człony, któ re uw zględniają lepkość oraz przew odzenie ciepła, otrzym am y postać rów nań zachow ania o k reślan ą ja k o rów nania Eulera. O p isu ją one przepływ czynnika n ielepkiego i nieprzew odzącego ciepła. Postać różniczkow a tych rów nań je st następująca:

- ^ 7 + v ■ f ) = 0 (2 .6)

Ot

W ektor strum ieni H w kartezjańskim układzie w spółrzędnych (x',y',zr) m a następujące składowe:

H = E \ X'+ F \ y ' + G \ z‘ (2.7)

W zw iązku z pom inięciem członów lepkich rów nania E ulera pozw alają na określenie na pow ierzchni ciała jed y n ie sił pochodzących od ciśnienia. W przypadku m odelow ania przepływ u je s t to często zadow alające przybliżenie szczególnie w tedy, gdy m am y do czynienia z przepływ em przy w yższych liczbach Reynoldsa.

R ów nania E ulera w postaci zależnej od czasu tw o rzą układ nieliniow ych rów nań różniczkow ych cząstkow ych typu hiperbolicznego, pierw szego rzędu. Postać równań zależna od czasu pozw ala na zachow anie tego sam ego typu rów nania różniczkow ego zarów no dla przepływ ów poddźw iękow ych, transonicznych i naddźw iękow ych.

P rz e z sw oje w łaściw ości polegające na tym , że m ożliw e je s t otrzym yw anie rozw iązań, w których param etry są ciągłe i nieciągłe, stanow ią one w ażne narzędzie w analizie zagadnień dynam iki gazów. W zagadnieniach przepływ ow ych, charakteryzujących się w ysokim i prędkościam i czynnika m o g ą się pojaw iać fale u d erzeniow e i rozrzedzeniow e, stanow iące nieciągłość param etrów przepływ owych.

M ożliw ość otrzym ania takiego rozw iązania rów nań E ulera zachodzi tylko w przypadku, gdy sform ułujem y je w tzw. form ie zachow aw czej. D la układu k artezjańskiego (x \y ',z r) m ożem y j e zapisać w rozw iniętej form ie następująco:

gdzie poszczególne w ektory m a ją postać (dla uproszczenia zapisu pom inięto nad w ektoram i kolum now ym i znak wektora):

c >

p

( >

pu p v r \

p w

pu p u 2 + p pu v pu w

p v £ = pvu F = p v 2 + p G - p vw

p w pwu p w v p w2 + p

. e > „ (\e + P > j

,

;

W p rzy p ad k u o bliczeń przepływ u w m aszynach w irnikow ych korzystne je s t tran sfo rm o w an ie ró w n a ń E ulera do w zględnego, kartezjańskiego układu w spółrzędnych, zw iązanego z w iru jącą ło p atk ą x,y,z, który obraca się w okół stałej osi x z p ręd k o śc ią k ą to w ą <3 = ( f i , 0 , 0). T ransform acja do układu w zględnego p o w o d u je pow stanie dodatkow ych członów źródłow ych w yrażających p rzy sp ie sz en ie C o rio lisa i przyspieszenie odśrodkow e. T ransform acja p rzep ro w ad zo n a je s t w ten sposób, że w układzie w zględnym w w ektorze k olum now ym Q składow e p rędkości w yrażone b ę d ą nadal w układzie bezw zględnym .

X

Rys.2.2. Oznaczenia układu współrzędnych Fig.2.2. Symbols of the coordinate system

Z w iązek m iędzy p o ch o d n ą p o czasie w ielkości skalarnej S w układzie b ezw zględnym i w w irującym u k ład zie w zględnym opisuje zależność:

( § L = ( f L +7-№>xi)>

D la w ielkości w ektorow ej V obow iązuje:

o

£2v

(2.10)

Z astosow anie tych zw iązków do rów nania (2.8) prow adzi do otrzym ania rów nań E ulera w kartezjańskim , w irującym układzie w spółrzędnych (x,y,z):

Ł Q + - ^ E ' + - ^ F ’ + 4z G ' + H= 0

dt dz (2.11)

/ *\

P ( pui ' f p vi ^{N V}

pu puu' + p puv' puw

p v E' = pvu' F' = p w ' + p G ‘ = pvw

p w pwu' pwv' pw w ' + p

. e . K e u + p u ' J ev' + p v ew' + p w ^

/ / = ( 0 , 0 , -pw C l, pvD., 0 ) 7

gdzie u',v',w' s ą prędkościam i w zględnym i

u' = u v ' = v + Q z w' = w - Cly (2 . 12)

P om iędzy poszczególnym i w ektoram i strum ieni w układzie w zględnym i b ezw zględnym istnieje następująca zależność:

E1 = E F' = F +C 1zQ G' = G - ClyQ (2.13)

U kład rów nań (2.11) je s t je d n ą z w ygodniejszych postaci rów nań zachow ania do obliczeń zagadnień zw iązanych z m odelow aniem przestrzennego przepływ u w m aszynach w irnikow ych.

2.3. Równania zachowania dla przepływu lepkiego

P on iew aż obliczenia przepływ u z zastosow aniem m odelu przepływ u czynnika lepkiego realizow ane były d la zagadnienia dw uw ym iarow ego, to dla zachow ania odpow iedniości w niniejszym punkcie przedstaw iona będzie postać tych rów nań zred u k o w an a do dw óch w ym iarów . O gólny układ rów nań N aviera-Stokesa (2.4) m ożem y ro zp isać w tym w ypadku do następującej formy:

20 Rozdział 2

dQ d ( E + E v) d ( F + F v)

dt 8x dy ⁼ 0 (2.14)

Q =

p pu p v

pu ■ E = ■ p u2 + p

■ F = ■ puv

p v puv p v2 + p

e u ( e + p ) v ( e + p )

0 0

E v = - ^T(xx) • F v = ■

T (xy) t(kk)

UZ(xx) + VT(jęy) + <J(X) _UT(*y) + VT(yy) + qty)

P rzy jm u jąc h ipotezę Stokesa, ok reśla ją cą zw iązek m iędzy w spółczynnikiem p rze w o d z en ia ciepła X i (j. jako:

3X + 2/u = 0

m ożem y zapisać składow e tensora naprężeń r następująco:

r = T(«) T(*y) l, Tto) TtoO J

(2.15)

(2.16)

*(xx) — — 3 n (2ux — vy )

T(xy) = —f l ( U y + Vx )

r ( y y ) = - jM (2vy - u x)

2.4. Normowanie wielkości obliczeniowych

W czasie przeprow adzania obliczeń num erycznych w skazane je s t posługiw anie się w ielkościam i bezw ym iarow ym i. Z astosow anie norm ow ania pow oduje, że używ ane w ielkości m a ją w artości zbliżonego rzędu. O trzym yw any w trakcie obliczeń b łąd zaokrągleń każdej z otrzym yw anych w ielkości je s t w tym w ypadku porów nyw alny. S tosow anie w ielkości bezw ym iarow ych ułatw ia rów nież porów nyw anie różnych przykładów obliczeniow ych.

P aram etry geom etryczne m ożna przykładow o odnieść do szerokości rozpatryw anego w ieńca łopatkow ego (lub długości cięciw y łopatki) /.

* = J ,

9= 7

Matematyczny opis zasad zachowania 21

P aram etry przepływ u są norm ow ane za p o m o cą gęstości spoczynkow ej po oraz prędkości dźw ięku a„. W ielkości stosow ane do opisu przepływ u lepkiego norm ow ane są za po m o cą odpow iadających im w ielkości odniesienia. N orm ow anie m oże się odbyw ać w dow olny sposób. N ajistotniejsze je s t zastosow anie odpow iedniej konsekw encji w doborze param etrów odniesienia. Jednym ze sp o so b ó w norm ow ania, który przyjęto w trakcie realizacji obliczeń, jest:

u V - w

u ~ ao v ~ a o W — ao

II

sh *

Ci = D.1

P — 2

Poa0 ao

- ta o

.

t = ~ T p o

(2.18)

L iczb y kryterialne R eynoldsa i P randtla s ą definiow ane przez odpow iednie w ielkości następująco:

D laopo D cpp.

x* = -jir> ?r=~T (219)

P raktycznie w dalszych rozw ażaniach w szystkie w ielkości należy rozum ieć jako w ielkości bezw ym iarow e. D la uproszczenia kreska nad w ielkością bezw ym iarow ą zostanie pom inięta.

2.5. Krzywoliniowy układ współrzędnych

D o obliczeń przepływ u w kanałach lub w okół ciał o złożonej geom etrii, kiedy to linie siatki w układzie kartezjańskim nie p okryw ają się z brzegiem , w ygodnie je st do k o n ać transform acji układu rów nań zachow ania do krzyw oliniow ego układu w spółrzędnych. R ozw iązyw anie zagadnień przepływ ow ych w krzyw oliniow ym u k ładzie w spółrzędnych pozw ala n a łatw iejsze postaw ienia w arunków brzegow ych, gdyż najczęściej w tedy linie siatki są ortogonalne do brzegu. S zczególnie m a to znaczenie p rzy obliczeniach przepływ u lepkiego, kiedy dopasow anie odległości pom iędzy liniam i siatki przy ściance do przew idyw anych gradientów param etrów przep ły w u je s t bardzo ważne.

T ransform acja rów nań zachow ania z kartezjańskiego układu w spółrzędnych do d ow olnego krzyw oliniow ego układu w spółrzędnych je s t prezentow ana szeroko w literaturze (np. H irsch 1990, Chm ielniak 1994).

W w yniku transform acji kartezjańskiego układu w spółrzędnych (x,y,z) do n ow ego krzyw oliniow ego układu w spółrzędnych ( |,T |,Q obszar fizyczny zostaje zastąpiony p rze z obszar o regularnych oczkach, który je s t tak dobrany, że pokryw a się z liniam i siatki, a odstęp pom iędzy liniam i wynosi A£ = A// = AC = 1 (rys.2.3).

Rys.2.3. Transformacja układu współrzędnych Fig.2.3. Transformation of the coordinate system

R ów nania zachow ania s ą transform ow ane do now ego krzyw oliniow ego układu w spółrzędnych, w którym s ą rozw iązyw ane.

Z astąp ien ie o peratorów różniczkow ania w przestrzeni p rzeprow adza się z w ykorzystaniem reguły łańcuchow ej, która dla ogólnej transform acji w postaci:

p o zw ala n apisać zależność:

Ś = £ ( t , x , y , z ) tj = r i ( r , x ,y ,z ) C = ( ( r , x , y , z )

T - t

(2.20)

f JL f 9 )

dx

JL ' f* 7* ęx ' di dy — iy Cy JL dr/

JL k Śz flz Cz d

\ dz , l 3C J

(2.21)

N ieznane człony m etryczne otrzym uje się p rzez odw rócenie m acierzy Jakobiego p rzekształcenia odw rotnego i zastosow anie w yznacznika:

J = l/y - 1 (2 .22)

Z w iązki m etryczne m a ją w tedy postać:

= J(ynZ( - Z,iy() >lx = Ązęyt; - y p t ) Cx = J ( y ^ n- Z tf n)

t y = J ( z , x c - x , z ( ) tjy = Ą x i z i ; - z i x c) Cy=J(z(x n - x i z n) (2 .2 3 ) i z = Ą x ny ę - y n xÓ n* =

J(y(xC - x<yÓ

J(xt y n- y ( xi)

W p rzypadku dyskretnym w ielkość J ' odpow iada objętości rozpatryw anej kom órki elem entarnej, a w yrażenia (2 .2 3 ) p rez en tu ją składow e w ektorów norm alnych do ścianek kom órki. M oduł w ektora norm alnego ścianki kom órki odp o w iad a p o lu pow ierzchni tej ścianki.

2.5.1. T ran sform acja rów nań E ulera

Z astosow anie reguły łańcuchow ej do rów nań zachow ania E ulera pozw ala otrzym ać p o transform acji zachow aw czą form ę tych rów nań w postaci:

<2 - 2 4 >

z poszczególnym i w ektoram i kolum nowym i:

Q = J l (p ,p u ,p v ,p w , e ) T

' P i f ' p V ' p W

p u l f + £ xp p u V + rjxp puW* + CxP

p v l f + iy P F - J ~ { p\>V + t]yp G = J - ‘ p v W + CyP

p w i f + i z p p w V + t]zP p w W 1 + ( zp

v e l f + p U j , e V + p V J e W + p W

H = J - ' ( 0 ,0 , -pw C l, p v d , 0) T

gdzie kontraw ariantne prędkości w zględne w ynoszą:

U1 = £ xu + i y(v + O z) + j(w - Oy)

V = t]xit + tjy(v + O z) + tj2(w - Q y) (2.25)

W 1 = CxU + Cy(v + O z) + C2(w - £2y)

a kontraw ariantne prędkości bezw zględne

U = £ xu + śyV + ę zw

V = t ] xu + t]y v + t] zw (2.26)

W = C x U + Cr v + ę2w

24 Rozdział 2

W p ro w ad zając do rów nania zachow ania energii w yrażenie na rotalpię (R + p ), gdzie:

R = j P + y ( w 2 + v 2 + w 2) - pęi(yvy - vz) (2.27)

m ożem y rów nania zachow ania (2.24) zapisać w form ie stosow anej często w o b liczen iach zagadnień przepływ ow ych w m aszynach w irnikow ych, w której w ektory kolum now e m a ją postać:

Q - J ^ ( p , p u , pv, pw , R) T

E - J ~ ]

' p i f p v ' P W '

p u l f + i xp p u V + t]xp p u W + ę xp

p v l f +ŹyP F = J~l p v V +rjyp G = J -' p v W + (yP

p w U1 + (,zp P w V + t]zp p w W + CzP

{ ( R + P ) l f ,

.

J

(2.28)

H = ( 0 ,0 ,-p w C l,p v ę i, 0 ) r

N ale ży zw rócić uw agę na fakt, że w szelkie zw iązki m etryczne otrzym ane w w yniku transform acji do układu w zględnego nie z a le żą od czasu, co pow oduje u pro szczen ie obliczeń.

2 .5.2. T ran sform acja rów nań N aviera-Stokesa

D la przepływ u dw uw ym iarow ego postać zachow aw cza rów nań N aviera-Stokesa m o że być zapisana w form ie bezw ym iarow ej, w krzyw oliniow ym układzie w spółrzędnych (4 ,t|) następująco:

M . M dF _

ar + l ź + -d j = Re

' dE v dF, l d i + dr,

g dzie transform ow ane strum ienie w yznacza się ze zw iązków :

E = J - {( i xE + £ yF ) F = J ~ l (tjxE + t]y F) E v ~ J (£XE V + i y F v) F v — J (tjxE v + i)vF v)

(2.29)

(2.30)

Po w prow adzeniu prędkości kontraw ariantnych U, V poszczególne w ektory kolum now e w ynoszą:

Matematyczny opis zasad zachowania_________________ 25

P p U p V

pu

p v E = J -' p u U + £ xp p v U + i yp

n p u V + tjxp

p v V + t] yp ■ (2.31)

. * . U(e + p ) V(e + p )

0 0

£ x ? ( x x ) + i y T (xy)

F V= J ~ X- t]x*(xjc) + riyX (xy)

£ x T (x y ) + i y t f y y ) t1xT(xy) + t f y t f y y )

i x E ev + £ y F ev t]xEev + rjyF ev

W ielkości E ev,F ev s ą składow ym i w ektorów E, F (rów nanie 2.14), dotyczącym i rów nania zachow ania energii. S kładow e tensora naprężeń m a ją po transform acji n a s tęp u jącą postać:

r (xx) = /i( 4 ( £ xu t + rjxun) - 2(źyV( + //yv ,))/3

T(xy) = p((Ś yU i + tlyitn + ę xv ( + rjxVq) (2.32) T(hv) = P(~ 2 (Ź*u( + Tl*un) + 4(ŚyV( + riyVnW 3

Z g o d n ie z h ip o te zą B oussinesąa (1897), zw iązek m iędzy tensorem naprężeń turbulentnych a tensorem prędkości deform acji, analogicznie ja k w przepływ ie lam inam ym , określa się za p o m o c ą w spółczynnika zwanego lepkością turbulentną.

E fektyw ny w spółczynnik lepkości składa się w tym przypadku z dw óch części:

H = p i + p t (2.33)

W a rto ść w spółczynnika lepkości lam inam ej w yznacza się z form uły Sutherlanda, któ ra dla pow ietrza m a postać:

l + 110.4/7o T + 110.4/7'o

Pi = Tm (2-34)

g dzie To o znacza tem peraturę spiętrzenia na wlocie.

W artość w spółczynnika lepkości turbulentnej p, w yznaczana je s t na podstaw ie je d n e g o przyjętego m odelu turbulencji. Jeżeli podobnie w prow adzim y turbulentną liczbę P randtla, to efektyw ny w spółczynnik przew odności cieplnej, uw zględniający zw ię k sz o n ą w ym ianę ciepła w w arstw ie turbulentnej, m ożem y określić przyjm ując:

( A S ] _ M L + J £ i (2 35 )

I P r J " Pr/ Pr, 1 ;

i składow e strum ienia ciepła m a ją w tedy postać:

< ? M - ( p r ) K l i d x T ( + rixTn)

r U \ 1 (Z .ó O )

<?(v) = [ p f ) K _ \ (£ yTi + ly ^ n )

D la p o w ietrz a p rzyjm uje się następujące w artości lam inam ej i turbulentnej liczby P ra n d tla /V ,=0.74, Pr,=0.9.

2.6. W arunki brzegowe

R ów n an ia N aviera-S tokesa z m atem atycznego punktu w idzenia s ą typu m ieszanego hiperboliczno-parabolicznego, natom iast rów nania E ulera są typu h iperbolicznego. R ozw iązanie tych rów nań w ym aga określenia skończonego obszaru całkow ania, a co się z tym w iąże - sprecyzow ania w tym obszarze w arunków p o cz ątk o w y c h i zdefiniow ania na je g o brzeg ach w arunków brzegow ych.

S form ułow anie odpow iednich w arunków na brzegu obszaru całkow ania, zap ew n iający ch fizycznie p o p raw n ą w ym ianę inform acji m iędzy obszarem c a łk o w an ia i obszarem zew nętrznym stanow i podstaw ę otrzym ania w łaściw ego ro zw ią zan ia rów nań zachow ania. Z drugiej strony w arunki brzegow e m u sz ą być tak skonstruow ane, aby w przypadku zagadnień, w których w ystępuje ograniczenie co do w ielkości o bszaru obliczeniow ego, nie generow ały niefizycznych zaburzeń pola przepływ u.

W p rzy p ad k u układu rów nań różniczkow ych typu hiperbolicznego, ja k i tw o rzą n iestacjo n arn e rów nania Eulera, m ożem y w yróżnić charakterystyczne kierunki ro zp rze strzen ian ia się zaburzeń. K ierunki te odpow iad ają w artościom w łasnym m acierzy Jacobiego otrzym anym po linearyzacji rów nań zachow ania. P rzykładow o d la kieru n k u % s ą to w artości w łasne m acierzy A = Identyfikacja kierunków

dQ

rozp rzestrzen ian ia się zaburzeń m a kluczow e znaczenie w określaniu w arunków brzegow ych n a pow ierzchniach ograniczających obszar obliczeniow y. P oniew aż o b szar obliczeniow y je s t transform ow any do pom ocniczego regularnego obszaru obliczeniow ego, w którym osie w spółrzędnych s ą ortogonalne do brzegów , identyfikację odpow iednich w arunków brzegow ych m ożna oprzeć na kierunku zgodnym z o sią krzyw oliniow ego układu w spółrzędnych. Jeżeli rozpatrzym y kierunek £,, to rów nanie zachow ania m ożem y zapisać w następującej quasi-liniow ej postaci:

dl 3 f ( ^

M a cie rz p o sia d a rzeczyw iste w artości w łasne oraz liniow o niezależne lew e i praw e w ektory w łasne. W zw iązku z tym m acierz A m ożna diagonalizow ać w następ u jący sposób:

(2.38)

gdzie A (o = d ia g (ł(Q i,..., A ^ s ) je s t m a cie rz ą d iagonalną w artości w łasnych, a L (() je s t m a c ie rz ą lewych w ektorów w łasnych m acierzy A . M nożąc rów nanie (2.37)

lew ostronnie p rzez L* i w prow adzając zależność (2.38) m ożem y zapisać:

dQ dQ

L (.0 ~df + L (i)L h A(0L (0 ~ ^ = 0 (2.39)

P o w prow adzeniu now ych zm iennych tzw. zm iennych charakterystycznych zdefiniow anych jako:

dW=L(i)dQ ^(2.40)

m ożem y rów nania zachow ania zapisać w postaci charakterystycznej:

d W d W n ^{,7 d n}

-qT + = o ( )

P ostać zm iennych charakterystycznych łatwiej je s t otrzym ać, w ykorzystując zam iast w ektora zm iennych zachow aw czych Q w ektor zm iennych prostych p = J ~ l (p, u, v, w , p ) T . W iąże się to z konieczn o ścią m odyfikacji m acierzy w ektorów w łasnych. M am y wtedy:

A o = L (0 dP_

dQ

W artości w łasne, tw orzące A (^ , p o siad ają tyle jednakow ych w artości, ile w ym iarów m a rozpatryw ane zagadnienie. O ile w przypadku zagadnienia jednow ym iarow ego istnieje w zw iązku z tym je d n a postać praw ego i lewego w ektora w łasnego, o tyle w przy p ad k u dwu- i trójw ym iarow ym m ożliw e s ą różne sposoby definicji.

W ykorzystując definicję p o d an ą przez M anna (1992):

0 - I * - l y - 6 1 l(pc)

1 0 0 0 - l / c2

o o - m i + z b Ш 2> + & ) 0 0 1 - 6 4 / ( 3 + £ ) - Ш ( ¥ у + & ) 0

О 6 l y Ь 1 ! { р с )

(2.43)

gdzie

<fw = ^w /K I

28 Rozdział 2

m o żem y w konsekw encji otrzym ać następujące postacie zm iennych charakterystycznych

d W = Ł ii)d P = J -'

% - d l A H |

df>

°P a2

-16/(3+62>dv+1,/(3 + ~e2)dw

Su - 66/(3+^)5v - 66/(3 +&№ » Js- + dL/|vcj;|

(2.45)

R ów nania zachow ania w postaci (2.41) tw o rzą układ rów nań różniczkow ych, które o b o w ią z u ją w zdłuż krzyw ych charakterystycznych.

dt

K ierunki te, określone ze zw iązku:

(2.46)

De/(/4 - -!({),/) = 0

p rz y jm u ją następujące w artości:

A«)i =W-aJiT+lf+if

^ (0 2 ,3 ,4 = C /

^(05 = l / + a j l T + ~ £ j + ę J

(2.47)

(2.48)

gdzie:

l f = ixU + ź y (v + Clz) + <f2(łv - Q y)

Z m ienne charakterystyczne w zdłuż odpow iadających im charakterystyk są stałe.

W ielkości W s ą określane ja k o niezm ienniki R iem anna. Przy w ykorzystaniu definicji prędkości dźw ięku:

a2 = (2.49)

w w yrażeniach (2.45) niezm ienniki R iem anna m ożna zapisać dla ro zprzestrzeniającego się zaburzenia, odpow iadajacego w artości w łasnej X{i)\ w postaci:

Matematyczny opis zasad zachowania 29

R i = p/QK = Const

Ri = -6/(3 + ^)v + 6/(3 + &)w = const

Rą — u —

66/(3

66/(3 +

rC 1

a dla rozprzestrzeniającego się zaburzenia, odpow iadającego w artości własnej w postaci:

R] = U\v£I - ~ ~ J a = const

R2= p /Q K = const ^ ^ (25Q a)

R i = - 6 / ( 3 + 3 ) v + 6 / ( 3 + 3 ) w = const Rą ~ u — 6 6 / ( 3 + 6 2)v - 6 6 / ( 3 + 6 )w = cons/

W artości w łasne A(f) określają kierunek przepływ u inform acji m iędzy obszarem obliczeniow ym a obszarem zew nętrznym . Liczba w arunków brzegow ych je st zw iązana z lic zb ą charakterystyk, w zdłuż których obszar zew nętrzny kom unikuje się z w nętrzem .

B rzeg obszaru obliczeniow ego m ożem y podzielić na trzy podstaw ow e podobszary: w lot, w ylot i ścianka nieprzenikalna.

W lot

W p rzypadku w lotu prędkość norm alna do brzegu je s t skierow ana do w nętrza obszaru. Jeżeli założym y, że w lot przyjm iem y dla £, = O, to m ożem y zapisać:

A(i)2,3,4 = t / > 0 (2.51)

i m ożem y m ieć do czynienia z następującym i przypadkam i:

N apływ p o ddźw iąkow y

G dy n apływ je s t poddźw iękow y, czyli:

A(ć)i = U1 + £ y + £z < O (2.52)

otrzym ujem y konfigurację charakterystyk przedstaw ioną na rys.2.4a. W idać z niej, że cztery z p ięciu charakterystyk skierow ane s ą do w nętrza obszaru i n io są ze so b ą inform acje z zew nątrz. W zw iązku z tym konieczne się staje sprecyzow anie czterech fizycznych w arunków brzegow ych na w locie. Jeden w arunek zw iązany z charakterystyką A ^ s określa się na podstaw ie param etrów z w nętrza obszaru.

N a p ływ n addźw iękow y

W p rzy p a d k u gdy spełniona je s t nierów ność:

A (f)i

= V - aJÎT T ïf+ ÎJ >

Rys.2.4. Konfiguracja fal na wlocie: a. wlot poddźwiękowy, b. wlot naddźwiękowy Fig.2.4. Waves configuration at the inlet: a. subsonic inlet, b. supersonic inlet

W y lo t

W przypadku w ylotu prędkość norm alna do brzegu je s t skierow ana na zew nątrz obszaru obliczeniow ego. Przy założeniu że w ylot je s t na linii £ = £ max> m ożem y zapisać:

A({)2,3,4 = l f > 0 (2.54)

P o dobnie ja k to m iało m iejsce dla w lotu, rozróżniam y dw a przypadki:

W ypływ pod d źw ięk o w y

Jeżeli spełniona je s t nierów ność

% ) i a j i î + i y + i z < 0 (2.55)

to na w ylocie m am y do czynienia z układem charakterystyk przedstaw ionym na rys.2.5a. W ynika z niego, że tylko w zdłuż jed n ej charakterystyki zaburzenie dociera

do w nętrza obszaru. K onieczne je s t w ięc postaw ienie jed n eg o w arunku fizycznego na brzegu.

W ypływ naddźw iękow y W przypadku gdy

% )! = v + ł l > 0 (2.56)

co przedstaw ione je s t na rys.2.5b, w szystkie zaburzenia rozp rzestrzen iają się w kieru n k u obszaru zew nętrznego; to w konsekw encji oznacza, że nie potrzeba zadaw ać żadnych w arunków brzegow ych na wylocie.

Rys.2.5. Konfiguracja fal na wylocie: a. wylot poddźwiękowy, b. wylot naddźwiękowy Fig.2.5. Waves configuration at the outlet: a. subsonic outlet, b. supersonic outlet

Ś cian ka n iep rzen ikalna

N a ściance w przypadku przepływ u nielepkiego zadaje się w arunek nieprzenikalności, tzn. składow a norm alna prędkości do ścianki musi być rów na zero lub, co je s t rów now ażne, prędkość płynu m usi być styczna do ścianki. O bow iązuje w ięc dla ścianki n a linii siatki r\ - 1 warunek:

A(,)2,3,4 = V = 0 (2.57)

P o zo stałe charakterystyki m a ją kierunki określone z następujących związków :

Aft)! - - a j n i + f ly + n t A (,)s= a j t j ł + rfi + rił

U kład charakterystyk dla ścianki nieprzenikalnej przedstaw ia rys.2.6.

32 Rozdział 2

Rys.2.6. Konfiguracja fal na ściance nieprzenilcalnej Fig.2.6. Waves configuration at the wali

W arunki brzegow e zostały podane na w ybranych płaszczyznach obszaru obliczeniow ego. W przypadku gdy brzeg usytuow any je s t ortogonalnie do innych k ie ru n k ó w krzyw oliniow ego układu w spółrzędnych, należy zam ienić prędkości i czło n y m etryczne transform acji odpow iednim i w ielkościam i w rozw ażanym kierunku.

2.7. Związki zamykające dla przepływu turbulentnego

M od ele turbulencji, stosow ane w przypadku m odelow ania przepływ ów z w ysokim i liczbam i R eynoldsa p rzez kanały o skom plikow anych kształtach, op ierają się w w iększości na uśrednionych po czasie rów naniach N aviera-Stokesa. R ów nania te o d ró żn ia ją się od nieuśrednionych rów nań N aviera-S tokesa w prow adzonym tensorem naprężeń R eynoldsa oraz turbulentnym strum ieniem ciepła. T e dodatkow e człony w rów naniu N aviera-Stokesa aproksym ują czasow o uśrednione fluktuacje p aram etrów przepływ u, w ykorzystując odpow iednie zw iązki zam ykające tzw.

m odele turbulencji. Zw iązki te najczęściej obow iązu ją dla konkretnych w arunków przepływ ow ych, gdyż ja k dotąd nie stw orzono m odelu uniw ersalnego. N ajprostsze zw iązki zam ykające op ierają się na hipotezie B oussinesąa (1897), która zakłada zw iązek m iędzy tensorem nap rężeń turbulentnych a tensorem prędkości deform acji za p o m o c ą lepkości turbulentnej, analogicznie ja k to m a m iejsce w przepływ ie lam inam ym . D o określenia w artości lepkości turbulentnej w ykorzystuje się szereg m odeli turbulencji. W zależności od typu rów nania zastosow anego do określenia lepkości turbulentnej rozróżnia się m etody: algebraiczne, je d norów naniow e i dw urów naniow e. M odele algebraiczne oparte na hipotezie "drogi m ieszania"

P ran d tla (1926) znalazły szerokie zastosow anie w obliczeniach przepływ u w kan ałach m aszyn przepływ ow ych. M odel algebraiczny opracow any przez B oldw ina i L om axa (1978) w ykorzystyw any je s t szczególnie często do analizy przepływ u w

Matematyczny opis zasad zachowania 33

złożonych geom etriach układów łopatkow ych (np. M erz 1995, C hm ielniak i W róblew ski 1995). W yróżnia się on tym , że do obliczenia w spółczynnika lepkości nie w ykorzystuje grubości w arstw y przyściennej. M odel ten, adaptow any z teorii w arstw y przyściennej, zakłada istnienie dw óch w arstw: w arstw y w ew nętrznej i w arstw y zew nętrznej. O bliczenia w w arstw ie turbulentnej od b y w ają się w zdłuż w spółrzędnej y norm alnej do ścianki, k tó rą w przypadku ortogonalnego układu w spółrzędnych m ożem y określić jako:

y Ą { ( ^ f > W di <2-59»

W spółczynnik lepkości turbulentnej określony je s t w tedy wzorem :

| ^{(Mi)wew y<ygr n}

= \ i \ ^ (2.60)

[ (Mt)zew y ^ y g r

gdzie y v o znacza najm n iejszą w artość y , dla której obow iązuje:

^Mt^wew ” ^Mt\ew (2.61)

W arto ść w spółczynnika lepkości w w arstw ie wew nętrznej określa się z zależności P randtla-V an D riesta:

(M,)w m = R e p l2\n \ (2.62)

gdzie Q ozn acza w irow ość zdefiniow anąjako:

(2 -63)

Z ależność (2.63) w krzyw oliniow ym układzie w spółrzędnych m ożna zapisać w postaci:

O = + tfxVq - iyUę - r,yUq) (2.64)

W artość drogi m ieszania / w (2.62) określa się w funkcji param etru y z zależności van D riesta (1956):

l = k y ( l - e (- ^ > ) (2.65)

przy przy jęciu bezw ym iarow ej odległości od ścianki:

W arto ść n ap rężeń na ściance odpow iada w artości naprężeń stycznych określo n y ch d la kom órki bilansow ej przy ściance, natom iast w spółczynnik lepkości p rzy zało żen iu podw arstw y lam inam ej przy ściance określa się z zależności (2.62).

D la w arstw y zew nętrznej w spółczynnik lepkości turbulentnej w yznacza się z zależności:

(M i\ew = R e CcipFwakeFki (2.67) g dzie w spółczynnik:

F W A K il — m in(y m a x / ł m a x » F m a x ) (2.68) w ym aga znajom ości wartości_>w . W yznacza się j ą j a k o argum ent funkcji:

F O O = > 4 f i l( l- e < - '+M,>) (2.69)

d la któ reg o funkcja m a sw oje m aksim um . W ielkość ud oznacza m aksym alną różnicę p rędkości w obszarze w arstw y przyściennej.

W spółczynnik korekcyjny K lebanow a wynosi:

F K l= l / ( + 5 . 5 ( ^ ) 6) (2.70)

a odpow iednie stałe p rzy jm u ją w artości: CKi = 0.3, A* =26, Ccl = 0.01688, £=0.41, C„*= 0.25.

2.8. Równanie stanu gazu

U kład rów nań zachow ania (2.4) je s t niekom pletny bez zdefiniow ania dodatkow ych zw iązków zam ykających. Z ależności m iędzy param etram i term odynam icznym i i kalorycznym i czynnika określa się z term icznego i kalorycznego rów nania stanu. R ów nania stanu s ą w ystarczające do zam knięcia u kładu rów nań zachow ania w przypadku rozpatryw ania przepływ u nielepkiego, a k o nieczne w p rzypadku rozpatryw ania p rzepływ u lepkiego.

R ów nania stanu p o zw a lają określić param etry term odynam iczne i kaloryczne w funkcji zm iennych zachow aw czych param etrów przepływ u przy założeniu stanu rów now agi term odynam icznej. D la w artości energii całkow itej m ożem y zapisać:

e = pE = p (E w + ^ - ) (2.71)

gdzie E v przed staw ia energię w ew nętrzną, któ ra m oże być w yrażona przez k alo ry czn e rów nanie stanu dla kalorycznie i term icznie doskonałego gazu w postaci zależności:

E w = c vT (2.72)

T erm iczn e rów nanie stanu p ozw ala n a w yznaczenie ciśnienia n a podstaw ie dw óch zm iennych. W przypadku term icznie doskonałego gazu m ożem y zapisać zw iązek:

p = pR T (2.73)

a przy uw zględnieniu (2.72) otrzym am y związek:

p = ( x - l )p E v (2.74)

co p ro w a d zi do zapisu rów nania (2.71) w postaci:

e = p E = 1£ T + f Ą - (2.75)

R ów nania stanu oraz w ykorzystanie zw iązków term odynam icznych pozw alają o kreślić inne w ym agane w obliczeniach w ielkości (np. prędkość dźwięku).

Rozdział 3

Obszar obliczeniowy i sposoby jego dyskretyzacji

3.1. W prowadzenie

Z agadnienia m echaniki płynów są opisyw ane nieliniow ym i rów naniam i różniczkow ym i, które najczęściej nie d a ją się rozw iązyw ać w sposób analityczny.

W tym przy p ad k u konieczne je s t zastosow anie dyskretnego sposobu rozw iązania przy w ykorzystaniu m etod różnego typu, ja k m .in. m etody różnic skończonych, o b jętości skończonych, elem entów skończonych. W ym ienione m etody w y d ają się być obecnie najbardziej rozpow szechnione i skuteczne w obliczeniach przepływ ow ych. W ym agają one określenia w ybranych p unktów obszaru fizycznego, dla któ ry ch konstruuje się rozw iązanie. W tym celu generuje się siatki num eryczne, k tó re w sw oich w ęzłach tw o rzą zbiór p unktów reprezentatyw nych do rozw iązania cząstkow ych rów nań różniczkow ych m echaniki płynów w rozpatryw anym zagadnieniu.

W o statnich latach w raz z postępem m etod obliczeniow ych w m echanice płynów obserw uje się dynam iczny rozw ój m etod generacji siatek num erycznych. M ożna m ów ić ju ż o w ykształtow aniu się autonom icznej dziedziny nauki z zakresu m atem atyki stosow anej - generacja siatek num erycznych (T hom pson i W arsi 1985).

Sposób przyjętej dyskretyzacji, czyli generacji siatki, m a istotny w pływ na ja k o ść otrzym yw anych rezultatów obliczeń. Istotnym elem entem przy generow aniu je s t takie dostosow anie siatki do rozpatryw anego zadania obliczeniow ego, aby otrzym ać popraw ne rezultaty. Przy konstruow aniu siatek num eiycznych dąży się do zgodności w ęzłów siatki z brzegiem obszaru obliczeniow ego, co pozw ala dokładnie postaw ić w arunki brzegow e, w przeciw nym razie należy się liczyć z m ożliw ością pojaw ienia się błędu w rozw iązaniu num erycznym rozw ażanego zagadnienia brzegow ego.

P oniew aż nie istnieje uniw ersalny przepis na stw orzenie siatki obliczeniow ej p ow stało w iele typów siatek. K ażda z siatek m a zarów no zalety, ja k i wady.

W obliczeniach m aszyn przepływ ow ych stosow ano, począw szy od prostej dyskretyzacji, polegającej na p odzieleniu obszaru na regularne o czka kw adratow e (K atsanis 1969), p oprzez regularne siatki dopasow ane do kształtu ciała o różnej

Obszar obliczeniow y i sposoby jego dyskretyzacji 37

strukturze oczek, w których linie m a ją gładki przebieg, do siatek nieregularnych o oczkach trójkątnych lub czw orokątnych (np. Irm ish 1994). W siatkach nieregularnych, w przeciw ieństw ie do siatek regularnych, punkty siatki nie są identyfikow ane z liniam i w spółrzędnych, lecz są num erow ane indyw idualnie w pew nym ustalonym porządku. Siatki te o d zn aczają się w ięk szą elastycznością w p rzypadku skom plikow anych geom etrii (m ają one także cechę łatw ego sam oadaptow ania się do zjaw isk panujących w przepływ ie). Z tego pow odu s ą coraz częściej stosow ane w zagadnieniach num erycznej m echaniki płynów , choć lokalny układ w spółrzędnych kom plikuje znacznie algorytm rozw iązania.

N ajw ażniejszym obszarem rozw oju m etod generacji siatek je st rozw ój system u dynam icznej adaptacji, w którym punkty siatki p rzem ieszczają się podczas obliczeń przepływ ow ych zależnie od w ystępujących zjaw isk fizycznych. G eneracja siatki w tym przy p ad k u odbyw a się rów nocześnie z obliczeniam i przepływ owym i.

Z astosow anie tych siatek pozw ala osiągnąć d u żą rozdzielczość w m odelow aniu n ieciągłości w przepływ ie. W ad ą tego typu m etod je s t złożony algorytm obliczeń i znacznie w ydłużony czas rozw iązania zagadnienia, co istotnie ogranicza m ożliw ości ich stosow ania.

W niniejszej pracy skupiono uw agę na zagadnieniu generacji siatek regularnych, najbardziej obecnie rozpow szechnionych w analizie przepływ u w m aszynach w irnikow ych.

3.2. Transformacja obszaru oraz podział metod generacji siatek

Z w ielu zagadnień fizyki m atem atycznej w ynika, że rozw iązanie danego p roblem u znacznie się upraszcza, gdy zastosujem y transform ację układu w spółrzędnych. Podstaw ow ym dążeniem przy transform acji je s t odw zorow anie b rzeg ó w obszaru fizycznego na obszar transform ow any, w którym linie układu w spółrzędnych pokryw ają się z brzegam i. T en zabieg znakom icie ułatw ia form ułow anie w arunków brzegowych.

P odczas generacji siatek num erycznych dokonuje się często transform acji złożonych fizycznych obszarów przepływ u do obszarów obliczeniow ych prostokątnych. S zczególnie korzystne je s t transponow anie geom etrii kanałów m aszyn przepływ ow ych i zadaw anie w arunków n a brzegach regularnych.

P od staw o w ą c e ch ą takiego przekształcenia je s t przyporządkow anie brzegom w rozpatryw anym obszarze fizycznym linii w spółrzędnych w przekształconym obszarze. Rys.3.1 obrazuje przekształcenie obszaru fizycznego zadanego w w spółrzędnych kartezjańskich (x,y) na obszar obliczeniow y w e w spółrzędnych krzyw oliniow ych (<!;, //).