VYSOKÁ ŠKOLA BÁŇSKÁ - TECHNICKÁ UNIVERZITA
OSTRAVA
Modelování proudění tekutin
FLUENT, CFX
Milada Kozubková
OSTRAVA 2008
i
Obsah
Obsah
SEZNAM POUŽITÝCH OZNAČENÍ ............................................................................................ VI
1.
ÚVOD DO NUMERICKÉHO MODELOVÁNÍ PROUDĚNÍ.................................................1
1.1.
ÚVOD DO PROUDĚNÍ TEKUTIN ...........................................................................................1
1.1.1.
Dělení podle fyzikálních vlastností tekutiny ..............................................................1
1.1.2.
Dělení podle kinematických hledisek.........................................................................2
1.2.
PŘENOS HMOTY, HYBNOSTI, TEPLA PŘI NEIZOTERMNÍM PROUDĚNÍ NESTLAČITELNÉ
TEKUTINY ........................................................................................................................................3
1.3.
PARCIÁLNÍ DIFERENCIÁLNÍ ROVNICE A OKRAJOVÉ PODMÍNKY...........................................5
1.4.
OKRAJOVÉ PODMÍNKY PROUDĚNÍ NESTLAČITELNÉ TEKUTINY ...............................................7
1.4.1.
Typy okrajových podmínek.........................................................................................7
1.4.2.
Podmínky vstupu a výstupu. ......................................................................................7
1.5.
ŘEŠENÍ ROVNICE VEDENÍ TEPLA .......................................................................................9
1.6.
ŘEŠENÍ PROUDĚNÍ PŘI NÁHLÉM ROZŠÍŘENÍ PRŮŘEZU .................................................... 12
2.
TVORBA GEOMETRIE A VÝPOČTOVÉ SÍTĚ ................................................................ 16
2.1.
POJEM „SÍŤ“ A JEHO VÝZNAM PRO MATEMATICKÉ MODELOVÁNÍ ..................................... 16
2.2.
GAMBIT, PRVKY SÍTĚ ....................................................................................................... 17
2.3.
KRITÉRIA PRO POSOUZENÍ KVALITY SÍTĚ ........................................................................ 19
2.4.
PŘÍKLAD VYTVOŘENÍ SÍTĚ ............................................................................................... 20
2.4.1.
Vytvoření sítě ............................................................................................................. 21
2.4.2.
Definování okrajových podmínek v Gambitu .......................................................... 23
3.
PROGRAMOVÝ SYSTÉM FLUENT .................................................................................. 24
3.1.
PŘEHLED METOD ŘEŠENÍ PARCIÁLNÍCH DIFERENCIÁLNÍCH ROVNIC................................ 24
3.2.
INTEGRACE METODOU KONEČNÝCH OBJEMŮ.................................................................. 24
3.3.
VÝBĚR INTERPOLAČNÍHO SCHEMATU.............................................................................. 27
3.4.
KONVERGENCE. .............................................................................................................. 28
3.4.1.
Residuály ................................................................................................................... 28
3.4.2.
Urychlení konvergence ............................................................................................. 30
3.4.3.
Relaxace .................................................................................................................... 30
4.
TURBULENTNÍ PROUDĚNÍ SKUTEČNÝCH KAPALIN................................................. 32
4.1.
KLASIFIKACE PROUDĚNÍ SKUTEČNÝCH KAPALIN ............................................................. 32
4.2.
TURBULENTNÍ PROUDĚNÍ ................................................................................................ 34
4.2.1.
Teorie turbulence ...................................................................................................... 35
4.2.2.
Metody matematického modelování turbulentního proudění ................................ 38
ii
Obsah
4.2.3.
Reynoldsova rovnice................................................................................................. 39
4.2.4.
Boussinesquova hypotéza o vírové (turbulentní) viskozitě.................................... 43
5.
STATISTICKÉ MODELY TURBULENCE ......................................................................... 45
5.1.
MODEL SMĚŠOVACÍ DÉLKY (NULAROVNICOVÝ MODEL) - PRANDTL ................................. 46
5.2.
JEDNOROVNICOVÝ MODEL .............................................................................................. 46
5.3.
DVOUROVNICOVÝ K-e MODEL ......................................................................................... 47
5.4.
RNG K-e MODEL ............................................................................................................. 48
5.5.
REYNOLDSŮV NAPĚŤOVÝ MODEL (RSM) ........................................................................ 48
5.6.
MODELOVÁNÍ PROUDĚNÍ V BLÍZKOSTI STĚNY, STĚNOVÉ FUNKCE .................................. 49
5.6.1.
Teorie stěnových funkcí dle Laundera a Spaldinga ............................................... 50
5.6.2.
Modelování proudění v blízkosti stěny ve Fluentu ................................................. 51
5.6.3.
Nerovnovážná stěnová funkce ................................................................................. 54
5.6.4.
Použití stěnových funkcí a omezení ........................................................................ 55
5.6.5.
Dvouvrstvový model (near-wall modelling) ............................................................. 55
5.6.6.
Vliv kvality sítě na volbu stěnové funkce pro různé modely turbulence ............... 56
5.6.7.
Vliv drsnosti na stěnovou funkci............................................................................... 57
6.
MATEMATICKÝ MODEL TURBULENCE PRO STLAČITELNÉ NEIZOTERMNÍ
PROUDĚNÍ ..................................................................................................................................... 58
6.1.
K-e DVOUROVNICOVÝ MODEL TURBULENCE .................................................................... 58
6.2.
OKRAJOVÉ PODMÍNKY PRO K-e TURBULENTNÍ MODEL .................................................... 59
6.2.1.
Hmotnostní průtok ..................................................................................................... 59
6.2.2.
Turbulentní veličiny ................................................................................................... 59
6.2.3.
Tlak na vstupu ........................................................................................................... 60
6.2.4.
Tlak na výstupu ......................................................................................................... 62
6.2.5.
Outflow ....................................................................................................................... 62
7.
ŘEŠENÍ PŘENOSU TEPLA (KONVEKCE, KONDUKCE) ............................................. 63
7.1.
ÚVOD DO PROBLEMATIKY MODELOVÁNÍ PŘENOSU TEPLA............................................... 63
7.1.1.
Kondukce – vedení tepla .......................................................................................... 63
7.1.2.
Konvekce ................................................................................................................... 65
7.2.
MATEMATICKÝ MODEL PŘENOSU TEPLA ......................................................................... 66
7.2.1.
Přenos tepla v tekutinách ......................................................................................... 66
7.2.2.
Přenos tepla ve vodivých stěnách ........................................................................... 67
7.3.
HUSTOTA ZÁVISLÁ NA TEPLOTĚ A TLAKU, BOUSSINESQOVA APROXIMACE ..................... 68
7.3.1.
Vyjádření hustoty pro stlačitelné médium ............................................................... 68
7.3.2.
Vyjádření hustoty pro nestlačitelné médium ........................................................... 68
7.3.3.
Boussinesquova aproximace ................................................................................... 69
iii
Obsah
7.4.
OKRAJOVÉ PODMÍNKY PRO NEIZOTERMNÍ PROUDĚNÍ ..................................................... 71
7.4.1.
Okrajové podmínky na tenké stěně ......................................................................... 73
7.4.2.
Okrajové podmínky na tenké dvoustranné stěně ................................................... 74
7.4.3.
Výpočet teploty a hustoty tepelného toku na stěně................................................ 75
8.
MODELOVÁNÍ PROUDĚNÍ PŘÍMĚSÍ ............................................................................... 80
8.1.
TRANSPORTNÍ ROVNICE PRO PŘENOS PŘÍMĚSÍ .............................................................. 80
8.2.
DEFINICE ZDROJE PŘÍMĚSI ............................................................................................. 81
8.2.1.
Zdroj-inlet ................................................................................................................... 82
8.2.2.
Objemový zdroj.......................................................................................................... 82
8.3.
FYZIKÁLNÍ VLASTNOSTI PLYNŮ A JEJICH SMĚSÍ ............................................................... 83
8.3.1.
Hustota ....................................................................................................................... 83
8.3.2.
Viskozita ..................................................................................................................... 84
8.3.3.
Měrná tepelná kapacita ............................................................................................ 85
8.3.4.
Tepelná vodivost ....................................................................................................... 85
8.3.5.
Standardní slučovací entalpie a entropie ................................................................ 85
8.4.
9.
OKRAJOVÉ PODMÍNKY PRO PŘÍMĚSI NA VSTUPU, VÝSTUPU A STĚNĚ.............................. 86
VÍCEFÁZOVÉ MODELY ..................................................................................................... 89
9.1.
VÍCEFÁZOVÉ MODELY OBECNĚ ....................................................................................... 89
9.1.1.
Euler-Lagrangeův přístup ......................................................................................... 90
9.1.2.
Euler-Eulerův přístup ................................................................................................ 90
9.2.
10.
VÍCEFÁZOVÉ MATEMATICKÉ MODELY .............................................................................. 91
9.2.1.
Vícefázový model směsi (mixture model) ............................................................... 93
9.2.2.
Trajektorie pevných částic v plynu ........................................................................... 97
ČASOVĚ ZÁVISLÉ ŘEŠENÍ ............................................................................................ 103
10.1.
DISKRETIZACE ČASOVĚ ZÁVISLÉ ROVNICE.................................................................... 104
10.2.
OKRAJOVÉ PODMÍNKY................................................................................................... 106
10.2.1.
Tabulka pro časovou okrajovou podmínku........................................................ 106
10.2.2.
UDF pro okrajovou podmínku............................................................................. 108
10.3.
11.
PŘÍKLAD VYHODNOCENÍ ČASOVĚ ZÁVISLÉ ÚLOHY ........................................................ 108
MĚŘENÍ A VÝPOČET TURBULENCE PŘI OBTÉKÁNÍ VÁLCE
V AERODYNAMICKÉM TUNELU ............................................................................................. 113
11.1.
TEORETICKÝ ROZBOR ÚLOHY OBTÉKÁNÍ VÁLCE............................................................ 113
11.2.
FYZIKÁLNÍ EXPERIMENT ................................................................................................ 113
11.3.
MATEMATICKÝ MODEL................................................................................................... 118
11.4.
VYHODNOCENÍ VÝSLEDKŮ FYZIKÁLNÍHO EXPERIMENTU A MATEMATICKÝCH MODELŮ .. 121
iv
Obsah
12.
LITERATURA .................................................................................................................... 128
13.
PŘÍLOHA - SČÍTACÍ PRAVIDLA .................................................................................... 131
13.1.
DEFINICE A PRAVIDLA ................................................................................................... 131
13.2.
PŘÍKLADY...................................................................................................................... 132
13.3.
SROVNÁNÍ S VEKTOROVÝM OZNAČENÍM ....................................................................... 133
13.4.
SUBSTANCIÁLNÍ DERIVACE ........................................................................................... 134
14.
PŘÍLOHA - ZÁKLADNÍ STATISTICKÉ METODY ......................................................... 136
14.1.
STŘEDNÍ HODNOTA ....................................................................................................... 136
14.2.
REYNOLDSOVA PRAVIDLA ............................................................................................. 137
14.3.
ROZPTYL, SMĚRODATNÁ ODCHYLKA, TURBULENTNÍ INTENZITA, KOVARIANCE, KORELACE
138
15.
PŘÍLOHA - ROZKLAD V TAYLOROVU ŘADU............................................................. 140
15.1.
DEFINOVÁNÍ PROBLÉMU................................................................................................ 140
15.2.
TAYLORŮV ROZVOJ ....................................................................................................... 140
v
Seznam použitých označení
SEZNAM POUŽITÝCH OZNAČENÍ
Poznámka: označení, u něhož není uveden rozměr, reprezentuje obecnou proměnnou.
a
okamžitá hodnota veličiny
a¢
fluktuační složka veličiny
a
r
a
časově středovaná složka veličiny
A, S
plocha
[m2 ]
A
Van Driestova konstanta
[1]
CD
konstanta
[1]
Cn
empirická konstanta
[1]
Cm
konstanta
[1]
C1e
empirická konstanta
[1]
C 2e
empirická konstanta
[1]
C 3e
empirická konstanta
[1]
c
rychlost zvuku
[ms-1]
cv
měrná tepelná kapacita při konstantním objemu
[J×kg-1×K-1]
cp
měrná tepelná kapacita při konstantním tlaku
[J×kg-1×K-1]
Dn,m
difúzní koeficient pro příměs n ve směsi
[m2s-1]
f
frekvence
[s-1]
f i , Fi
síla
[N]
fc
Coriolisův parametr
[s ]
E
energie
[J×kg-1]
E
empirická konstanta
[1]
Fi ( J )
i-tá složka složka vektoru hustoty toku veličiny J
G
termická produkce turbulentní kinetické energie
[m2s-3]
Gr
Grashofovo číslo
[1]
g
tíhové zrychlení
[ms-2]
hf
součinitel přestupu tepla
[W×m-2×K-1]
obecný vektor
-1
vi
Seznam použitých označení
h
statická entalpie
[J×kg-1]
ho
celková entalpie
[J×kg-1]
I
r r r
i , j, k
intenzita turbulence
[%]
jednotkový vektor ve směru os x, y, z
[1]
J
bilancovaná veličina
J e, J p , J w
hmotnostní tok přes stěny konečných objemů
[kg×m-2s-1]
k
turbulentní kinetická energie
[m2s-2]
kP
turbulentní kinetická energie v logaritmické vrstvě
[m2s2]
I, L
délkové měřítko turbulence
[m]
lm
směšovací délka
[m]
M
Machovo číslo
[1]
nj
j-tá složka normálového vektoru
[1]
r
n
vektor vnější normály k ploše
[1]
p
tlak
[Pa]
pop
operační tlak
[Pa]
ps
statický tlak
[Pa]
P
mechanická produkce turbulentní kinetické energie
[m2s-3]
P( J )
hustota produkce veličiny J
Pr
molekulové Prandtlovo číslo (stěnové funkce)
[1]
Prt , s h
turbulentní Prandtlovo číslo
[1]
Q
průtok
[ m3s-1]
r
měrná plynová konstanta
[J.kg-1.K-1]
q ¢¢
tepelný tok
[J×m -2×s-1]
R
reziduál
R
normalizovaný reziduál
[1]
Ra
Rayleighovo číslo
[1]
Re
Reynoldsovo číslo
[1]
Re y
turbulentní Reynoldsovo číslo buňky
[1]
rw
polohový vektor na stěně
[m]
S
modul tensoru střední rychlosti deformace
[s-1]
vii
Seznam použitých označení
S i, j
tensor rychlosti deformace
[s-1]
Sc
Schmidtovo číslo
[1]
t
čas
[s]
T
absolutní teplota
[K]
ui
i-tá složka rychlosti
[ms-1]
ui
i-tá složka střední rychlosti
[ms-1]
u i¢
i-tá složka fluktuační rychlosti
[ms-1]
u + , u*
rychlost definovaná stěnovou funkcí
[ms-1]
ut
třecí rychlost
[ms-1]
V
objem
[m3]
xi
souřadnice v kartézském systému [x1, x2, x3] nebo [x, y, z]
y
kolmá vzdálenost od stěny
[m]
y + , y*
bezrozměrná veličina při odvozování stěnových funkcí
[1]
bezrozměrná tloušťka podvrstvy
[1]
yv
tloušťka vazké podvrstvy
[m]
Dy P
vzdálenost bodu P od stěny ve směru normály
[m]
a
relaxační faktor
[1]
a
teplotní vodivost
[m2s-1]
a k ,a e ,a z
inverzní turbulentní Prandtlova čísla
[m2s-1]
b
součinitel teplotní roztažnosti
[K-1]
d ij
Kroneckerovo delta
[1]
e
rychlost disipace
[m2s-3]
eP
rychlost disipace v logaritmické vrstvě
[m2s-3]
F
disipační funkce
Gt
turbulentní difusivita
[m2s-1]
k
von Kármánova konstanta, poměr měrných tepelných kapacit
[1]
l
součinitel tepelné vodivosti
[W×m-1×K-1]
lt
součinitel efektivní tepelné vodivosti
[W×m-1×K-1]
m
dynamická viskozita
[Pa×s]
yv
*
viii
Seznam použitých označení
m eff
efektivní viskozita
[Pa×s]
mt
turbulentní viskozita
[Pa×s]
n
kinematická viskozita
[m2s-1]
nt
turbulentní viskozita
[m2s-1]
p ij
celkový tenzor napětí
[Pa]
r
hustota
[kg×m-3]
r ref
referenční hustota
[kg×m-3]
sk
empirická konstanta
[1]
se
empirická konstanta
[1]
sh
turbulentní Prandtlovo číslo
[1]
t
časová perioda
[s]
t
vazké napětí
[Pa]
t ij
tenzor vazkých napětí
[Pa]
tw
vazké napětí na stěně
[Pa]
tt
turbulentní napětí
[Pa]
x
druhá viskozita
[Pa×s]
z
obecná proměnná
Indexy:
i
index složky rychlosti, index iterace
[1]
i
sčítací index
[1]
j, k , l
sumační Einsteinův index
[1]
w , e, p
index stěny konečného objemu
[1]
W , E, P , N , S , F , B, NB
[1]
index konečného objemu
ix
Předmluva
Předmluva
Tato skripta jsou určena pro posluchače všech fakult magisterských a doktorských
studijních proramů, kteří se chtějí seznámit
se základy numerického
modelování
přenosových jevů v tekutině, tj. přenosu hmoty, hybnosti (momentu), tepla, příměsí atd. při
laminárním a turbulentním proudění. Navazují na předchozí vydání s názvem „Matematické
modely nestlačitelného a stlačitelného proudění. Metoda konečných objemů.“ autorů M.
Kozubková, S. Drábková, P. Šťáva, které bylo přepracováno dle zkušeností při aplikacích ve
výuce a při užití software Fluent a CFX.
Numerické modelování mnoha fyzikálních jevů je úzce spojeno s modelováním určité
formy pohybu matematickými prostředky. Pohyb tekutin souvisí s řešením nejrůznějších
problémů, daných fyzikálním modelem:
·
laminární a turbulentní proudění v jednoduchých i složitých geometriích
·
stlačitelné a nestlačitelné proudění
·
stacionární, nestacionární a přechodové proudění
·
přenos tepla, přirozená a smíšená konvekce, radiace
·
přenos chemické příměsi včetně chemických reakcí
·
vícefázové proudění, proudění s volnou hladinou, proudění s pevnými částicemi,
bublinami, resp. kapkami
·
hoření a chemické reakce
·
proudění porézním prostředím, atd.
Matematický model spočívá v definici rovnic výše uvedené děje popisujících. Vzhledem
k tomu, že se jedná o děje rovinné dvourozměrné, osově symetrické nebo obecně
trojrozměrné a časově závislé, jsou popsány soustavou parciálních diferenciálních rovnic,
kterou je nutné řešit numerickými metodami. Jejich využívání je podmíněno rozšířením
znalostí z oblasti proudění, turbulence, numerických metod, výpočetní techniky.
K řešení proudění je možno využít komerční programové systémy, jako je Fluent,
CFX a další. Úkolem uživatele je sestavní správného výpočtového modelu, což obsahuje
některé matematické, fyzikální a technické principy. Pro takový model je nutné najít všechny
vstupní údaje v platných normách, sestavit vstupní data pro program, kterým lze výpočtový
model řešit, provést řešení u terminálu, správně interpretovat výsledky pro další použití a ve
všech fázích provádět účinné kontroly všech vstupů a výstupů. Uživatel musí bezpečně
rozčlenit všechny informace na údaje geometrické (dvourozměrné nebo třírozměrné útvary,
topologie), údaje o působení vnějších sil a fyzikální údaje (informace o proudícím médiu,
x
Předmluva
jeho
fyzikálních
vlastnostech).
Tedy
nezastupitelnou
úlohou
uživatele
je
znalost
hydromechaniky, termomechaniky a dalších věd podle složitosti problému.
Pokud jde o výpočetní metody, na nichž jsou založeny užívané programy, měl by
projektant znát jejich podstatu v rozsahu potřebném pro spolehlivé použití ve standardních
případech. U programu Fluent resp. CFX je třeba vědět, s jakými tvary konečných objemů se
bude pracovat, z toho vyplývá volba hustoty sítě, jaká aproximační schémata bude vhodné
použít, u dynamiky mít představu o charakteru časové závislosti jednotlivých veličin a z toho
vyplývající velikosti časového kroku, apod.
Neméně významnou částí je vyhodnocení výsledků, které je obzvlášť obtížné u
trojrozměrných úloh. Je optimální mít k dispozici alespoň orientační hodnoty počítaných
veličin, ideální je srovnání výsledků s experimentem.
Skripta jsou členěna do tří základních oblastí. V první oblasti je velmi stručně
definován matematický model nejjednoduššího typu proudění, popsána tvorba sítě a příklady
ilustrující zjednodušeně numerické řešení. Druhá část je věnována modelům turbulence,
stěnovým funkcím a okrajovým podmínkám. Pak následují aplikace na řešení proudění
s příměsí, přenosem tepla a vícefázové proudění. Tyto kapitoly jsou doplněny příklady
v jednoduchých geometriích, aby je bylo možno snadno testovat bez vytváření nové
geometrie a věnovat se novým matematickým modelům a způsobům vyhodnocení, které
v trojrozměrné geometrii není jednoduché. V přílohách jsou vysvětleny použité matematické
pojmy a označení.
V označení veličin se mohou vyskytnou jisté nejednotnosti, které jsou způsobené
čerpáním podkladů z literatury české a zahraniční. Úplné sjednocení by jistě bylo možné, ale
vzhledem k tomu, že tato skripta jsou jen určitým vodítkem pro numerické modelování a jistě
bude nutné doplňovat znalosti z doporučené literatury a především manuálu Fluentu, bylo
někde ponecháno označení veličin v souladu s tímto manuálem.
Skripta jsou doplněna dvěma skripty určenými pro cvičení s názvy:
Bojko, M.: Návody do cvičení „Modelování proudění“ Fluent
Blejchař, T.: Návody do cvičení „Modelování proudění“ CFX.
xi
Úvod do numerického modelování proudění
1. Úvod do numerického modelování proudění
1.1.
Úvod do proudění tekutin
Proudění kapalin je možno rozdělit podle několika hledisek [25]:
1.1.1. Dělení podle fyzikálních vlastností tekutiny
PROUDĚNÍ TEKUTINY
proudění ideální
(nevazké) tekutiny
potenciální
proudění skutečné
(vazké) tekutiny
laminární
vířivé
turbulentní
Proudění ideální (dokonalé) tekutiny
·
Potenciální proudění (nevířivé)
Částice tekutiny se pohybují přímočaře nebo křivočaře po dráhách tak, že vůči pozorovateli
se neotáčejí kolem vlastní osy viz obr. 1.1. Natočení částice na křivé dráze je kompenzováno
stejně velkým natočením částice kolem vlastní osy, ale v opačném smyslu. Mezi potenciální
proudění patří rovněž potenciální vír, u něhož částice krouží kolem vírového vlákna
potenciálně s výjimkou částice, která tvoří vlákno viz obr. 1.2.
·
Vířivé proudění
Částice tekutiny se vůči pozorovateli natáčejí kolem vlastních os viz obr. 1.3
obr. 1.1 Potenciální proudění
obr. 1.2 Potenciální vír
obr. 1.3 Vířivé proudění
Proudění skutečné (vazké) tekutiny
·
Laminární proudění
Částice tekutiny se pohybují v tenkých vrstvách, aniž se přemísťují po průřezu, viz obr. 1.4.
·
Turbulentní proudění
Částice tekutiny mají kromě podélné rychlosti také turbulentní (fluktuační) rychlost, jíž se
1
Úvod do numerického modelování proudění
přemísťují po průřezu viz obr. 1.5.
obr. 1.4 Laminární proudění
obr. 1.5 Turbulentní
proudění
1.1.2. Dělení podle kinematických hledisek
PROUDĚNÍ TEKUTINY
uspořádání v prostoru
1D
2D
závislost na čase
3D
stacionární
nestacionární
Dělení proudění dle uspořádání v prostoru
·
Proudění je třírozměrné neboli prostorové (3D), jestliže veličiny, např. rychlost, závisí
na poloze v prostoru v = v (x, y , z )
·
Proudění dvourozměrné neboli rovinné (2D) je charakterizované veličinami, jako je
např. rychlost, závisí na poloze v rovině (příkladem je osově symetrické proudění
v potrubí) v = v (x, y )
·
Proudění jednorozměrné (1D) předpokládá závislost počítaných veličin na poloze na
křivce (příkladem je proudění v potrubních systémech) v = v (s )
Dělení proudění podle závislosti na čase
¶
=0
¶t
·
Proudění ustálené (stacionární) nezávisí na čase v ¹ v (t ) ;
·
Proudění neustálené (nestacionární) je proudění, u něhož veličiny jsou závislé na
čase, v = v (x, y , z, t ) ; v = v (s, t ) ; v = v (t ) .
2
Úvod do numerického modelování proudění
V nejobecnějším případě je poloha
y (x2)
bodu
X=
definována
(x, y , z )
Vektor
v (u2)
w (u3)
složkami
z (x3)
u=
souřadnicemi
resp. X = (x1, x 2 , x3 ) .
rychlosti
u=
(u1, u2, u3 ) .
je
(u,v ,w )
definován
resp.
Označení je patrné
z obr. 1.6
u (u1)
x (x1)
obr. 1.6 Souřadný systém
1.2.
Přenos hmoty, hybnosti, tepla při neizotermním proudění
nestlačitelné tekutiny
Základní fyzikální zákony popisující proudění jsou zákony zachování hmotnosti,
hybnosti, tepla případně dalších skalárních veličin. Jsou vyjádřeny Navierovými Stokesovými
rovnicemi spolu s rovnicí kontinuity a popisují laminární i turbulentní režim proudění.
V případě nestacionárního nestlačitelného neizotermního proudění mají následující tvar
([20]):
Rovnice kontinuity:
¶u ¶v ¶w
+
+
=O
¶x ¶y ¶z
( 1.2.1)
Navier-Stokesovy rovnice:
æ ¶ 2 u ¶ 2u ¶ 2u ö
1¶p
¶ u ¶ (uu ) ¶ (uv ) ¶ (uw )
÷+f
+
+
+
=+n ç
+
+
ç ¶ x 2 ¶ y 2 ¶ z2 ÷ x
¶t
¶x
¶y
¶z
r ¶x
ø
è
æ ¶ 2v ¶ 2v ¶ 2v ö
1 ¶p
¶ v ¶ (vu ) ¶ (vv ) ¶ (vw )
÷+f
+
+
+
=+n ç
+
+
ç ¶ x 2 ¶ y 2 ¶ z2 ÷ y
¶t
¶x
¶y
¶z
r ¶y
ø
è
( 1.2.2)
æ ¶ 2w ¶ 2w ¶ 2w ö
1¶p
¶ w ¶ (wu ) ¶ (wv ) ¶ (ww )
÷+f
+
+
+
=+n ç
+
+
ç ¶ x2 ¶ y 2 ¶ z2 ÷ z
¶t
¶x
¶y
¶z
r ¶z
ø
è
kde podle schématu na obr. 1.6 jsou u , v a w složky rychlosti, p tlak, r hustota, n
kinematická viskozita a f x ,y ,z označuje složky vnější objemové síly (gravitační, odstředivé
sily).
Rovnice pro přenos tepla, tj. zákon zachování energie je ve tvaru
3
Úvod do numerického modelování proudění
æ ¶ 2T ¶ 2T ¶ 2T
¶ T ¶ (uT ) ¶ (vT ) ¶ (wT )
+
+
+
= a çç 2 +
+
¶t
¶x
¶y
¶z
¶ y 2 ¶ z2
è¶ x
ö
÷÷ + af
ø
æ æ ¶u ö 2 æ ¶v ö 2 æ ¶w ö 2 ö
÷ ÷+
÷ +ç
÷ +ç
f = 2ç çç
ç è ¶ x ÷ø çè ¶ y ÷ø çè ¶ z ÷ø ÷
ø
è
( 1.2.3)
æ æ ¶u ¶v ö 2 æ ¶u ¶w ö 2 æ ¶v ¶w ö 2 ö
÷ ÷
÷ +ç
÷ +ç
+
+
+ ç çç
+
ç è ¶ y ¶ x ÷ø çè ¶ z ¶ x ÷ø çè ¶ z ¶ y ÷ø ÷
è
ø
kde a =
l
je teplotní vodivost, l je molekulová tepelná vodivost a c p je měrné teplo.
rc p
Při vyjádření proměnných o třech případně devíti složkách (složky rychlostí, napětí
apod.) je vhodné využít speciální zkrácené označení s přesně definovanými pravidly, známé
jako Einsteinova sumace, viz kap. 13, kdy pouze jedním členem lze vyjádřit všechny tři
složky rychlostí resp. devět napětí. Totéž lze pro přehlednost vyjádřit matematicky užitím
znaku sumy. Tedy rovnice kontinuity se zapíše zjednodušeně:
¶ u1 ¶ u2 ¶ u3
+
+
= O resp.
¶ x1 ¶ x2 ¶ x3
n
¶uj
å¶x
j =1
= O resp.
j
¶ uj
=O
¶ xj
( 1.2.4)
Navier-Stokesovy rovnice podobně:
n
1 ¶p
¶ ui n ¶ (ui u j )
¶ 2ui
+å
=+n å
+ fi
2
¶ t j =1 ¶ x j
r ¶ xi
j =1 ¶ x j
resp.
1 ¶p
¶ ui ¶ (ui u j )
¶ 2ui
+
=+n
+ fi , i = 1,..., n
¶t
¶ xj
r ¶ xi
¶ x 2j
( 1.2.5)
kde důsledně index i vyjadřuje složku vektoru a index j (případně další podle abecedy)
vyjadřuje sčítací index ( j = 1,2 resp. 3 ).
Rovnice pro přenos tepla lze podobně zapsat:
n
n
¶ (u jT )
¶T
¶ 2T
+å
= aå
+ af ,
2
¶ t j =n ¶ x j
j =n ¶ x j
1 n n æ ¶u
¶u
f = åå çç j + l
2 j =1 l =1 è ¶ xl ¶ x j
ö
÷
÷
ø
( 1.2.6)
2
resp.
¶ T ¶ (u jT )
¶ 2T
+
=a
+ af
¶t
¶ xj
¶ x 2j
4
Úvod do numerického modelování proudění
¶u ö
1 æ ¶u
f = çç j + l ÷÷
2 è ¶ xl ¶ x j ø
1.3.
2
Parciální diferenciální rovnice a okrajové podmínky
Rovnice pro zachování hmotnosti, hybnosti, rovnice energie a rovnice pro transport
chemické příměsi v obecné konzervativní formě tvoří systém parciálních diferenciálních
rovnic, přitom všechny rovnice lze formálně vyjádřit zápisem
¶ ( rz )
¶
=(ru j z ) + ¶
¶t
¶ xj
¶ xj
akumulace
konvekce
é
¶z ù
êa z
ú + Sz
êë ¶ x j úû
difúze
( 1.3.1)
zdroj
kde z je proměnná a členy na pravé straně jsou postupně konvektivní, difúzní a zdrojový
člen, proto se rovnice nazývá také konvekčně-difúzní rovnice. Pokud z představuje teplotu,
příměs nebo jinou skalární veličinu, pak se jedná o lineární rovnici druhého řádu, pokud z
představuje složku rychlosti, lze tuto rovnici lze považovat za nelineární rovnici druhého
řádu.
Převládá-li vliv difúzního členu, jedná se o rovnice eliptické, u parabolických rovnic
převládá vliv konvektivního transportu a u hyperbolických rovnic jsou významné tlakové
změny. Rovnice ( 1.3.1) však sama nestačí k určení funkce z . Aby byla funkce určena, je
třeba znát ještě doplňující podmínky.
Nechť je definován problém najít rozložení teploty v daném tělese. Předpokládejme,
že hledaná funkce T (tedy teplota) je funkcí následných prostorových nezávisle proměnných
a čase
T=T(t,x)
¶T
¶ 2T
=a 2
¶t
¶x
( 1.3.2)
T=T(t,x, y)
¶T
¶ 2T
¶ 2T
=a
+
a
¶t
¶x 2
¶y 2
T=T(t, x, y, z)
¶T
¶ 2T
¶ 2T
¶ 2T
=a 2 +a
+a 2
¶t
¶x
¶y 2
¶z
Pokud je funkce T závislá pouze na prostorových souřadnicích a nezávislá na čase,
úloha se nazývá stacionární. Problém se zkomplikuje, jestliže je funkce T závislá na čase a
5
Úvod do numerického modelování proudění
eventuelně na prostorových souřadnicích, jak bylo uvedeno výše. Např. těleso je ohříváno
nebo ochlazováno s rostoucím časem, pak se úloha nazývá nestacionární a funkce T je
funkcí následných nezávisle proměnných.
NEZÁVISLE PROMĚNNÉ
stacionární proudění
nestacionární proudění
t
x
t, x
x, y
t, x, y
x, y, z
t, x, y, z
Fyzikální děje (tj. rozložení teplot), jsou popisovány diferenciálnímí rovnicemi dvojího typu
·
obyčejné diferenciální rovnice,
resp.
tj. T = T (t ) , resp. T = T (x ) , např.
dT
= g (x ) , tj. T je funkcí pouze jedné nezávislé proměnné
dx
parciální diferenciální rovnice, tj. T = T ( x, y ) ,
·
dT
= f (t ) ,
dt
T = T ( x, y , z ) , T = T ( t , x ) ,
T = T ( t , x, y ) , T = T ( t , x, y , z )
Např. rovnice vedení tepla v tyči je
t
¶T
¶ 2T
=a
¶t
¶x 2
D
T=T1
Pro rovnici vedení tepla v tyči se řešení
T=T2
hledá v obdélníku D(x, t ) , viz obr. 1.7, a
musí splňovat tzv. počáteční podmínku
T(x)=j(x)
x
L
(jistá
analogie
s řešením
obyčejných
diferenciálních rovnic)
obr. 1.7 Definice oblasti
Pro teplotu T je počáteční podmínka definována následovně:
T (x,0 ) = j (x )
0á xáL
( 1.3.3)
V případě, že se v rovnici vyskytuje i druhá derivace funkce T podle t, bude definovaná i
druhá počáteční podmínka jako
¶T (x,0 )
= y (x ) .
¶t
Okrajová podmínka definuje teplotu na počátku a konci oblasti (okraj tyče)
T (0, t ) = T1 (t )
( 1.3.4)
T (L, t ) = T2 (t )
6
Úvod do numerického modelování proudění
Úloha najít řešení rovnice ( 1.3.1) splňující okrajové i počáteční podmínky se nazývá
smíšenou úlohou. Jsou-li okrajové podmínky rovny nule, nazývají se homogenní okrajové
podmínky, podobně jsou-li počáteční podmínky rovny nule, nazývají se homogenní
počáteční podmínky. Místo okrajových podmínek ( 1.3.4) mohou být dány podmínky jiného
typu, které se též nazývají okrajové.
Úvaha o okrajových a počátečních podmínkách pro teplotu je platná pro obecnou
proměnnou z .
1.4.
Okrajové podmínky proudění nestlačitelné tekutiny
1.4.1. Typy okrajových podmínek
Okrajové podmínky nemusí být jen konstantní veličiny, ale mohou nabývat hodnot
definovaných funkcí, tabulkou atd. ([15]):
·
konstanty y = konst.
·
polynom. funkce y ( x ) = A0 + A1x + A2x 2 + ... ,
0.05
0.045
kde koeficienty se zadávají pouze na 5 platných
0.04
0.035
·
y [m]
cifer
derivace podle normály (OUTLET, teplotní
¶y (x )
tok)
= konst1.
¶x
·
u-polynom
0.025
u-po částech
lin. funkce
0.02
0.015
0.01
0.005
0
po částech lineární funkce (piecewice linear)
(x1, y1 ), (x2, y 2 ), (x3, y 3 ),
·
u-konst
0.03
...
0
(xN , y N )
0.5
1
1.5
2
u [m.s-1]
obr. 1.8 Profily rychlostí
kombinace polynom. a po částech lin.
funkce
1.4.2. Podmínky vstupu a výstupu.
Pro dvě průtočné hranice mohou
nastat
pouze
kombinace
následující
okrajových
o u tlet
základní
podmínek,
rych lo st
tla k stat.
rych lo st
tla k stat.
(kombinace vstupní rychlosti a výstupní
rychlosti nemůže nastat, protože rychlost
na druhém vstupu se počítá z rovnice
tlak totál.
spojitosti). Při uvažování rovnice energie
se zadává navíc hodnota teploty, viz.obr.
obr. 1.9 Kombinace vstupních a výstupních
okrajových podmínek
1.9.
Podmínky na průtočných hranicích - na vstupních a výstupních průtočných hranicích
7
Úvod do numerického modelování proudění
lze definovat tři typy okrajových podmínek
VSTUP
VÝSTUP
rychlost u
podmínka ustáleného proudu
¶u
¶p
¶T
= 0,
= 0,
=0
¶n
¶n
¶n
totální (celkový) tlak ptot
ptot = pstat + pdyn = pstat +
statický tlak pstat
1 2
ru
2
teplota T
Podmínky na stěně - stěna může být nepohyblivá (rychlost tekutiny je rovna nule) nebo
pohyblivá (např. rotující nebo klouzající), se třením nebo bez tření, hladká nebo drsná.
Teplota je dána formou hodnoty nebo derivace podle nomály ke stěně, např.
¶T
= 0 je
¶n
izolovaná stěna.
Podmínky symetrie - nulová normálová rychlost a nulové normálové gradienty všech
hledaných veličin, viz obr. 1.11.
Podmínky osové symetrie – definují osu při osově symetrických dvourozměrných
úlohách, viz obr. 1.10.
řešená
oblast
řešená
oblast
osa
symetrie
rovina
sym etrie
obr. 1.10 Válec s osou symetrie a rovinou, pro
obr. 1.11 Válec s rovinou symetrie a
kterou je řešeno proudění.
oblast, pro kterou je řešeno proudění.
Periodické (cyklické) podmínky – používají se v případě, kdy se opakují proudové
útvary, mohou být rotačního typu a translačního typu, kdy se umožňuje definování tlakového
spádu ve směru proudící tekutiny po celé délce oblasti.
Všechny typy podmínek mohou být časově závislé, pokud to vyžaduje jejich charakter.
Další okrajové podmínky se netýkají proudění jako takového, ale dalších veličin vyplývajících
ze složitosti matematického modelu, jako je skalární veličina teplota, teplotní toky, radiace,
hmotnostní zlomky (resp. molové zlomky) příměsí apod.
8
Úvod do numerického modelování proudění
periodické
periodické
podmínky
podmínky
řešená
řešená
oblast
oblast
periodické podmínky
periodické podmínky
rotačního typu
translačního typu
obr. 1.12 Periodické podmínky rotačního a translačního typu
1.5.
Řešení rovnice vedení tepla
Cílem numerických metod pro řešení parciálních diferenciálních rovnic je hledat
diskrétní řešení definované v dostatečně malých podoblastech základní oblasti pomocí
systému tzv. diferenčních (algebraických) rovnic v základních bodech
·
dělení oblasti na diskrétní geometrické elementy – vytvoření sítě
·
bilancování neznámých veličin v konečných objemech nebo uzlech a diskretizace
·
numerické řešení diskretizovaných rovnic v obecném tvaru
přitom diskretizační chyba se definuje jako rozdíl mezi řešením diferenciálních a diferenčních
rovnic. Základní vlastnosti numerických metod jsou:
·
míra přesnosti diskretizační chyby a residuálu
·
míra stability
Existuje určitý vývoj v numerickém řešení rovnic definujících proudění tekutin.
Nejstarší klasickou metodou je diferenční metoda. Princip diferenční metody pro řešení
diferenciálních rovnic lze popsat následovně
·
oblast, ve které se hledá řešení, se pokryje sítí složenou z konečného počtu
nepřekrývajících se elementů. Nejjednodušší sítě jsou
9
Úvod do numerického modelování proudění
Dx
úsečky v jednorozměrném případě
Dx
obdélníky ve dvourozměrném případě
Dy
Dy
šestistěny ve trojrozměrném případě
Dz
Dx
·
v těchto bodech se nahradí derivace diferencemi o různých přesnostech (např.
æ ¶T ö æ DT ö Ti +1 - Ti
), vztahy pro potřebné derivace se odvodí z Taylorova rozvoje
ç
÷ »ç
÷ =
Dx
è x ø i è Dx ø i
se specifickým označením souvisejícím s prouděním, viz kap. 15,
·
diferenciální rovnice přejde na soustavu algebraických rovnic o neznámých, které určují
přibližné hodnoty neznámé funkce ve všech uzlech sítě,
·
soustava algebraických rovnic se řeší numericky.
Příklad
Řešte rovnici vedení tepla v tyči, danou parabolickou diferenciální rovnicí
¶T
¶ 2T
=a
¶t
¶x 2
Řešení se hledá v obdélníku D (t, x ) , viz obr. 1.13, musí splňovat podmínky:
počáteční podmínka
okrajová podmínka
T (x,0 ) = T0 (x ) = 20 O C
T (0, t ) = T1(t ) = 80O C
T (L, t ) = T2 (t ) = 200 C
10
0á xáL .
Úvod do numerického modelování proudění
Diferenční rovnice vedení tepla má
D
i-1
t
i
i
Ti,n+1
tvar
i+1
Ti ,n +1 - Ti ,n
T
+ Ti -1,n - 2Ti ,n
= a i +1,n
Dt
Dx 2
i
Dx
n+1
Dt
a po úpravě platí
n
T=T1
n-1
Ti ,n +1 = Ti ,n + aDt
T=T2
Ti +1,n + Ti -1,n - 2Ti ,n
Dx 2
Tedy lze explicitně vyjádřit Ti ,n +1
pomocí
T=T0(x)
t=0
x
časovém
L
obr. 1.13 Geometrie oblasti, okrajové podmínky, síť
hodnot
v
kroku
jednoduchém
předchozím
n.
případě
V tomto
lze
najít
0.45
0.50
řešení v Excelu
Tabulka zadání parametrů pro iterační výpočet
a=
0.1
T(x=0)= 80
koef=
0.5
L=
1
T(x=L)= 20
Dx=
0.1
n=
10
T(t=0)= 20
Dt=
0.05
0.20
0.25
time
BC
x
BC
0.00
0.05
0.10
0.15
0.30
0.35
0.40
0
80.00 80.00
80.00 80.00 80.00
80.00 80.00 80.00 80.00 80.00 80.00
0.1
20.00 50.00
50.00 57.50 57.50
61.25 61.25 63.59 63.59 65.23 65.23
0.2
20.00 20.00
35.00 35.00 42.50
42.50 47.19 47.19 50.47 50.47 52.93
0.3
20.00 20.00
20.00 27.50 27.50
33.13 33.13 37.34 37.34 40.63 40.63
0.4
20.00 20.00
20.00 20.00 23.75
23.75 27.50 27.50 30.78 30.78 33.59
0.5
20.00 20.00
20.00 20.00 20.00
21.88 21.88 24.22 24.22 26.56 26.56
0.6
20.00 20.00
20.00 20.00 20.00
20.00 20.94 20.94 22.34 22.34 23.93
0.7
20.00 20.00
20.00 20.00 20.00
20.00 20.00 20.47 20.47 21.29 21.29
0.8
20.00 20.00
20.00 20.00 20.00
20.00 20.00 20.00 20.23 20.23 20.70
0.9
20.00 20.00
20.00 20.00 20.00
20.00 20.00 20.00 20.00 20.12 20.12
1
20.00 20.00
20.00 20.00 20.00
20.00 20.00 20.00 20.00 20.00 20.00
Konvergence úlohy závisí na volbě časového a prostorového kroku. Dalším problémem je
efektivní řešení této soustavy algebraických rovnic.
11
Úvod do numerického modelování proudění
80.00
70.00
60.00
50.00
T 40.00
30.00
10.00
0.50
0.60
0.70
0.80
0.90
1.00
20.00
0.30
0.20
0.10
1
0.00
0 0.1
0.2 0.3
0.4 0.5
0.6 0.7
length
0.8 0.9
0.40
time
0.00
obr. 1.14 Grafické zobrazení řešení v Excelu.
1.6.
Řešení proudění při náhlém rozšíření průřezu
Tato kapitola ilustruje, jak zadávat a řešit nestlačitelné proudění při náhlém rozšíření
průtočného průřezu ve Fluentu, což je jednou ze základních testovacích úloh proudění,
zkoumaných jak experimentálně tak numericky. Úkolem je:
· definovat fyzikální model, fyzikální vlastnosti proudícího média a stěn,
· definovat matematický model, okrajové podmínky,
· vytvořit geometrii a sítě,
· zadat okrajové a počáteční podmínky ve Fluentu, výpočet
· vyhodnotit vypočtené veličiny
Příklad
Fyzikální model je dán tvarem oblasti, typem proudění a hydraulickými parametry
proudění. Schéma oblasti je zobrazeno na obr. 1.15 a rozměry s fyzikálními vlastnostmi
v tabulce.
12
Úvod do numerického modelování proudění
Ls
ds
d
u
Oblast proudění
d-ds
z
L
obr. 1.15 Geometrie problému
Geometrie oblasti
délka oblasti L [m]
3.5
výška oblasti d [m]
0.5
rozšíření Ls xd s [mxm]
0.7 x 0.1
Fyzikální vlastnosti tekutiny
Předpokládá se proudění nestlačitelné tekutiny o parametrech
hustota tekutiny r [kg.m-3]
1.2
dynamická viskozita m resp. h [kg.(m.s)-1]
0.0000171
Okrajové podmínky
Okrajové podmínky na vstupu jsou definovány rychlostí a na výstupu je dána podmínka
ustáleného proudu (OUTFLOW). Na stěnách se předpokládá nulová rychlost proudění.
vstup - střední rychlost u s [m.s-1]
3
výstup
OUTFLOW
Reynoldsovo číslo Re =
ru s d
[1]
m
84705
Matematický model
Z okrajových a fyzikálních podmínek je možno určit Reynoldsovo číslo jako
bezrozměrné kritérium pro rozhodnutí, zda je proudění laminární nebo turbulentní. Hodnota
Reynoldsova čísla charakterizuje proudění v přechodové oblasti mezi laminárním a
turbulentním prouděním. Bude se tedy řešit na úvod jako laminární.
13
Úvod do numerického modelování proudění
Vytvoření geometrie a sítě
V prostředí GAMBIT se vytvoří přesná geometrie metodou podobnou prostředí CAD
programů. Navíc se využije možností tohoto programu tvořit sítě. Výsledná síť by měla mít
tvar zohledňující oblasti zavíření a vývoj rychlostního profilu v blízkosti stěny, viz obr. 1.16.
obr. 1.16 Výpočetní síť v podélném řezu
Výsledky výpočtu
Pro přehlednost se uvádějí možnosti vyhodnocení, jako jsou vektory rychlosti, profily
rychlosti a vyplněné izočáry.
obr. 1.17 Vektory rychlosti u
14
Úvod do numerického modelování proudění
obr. 1.18 Detail oblasti zavíření (vektory rychlosti u).
obr. 1.19 Vyplněné izočáry statického tlaku
15
Tvorba geometrie a výpočtové sítě
2. Tvorba geometrie a výpočtové sítě
2.1.
Pojem „síť“ a jeho význam pro matematické modelování
Síť představuje systém rozdělení výpočtové oblasti na dílčí na sebe navazující 2D
buňky ve dvoudimenzionálním prostoru nebo 3D buňky ve třídimenzionálním prostoru ([23]).
Lze říci, že výpočtová oblast pokrytá sítí je základem matematického modelování. Neboť
samostatný matematický model (systém matematických vztahů) je pouze „pasivním“
nástrojem, který nabývá smyslu až ve chvíli, kdy je aplikován na konkrétní problém
(výpočtovou oblast pokrytou sítí).
Pokud se hovoří o matematických modelech, které jsou založeny na numerickém
řešení systému parciálních diferenciálních rovnic a vyžadují takto i zadání okrajových
podmínek, lze konstatovat, že možnosti realizování úlohy jsou silně limitovány výkonem
počítačové techniky. Platí zde několik zásad:
·
výpočet je o to náročnější (pomalejší), čím více rovnic je v rámci matematického
modelu do výpočtu zahrnuto (podle náročnosti a komplexnosti modelu);
·
výpočet je o to náročnější, čím více má výpočtová oblast buněk;
·
výpočet je o to náročnější, čím méně kvalitní je síť výpočtové oblasti.
V zájmu přesnosti matematické simulace je nutné provést tomu odpovídající
nastavení matematického modelu. Do různých modelovaných fyzikálních jevů mohou svým
vlivem zasahovat mnohé jevy další. Toto všechno je třeba v nastavení zohlednit. Ovšem
s každým dalším vlivem vstupujícím do výpočtu přibývají také další rovnice, které
matematický model musí řešit. Proto se mohou i při stejně definované výpočtové oblasti i síti
časy výpočtu u různých úloh značně lišit.
Počet buněk patří k hlavním limitujícím faktorů současného matematického
modelování. U mnohých praktických úloh se počty buněk výpočtové oblasti pohybují v řádu
milionů či mnohdy i desítek milionů. Nejsou to zanedbatelná čísla, neboť v každé z buněk je
počítáno mnoho různých veličin. Proto je cílem každého řešitele s ohledem na budoucí čas
výpočtu redukovat počet buněk na nutné minimum. Z hlediska počtu buněk představuje
obrovský nárůst například vytváření tzv. mezních vrstev (viz dále).
Minimalizování počtu buněk by však nemělo být prováděno na úkor kvality sítě.
Kvalitní síť je taková, která se skládá z na sebe navazujících geometricky pravidelných
přibližně stejně velikých a pravidelně po celé výpočtové oblasti rozložených elementů
(buněk). Elementy by měly mít rovněž přiměřenou velikost, aby bylo možné jimi zachytit
v dostatečné míře modelovaný fyzikální děj (například turbulentní vírové struktury a jevy
16
Tvorba geometrie a výpočtové sítě
související s šířením tepla). Z hlediska reálného možného počtu buněk však v praxi dodržení
všech ideálních předpokladů pro tvorbu sítě není většinou možné. Proto se používá
zhušťování sítě v místech, která jsou z hlediska proudění tekutin nebo sdílení tepla pro
řešitele zajímavá nebo pro výpočet stěžejní a naopak použití řidší sítě v místech jiných.
Zvláštním případem zhuštění buněk je vytvoření tzv. mezní vrstvy v blízkosti stěn, která má
za úkol zachytit velké změny fyzikálních veličin u stěny. Zhušťování buněk by mělo být
plynulé. Pokud by byla změna ve velikosti buněk provedena příliš velikou skokovou změnou,
projevilo by se to znatelně na průběhu výpočtu (problémy s konvergencí úlohy) i konečném
výsledku výpočtu (chybný výsledek v daném místě výpočtové oblasti).
2.2.
Gambit, prvky sítě
Funkce Gambitu je mnohem širší, než jen fyzické kreslení oblasti a síťování. Na obr. 2.1
je schématicky znázorněna základní struktura zapojení ostatních programů do tvorby sítě.
FLUENT
Základní struktura programu
FLUENT
GAMBIT
-Definice geometrie
geometrie
síť
jiné CAD/CAE systémy
-2D/3D síť
síť na hranicích
oblasti
2D/3D síť
FLUENT
síť na hranicích
oblasti/v objemu
1. Geometrie, tvorba sítě
2.
3.
5.
6.
7.
8.
Fyzikální model
Okrajové podmínky
Fyzikální vlastnosti
Parametry výpočtu
Výpočet
Vyhodnocení,
interpretace výsledků
TGRID
síť
-2D/3D síť
-hybridní sítě
obr. 2.1 Funkce programu Gambit.
Numerická
metoda
konečných
objemů
je založena
na
vytvoření
systému
nepřekrývajících se elementů, konečných objemů. Původně byla metoda konečných objemů
postavena na konečných objemech tvaru obdélníků a křivočarých čtyřúhelníků ve
17
Tvorba geometrie a výpočtové sítě
dvourozměrném případě a kvádrů nebo obecných šestistěnů v trojrozměrných úlohách (viz
obr. 2.2). Takto vytvořená síť se nazývá strukturovaná síť. Zásadním pravidlem je, že
hranice prvků musí sousedit s jedinou hranicí sousedního elementu, nelze tedy libovolně
zhušťovat síť (je analogií pro metodu konečných diferencí včetně možnosti použití
indexování). Také výsledná výpočtová oblast je pak kvádr nebo obdélník. V současné době
se začíná prosazovat nový přístup, kdy se buduje tzv. nestrukturovaná síť. Konečným
objemem je ve 3D kvádr, čtyřstěn, prizmatický a pyramidový
prvek, jehož výhody byly
ověřeny v úlohách pružnosti, řešených metodou konečných prvků.
kvádr
prizmatický
prvek
čtyřstěn
pyramidový
prvek
obr. 2.2 Tvar konečného objemu
Výše vyjmenované prvky se v současné době mohou kombinovat, čímž se získá optimální
síť, kde v okolí stěny jsou použity čtyřúhelníky a kvádry (pro výpočet z hlediska přesnosti
jsou optimální) a v dalších oblastech, kde nedochází z důvodu existence mezní vrstvy
k velkým gradientů řešených veličin, se použijí zbývající prvky. Ty zajistí snadnou změnu
hustoty sítě, viz obr. 2.3.
obr. 2.3 Použití různých typů prvků
18
Tvorba geometrie a výpočtové sítě
2.3.
Kritéria pro posouzení kvality sítě
Kvalita sítě se posuzuje podle:
·
velikosti buněk (s ohledem na modelovaný děj a požadavek na přesnost výpočtu)
·
vhodnosti uspořádání buněk v prostoru (například zhuštění v místech zajímavých
z hlediska proudění tekutin) s ohledem na konkrétní typ úlohy
·
kvality buněk (nesouměrnost – Skewness, poměr hran (ploch) prvků - Acpect Ratio,
atd.)
Nejvýznamnějším kriteriem pro posouzení kvality buňky je nesouměrnost, kdy se
posuzuje, jak hodně se buňka svým tvarem blíží ideálnímu pravidelnému geometrickému
tvaru v souladu s odpovídajícím schématem sítě. Pokud je buňka jakkoliv deformována, je
její kvalita horší. Obecně se kvalita každé buňky vyjadřuje bezrozměrným číslem v rozsahu 0
– 1, kde 0 znamená výsledek nejlepší a naopak 1 výsledek nejhorší, tedy problematickou
buňku pro výpočty. Tato hodnota se nazývá „míra skosení buňky“ (angl. „skewness
measure“) neboli také míra deformace.
Optimální plocha
(rovnostranná)
Optimální buňka
(rovnostranná)
Aktuální buňka
Kružnice
opsaná
Teoretická obalová
plocha koule
Aktuální plocha
obr. 2.4 Princip posuzování kvality 2D buňky
obr. 2.5 Princip posuzování kvality 3D buňky
pro trojúhelníkové prvky (angl. „tri“)
pro sítě tvořené čtyřstěny (angl. „tet“)
Pro určení kvality 2D buňky, resp. míry její deformace slouží následující vztah:
Skewness measure (TRI ) =
Soptimal - Sreal
Soptimal
( 2.3.1 )
kde Soptimal představuje „optimální plochu buňky“, Sreal „reálnou plochu buňky“ (která může,
ale také nemusí být optimální) a Skewness measure (TRI ) je anglický výraz pro „míru
deformace buňky“ vztahující se ke 2D trojúhelníkovému schématu sítě. U jiných schémat se
vychází z obdobné logiky. Výsledná hodnota by neměla přesáhnout 0.85. Pokud by se tak
stalo, je třeba buňku nebo schéma sítě upravit, aby nebyla ohrožena realizovatelnost a
19
Tvorba geometrie a výpočtové sítě
přesnost výpočtu.
Pro určení kvality 3D buňky odpovídající schématu sítě tvořené čtyřstěny platí vztah
obdobný:
Skewness measure (TET ) =
kde Voptimal
Voptimal - Vreal
Voptimal
,
představuje „optimální objem buňky“, Vreal
( 2.3.2 )
„reálný objem buňky“ a
Skewness measure (TET ) je anglický výraz pro „míru deformace buňky“ vztahující se ke
schématu 3D sítě tvořené čtyřstěny. Výsledná hodnota by neměla přesáhnout 0.9. Pokud by
se tak stalo, je opět třeba buňku nebo schéma sítě upravit.
obr. 2.6– Barevná orientační škála pro vyhodnocení kvality buněk sítě
(0 – nejkvalitnější buňky; 1 – nejméně kvalitní buňky)
Kvalita sítě se dá testovat a počítat v preprocesoru Gambit 2.2.30 pomocí příkazu
„TESTUJ SÍŤ“ (angl. „Examine Mesh“). Po zadání tohoto příkazu je k dispozici
okno poskytující několik možností, jak kontrolovat lokálně či globálně kvalitu sítě v celé
výpočtové oblasti. Výsledkem testu sítě je zbarvení buněk sítě do barevného odstínu, jež
odpovídá podle barevné orientační škály stupni kvality buňky.
2.4.
Příklad vytvoření sítě
Pro ilustraci je prezentováno
vytvoření geometrického modelu
a sítě pro řešení proudění ve
ventilu při maximálním využití
dostupných
prostředků
pro
konstrukci prvků ([24]). Model je
vytvořen v CAD software Solid
Edge (obr. 2.7), který jako takový
nemá na program Fluent, lépe
řečeno
na
program
žádnou
návaznost.
Gambit,
Proto
tento model musel exportovat
se
obr. 2.7 Řez 3D modelem otevřeného pojistného ventilu
(čáry se šipkou značí směr proudu kapaliny)
(uložit) ve formátu, který je pro Gambit srozumitelný a dokáže načíst potřebná data o
geometrii modelu. Proto byl proveden export modelu ze Solid Edge do souboru o formátu
20
Tvorba geometrie a výpočtové sítě
IGES a STEP, což jsou dva nejpoužívanější soubory pro obecné exportování dat, tedy
především geometrie modelu. Po importu do Gambitu překladač tohoto programu bohužel
nebyl schopen načíst celou geometrii. Proto se geometrie přenesla jako mezioperace do
programu Ideas, ze kterého se dále exportovala do IGES formátového souboru a následně
opět do programu Gambit.
Model je řešen jako konstruktérská
součást, což znamená z hlediska
proudění jako „skořápka“ kolem
fluidního jádra, viz obr. 2.8, se
kterým
pracuje
matematický
program Fluent. Nejdříve je nutno
provést separaci tohoto jádra od
skořápky pomocí booleanovských
příkazů. Odstranění všech objemů
součásti,
které
nejsou
obr. 2.8 Řez modelem s vytvořenou
v bezprostředním kontaktu
oblastí proudění
s kapalinou (pružiny, těsnění, plastový obal atd.) se provede opět pomocí booleanovských
funkcí (průnik) a zůstane jen objem vnitřní, tedy objem fluidní (obr. 2.9), ve kterém bude
řešeno proudění.
obr. 2.9 Vnitřní objem pojistného ventilu
obr. 2.10 Vytvořený zjednodušený objem
s naznačenou oblastí pro zhuštění
Geometrie je periodická, je tedy možno vytvořit zjednodušený objem (¼ původního objemu),
dále se zkrátí vstupní přívod a provede se „vyčištění“ geometrie o velmi skloněné hrany, viz
obr. 2.10. Toto zjednodušení pomáhá usnadnit a především urychlit výpočet.
2.4.1. Vytvoření sítě
Vzhledem ke složitosti 3D geometrie byla zvolena síť s tetra/hybridními prvky (3D
čtyřstěny), která je nestrukturovaná s velikostmi hran jednotlivých elementů 0.17 a bylo
21
Tvorba geometrie a výpočtové sítě
vytvořeno zahuštěním kolem ostrých hran užitím size function v místě, jak je naznačeno na
obr. 2.11. S takto vytvořenou sítí se vygenerovalo 1 500 000 buněk.
obr. 2.11 Vytvořená síť s naznačeným zhuštěním kolem zvolených hran
Kvalita vytvořené sítě v této oblasti byla posouzena užitím vyhodnocení kriteria
nesouměrnosti (skewness) a zobrazena v detailu graficky, viz obr. 2.12.
obr. 2.12 Detail sítě s vyhodnocením kvality sítě.
V tomto případě je většina buněk síťované oblasti málo nebo středně deformovaná, avšak
několik buněk se pohybuje v zóně velké deformace. To však nevadí, protože při práci ve
Fluentu bude tato síťová doména převedena na mnohostěnné buňky s jinou kvalitou sítě z
důvodu výrazného urychlení výpočtu.
22
Tvorba geometrie a výpočtové sítě
2.4.2. Definování okrajových podmínek v Gambitu
Gambit je také nástrojem k definování vlastností stěn v souvislosti s okrajovými
podmínkami a definování vnitřních oblastí pro případně vícefázové proudění nebo pro stěny
s přestupem tepla. Tedy lze definovat dva typy okrajových podmínek :
·
okrajové podmínky na hranici
·
podmínky pro oblasti kontinuity (continuum types), specifikace oblasti proudění nebo
pevné látky.
Vytvoření okrajových podmínek se realizuje zadáním vlastností vybraným oblastem modelu
pro další práci, tedy vytvoření tlakového vstupu, výstupu, vytvoření dvou symetrických stěn
(rotační součást) a ostatních stěn, které nemají jinou funkci než jen ohraničit „fluidní“ objem,
viz obr. 2.13.
Tlaková výpusť
Symetrická stěna
vstupní podmínka
Pressure inlet
výstupní podmínka
Pressure outlet
podmínka symetrie
Symmetry 1,2.
zbylé stěny
Wall
Symetrická stěna
obr. 2.13 Okrajové podmínky
Tlakový vstup
Takto vytvořený vysíťovaný model je vyexportován do souboru s příponou .msh a tím je
ukončena všechna práce v programu Gambit.
23
Programový systém Fluent
3. Programový systém Fluent
3.1.
Přehled metod řešení parciálních diferenciálních rovnic
Diferenční metoda je nejstarší všeobecně známou metodou řešení diferenciálních
rovnic, která je použita v ilustrativním příkladu vedení tepla v tyči. Spočívá v nahrazení
derivací diferenčními podíly použitím Taylorova rozvoje, odvozením diferenčních rovnic a
jejich řešením.
Metoda konečných objemů spočívá stručně řečeno ve třech základních bodech
·
dělení oblasti na diskrétní objemy užitím obecné křivočaré sítě
·
bilancování neznámých veličin v individuálních konečných objemech a diskretizace
·
numerické řešení diskretizovaných rovnic
Fluent definuje diskretní konečné objemy užitím non-staggered schematu, kdy všechny
proměnné jsou uchovávány ve středech konečných objemů.
V současné době se začíná prosazovat v řešení proudění také metoda konečných
prvků, která spočívá v těchto bodech
·
násobení diferenciální rovnice bázovými funkcemi
·
dělení oblasti na trojúhelníkové nebo čtyřúhelníkové prvky ve dvourozměrné oblasti (2D)
a čtyřstěny resp. šestistěny ve třírozměrné oblasti (3D)
·
integrace přes konečné elementy založená na variačním principu
·
minimalizace reziduálů
Speciální
metodou
je
spektrální
metoda
vhodná
pro
periodické
proudění
v jednoduchých oblastech (Taylorovy víry vznikající v mezeře mezi koncentrickými válci,
z nichž jeden rotuje).
Další kapitoly jsou výhradně věnovány metodě konečných objemů.
3.2.
Integrace metodou konečných objemů
Integrace diferenciálních rovnic je jednoduše vysvětlena při užití kartézských souřadnic
a pro jednoduchost na rovnicích o jedné prostorové nezávisle proměnné, které si lze
představit jako proudění v trojrozměrném prostoru, kde všechny derivace proměnných ve
směru y a z jsou nulové. Proudění je navíc stacionární (nezávislé na čase).
Rovnice kontinuity je definována takto:
24
Programový systém Fluent
¶u
=0
¶x
( 3.2.1)
rovnice zachování hybnosti
¶
(uu ) = - 1 ¶ p + ¶
¶x
r ¶x ¶x
é ¶uù
ên ¶ x ú + S
ë
û
( 3.2.2)
a rovnice pro přenos skalární veličiny
¶
(uz ) = ¶
¶x
¶x
é ¶z ù
êa z ¶x ú + Sz
ë
û
( 3.2.3)
obr. 3.1 Souřadnicové schéma se speciálním značením buněk pro 1D a 3D model místo
indexů, kde N – sever (North), S – jih (South), E – východ (East), W – západ (West), F –
vpřed (Front), B – vzad (Back)
Integrací těchto rovnic přes konečné objemy se převedou výchozí diferenciální rovnice na
objemový
(
integrál
æ ¶a x ¶a y ¶az
çç
+
+
òòò
¶
¶
¶z
x
y
è
V
(dV=dx.dy.dz,
ö
÷÷dxdydz=
ø
dA=dy.dz),
užitím
divergenčního
òò (a dydz + a dxdz + a dxdz ) ,
x
y
z
teorému
[25]) na plošný (v
S
geometrii velká písmena označují středy konečných objemů a malá písmena hranice, tj.
stěny mezi konečnými objemy, viz obr. 3.1) a diskretizací na výsledný algebraický tvar
následujícím způsobem:
¶u
¶u
ò ¶ x dV = ò ¶ x dxdydz = ò (u )dA = (uA) - (uA)
e
V
V
w
( 3.2.4)
A
Integrace rovnice kontinuity vede na tvar
(uA)e - (uA)w
=0
Fyzikálně výrazy na levé straně označují rozdíl objemových průtoků
Qe - Qw = 0
( 3.2.5)
Integrací rovnice zachování hybnosti se získá
25
Programový systém Fluent
¶
(uu ) = - 1 ¶ p + ¶ éên ¶ u ùú + S Þ
¶x
r ¶ x ¶ x ë ¶ xû
æ u - uP
1
Qeue - Qw uw = - (pe - pw )A + ççn e E
Dxe
r
è
ö
æ u - uW
÷÷ A - ççn w P
Dxw
ø
è
ö
÷÷ A + SDV
ø
( 3.2.6)
a rovnice pro skalární veličinu se upraví shodným postupem na tvar
æ z - zP
z - zW
Qez e - Qw z w = çç a e E
- aw P
Dxe
Dxw
è
ö
÷÷ A + Sz DV
ø
( 3.2.7)
V předchozích rovnicích se používají veličiny a koeficienty jednak definované ve středech
konečných objemů a jednak na stěnách těchto objemů (např. rychlost v rovnici ( 3.2.6)). To
je určitá nevýhoda a je v dalším nezbytně nutné ujednotit ukládání veličin pouze ve středech
konečných objemů. Pokud tato veličina bude určena na stěně, použije se interpolační
schéma pro interpolaci této veličiny do středu buňky. Pro ilustraci je použito nejjednodušší
schéma, tj. aritmetický průměr a diferenční schéma se zjednoduší. Např. rovnici ( 3.2.6) lze
upravit následovně:
uE + uP
u + uW
=
- Qw P
2
2
æ u - uE
u - uW
1
= - (p e - pw )A + ççn e P
-n w P
Dx e
Dx w
r
è
Qe
( 3.2.8)
ö
÷÷ A + SDV
ø
a rovnici pro obecnou proměnnou takto:
Qe
zP +zE
z + zW æ z P -z E
z -zW
- Qw P
= çç a e
- aw P
Dx e
Dx w
2
2
è
ö
÷÷ A + Sz DV
ø
Pak lze pro tuto obecnou rovnici v jednorozměrném případě vyjádřit zP pomocí hodnot
v sousedních konečných objemech následujícími úpravami
æ Qe Qw
æ
æ
Q ö
Q ö
A
A
A ö
A
çç
÷÷z P = çç - a e
+ ae
- aw
- e ÷÷z E + çç a w
+ w ÷÷z W + Sz DV
2
Dx e
Dx w ø
2 ø
2 ø
Dx e
Dx w
è
è 2
è
AP z P = AE z E + AW z W + SC
( 3.2.9)
kde AP = - AE - AW - S P . Předešlou rovnici lze ještě obecně upravit
z P å (- Ai - SP ) = å Ai z i + SC
i
( 3.2.10)
i
kde součet se provede přes sousední buňky (v jednorozměrném případě je i=E, W; v
trojrozměrném případě i=N, S, E, W, F, B,). Ai jsou koeficienty, které obsahují příspěvky od
konvektivních, difúzních a zdrojových členů a SC a SP jsou složky linearizovaných
zdrojových členů a Sz = SC + SP.z P. Použité označení je patrné z obr. 3.1.
26
Programový systém Fluent
Rovnice řešené ve Fluentu jsou rozšířením předchozích na třídimenzionální křivočarý
souřadný systém. Každá iterace sestává z kroků, které jsou zobrazeny diagramem na obr.
3.2.
END
START
řešení rovnice pro zachování
hybnosti
kontrola
konvergence
řešení rovnice kontinuity
(tlaková oprava)
aktualizace rychlosti tlaku
aktualizace vlastností
tekutiny
řešení rovnic pro skalární veličiny
aktualizace turbulentních veličin
aktualizace skalárů
obr. 3.2 Diagram algoritmu řešení Fluentem [15]
a jsou popsány následovně
·
pohybové rovnice pro neznámé složky rychlosti jsou řešeny s užitím hodnot tlaků tak, aby
se aktualizovalo rychlostní pole
·
rychlosti určené v předchozím bodě nemohou splňovat rovnici kontinuity, proto se určují
tzv. tlakové korekce a následně korekce i rychlostního pole
·
pomocí nových hodnot rychlostí se řeší rovnice pro turbulentní energii k a disipaci e
·
řeší se další rovnice pro určení teploty a dalších skalárních veličin
·
aktualizují se fyzikální vlastnosti kapalin (např. viskozita)
·
kontrola konvergence
3.3.
Výběr interpolačního schematu
FLUENT ukládá složky rychlosti a skalární veličiny v geometrických středech konečných
objemů definovaných sítí. Z důvodu výpočtového procesu jsou potřebné hodnoty těchto
veličin na hranicích konečných objemů. Tyto hodnoty jsou získány interpolací, přitom si lze
27
Programový systém Fluent
vybrat mezi následujícími třemi variantami lišícími se řádem přesnosti (vzestupně)
·
mocninová interpolace
·
kvadratická upwind interpolace (QUICK)
·
interpolace druhého řádu/centrální diference
·
QUICK
Při velkých změnách tlaků a průtoků je vhodné rozpočítat úlohu s nejnižším řádem přesnosti
(což je předdefinováno) a po několika iteracích využít vyšší řád přesnosti (pro proudění se
zavířením, s přenosem tepla, disipací apod.)
3.4.
Konvergence.
3.4.1. Residuály
Při simulaci proudění pomocí programu Fluent je velmi důležité získat konvergentní
řešení. Mírou konvergence jsou reziduály, které představují maximum rozdílu dvou
odpovídajících si veličin ve stejném bodě sítě ve dvou po sobě následujících iteracích.
Residuály jsou vyhodnocovány pro všechny počítané veličiny v každém kroku iterace a
zobrazovány pro vybrané veličiny.
i+1-tá iterace
Pi+1
i-tá iterace
Pi
obr. 3.3 Iterace při numerickém stacionárním výpočtu
Měřítkem je přesněji řečeno součet změn počítané veličiny v rovnici pro všechny buňky
v oblasti. Po diskretizaci je odvozena např. rovnice v konzervativním tvaru v jednorozměrné
oblasti pro obecnou proměnnou z ( 3.2.9)
AP z P = AE z E + AW z W + SC
Reziduál R je pak definován součtem přes všechna P
R = å AE z E + AW z W + SC - AP z P
( 3.4.1)
P
28
Programový systém Fluent
R je nenormalizovaný reziduál, který má fyzikální rozměr odpovídající rozměru každého
členu rovnice a číselně se reziduály, např. pro tlak a rychlost, mohou lišit o řády. Proto se
obvykle používá normalizovaný reziduál definovaný následovně:
R =
åAz
E
E
+ AW z W + SC - APz P
P
åAz
P
( 3.4.2)
P
P
Reziduály lze vyhodnocovat graficky v každém kroku iterace, viz obr. 3.4, a snižující se
hodnota reziduálu svědčí o dobře konvergující úloze. Číselně lze sledovat přesné hodnoty
v tab. 3.1.
obr. 3.4 Závislost reziduálu tlaku a entalpie na počtu iterací.
Normalizovanéreziduály
počet iterací
tlak
rychlost u
rychlost v
entalpie
5.00000E+00
2.647l0E-01
3.62767E-01
4.01625E-01
3.49503E-01
1.00000E+01
4.21854E-02
7.31087E-02
6.09502E-02
1.15607E-01
1.50000E+01
1.61787E-02
5.57186E-02
6.77023E-02
6.28091E-02
2.00000E+01
9.91924E-03
4.11899E-02
5.52667E-02
4.24032E-02
2.50000E+01
7.78245E-03
3.71804E-02
5.02612E-02
3.19044E-02
3.00000E+01
6.71127E-03
3.33559E-02
4.61688E-02
2.55360E-02
3.50000E+01
4.96045E-03
3.13033E-02
4.34564E-02
2.11992E-02
4.00000E+0l
6.07668E-03
3.01096E-02
4.09786E-02
1.80783E-02
4.50000E+01
5.21358E-03
2.85215E-02
3.89507E-02
1.56768E-02
5.00000E+01
6.70681E-03
2.67667E-02
3.67708E-02
1.38577E-02
5.50000E+02
5.67326E-04
1.93322E-03
1.06027E-03
1.25243E-04
5.55000E+02
6.90138E-04
1.85336E-03
9.84238E-04
1.18761E-04
5.60000E+02
5.32427E-04
1.78220E-03
9.68897E-04
1.12885E-04
atd.
29
Programový systém Fluent
5.65000E+02
4.20846E-04
1.68717E-03
1.02076E-03
1.07433E-04
5.70000E+02
3.76113E-04
1.62817E-03
1.07420E-03
1.02499E-04
5.75000E+02
3.26542E-04
1.52597E-03
1.07393E-03
9.78981E-05
5.80000E+02
3.00249E-04
1.47025E-03
1.08095E-03
9.36884E-05
5.85000E+02
2.94075E-04
1.31286E-03
1.07175E-03
8.96829E-05
5.90000E+02
2.54086E-04
1.16047E-03
l.05932E-03
8.59168E-05
5.95000E+02
2.33645E-04
1.02260E-03
1.05061E-03
8.23961E-05
6.00000E+02
2.12518E-04
9.10742E-04
1.04899E-03
7.90365E-05
tab. 3.1
Dále lze vyhodnotit, ve kterém bodě oblasti je nejvyšší hodnota reziduálu. Reziduály slouží
k vyhodnocení konvergence. Obecně řešení velmi dobře konverguje, když se normalizované
reziduály snižují řádově k hodnotě 1.10-3 a entalpii k hodnotě řádově 1.10-6.
3.4.2. Urychlení konvergence
Konvergence je ovlivněna mnoha faktory, jako je počáteční odhad, velký počet
buněk, relaxační faktor atd.
Pro urychlení konvergence se navrhuje využít počátečního odhadu proměnných
významných pro proudění, což je nejlepší způsob, jak začít řešit úspěšně úlohu. V opačném
případě jsou všechny veličiny definovány inicializací, často jsou pokládány rovny nule na
počátku výpočtu. Nejvýznamnější příklady nastavení počátečních podmínek jsou:
·
teplota pro problémy řešící přenos tepla při užití stavové rovnice
·
rychlost při velkém počtu buněk
·
teplota i rychlost při řešení přirozené konvekce
·
proudění s reakcí, kdy je dobré nastavit teplotu i hmotnostní podíly.
Důležitou technikou k urychlení konvergence je technika step by step (postupně od
jednoduché úlohy ke složitější). Při řešení problému s přenosem tepla je dobré začít výpočet
z izotermního proudění, při řešení reagujícího proudění z proudění bez reakce se zahrnutím
příměsí. Problém se nadefinuje nejprve celý a teprve potom se vyberou proměnné, pro které
se vyřeší počáteční stav.
3.4.3. Relaxace
Z důvodu nelinearity diferenciálních rovnic není obecně možné získat hodnoty všech
proměnných
řešením
původně
odvozených
aproximačních
diferenčních
schémat.
Konvergence lze však dosáhnout užitím relaxace, která redukuje změny každé proměnné
30
Programový systém Fluent
v každé iteraci. Jednoduše řečeno, nová hodnota z P,i +1 v konečném objemu obsahujícím
bod P závisí na staré hodnotě z předešlé iterace z P,i , nové hodnotě z aktuální iterace
z P ,i +1,vyp (resp. vypočtené změně Dz P = z P ,i +1,vyp - z P ,i a relaxačním parametru a Î 0,1
následovně
z P ,i +1 = az P .i +1,vyp + (1 - a )z P .i
( 3.4.3)
resp.
z P ,i +1 = z P ,i + aDz P
z
z P ,i +1,vyp
Dz P = z P ,i +1,vyp - z P ,i
z P ,i + 1
z P,i
0
a Î 0,1
1
a
obr. 3.5 Vysvětlení relaxačního parametru
Tyto relaxační parametry se mohou nastavit pro všechny počítané proměnné. Zvláště pro
rychlosti se nastavují velmi malé, řádově desetiny až setiny. Přitom je vhodné během
výpočtu tyto hodnoty měnit a tím urychlovat konvergenci, tzn. jestliže změny reziduálů jsou
velké při přechodu od jedné iterace k druhé, nastaví se malý relaxační faktor a tím se tlumí
nelinearity, pokud se změny reziduálů stávají konstantní, je vhodné relaxační faktory zvětšit.
31
Turbulentní proudění skutečných kapalin
4. Turbulentní proudění skutečných kapalin
4.1.
Klasifikace proudění skutečných kapalin
Proudění skutečných kapalin může být klasifikováno jako laminární nebo turbulentní
proudění.
U turbulentního proudění bylo na základě experimentálních měření zjištěno, že na
stěnách potrubí nebo obtékaného tělesa vzniká vrstva kapaliny s laminárním pohybem, tzv.
laminární podvrstva, jejíž tloušťka je několik desetin milimetrů. Těsně za laminární
podvrstvou je přechodová vrstva mezi laminární podvrstvou a turbulentním jádrem, které
tvoří další oblast turbulentního proudu. Laminární podvrstva a přechodová vrstva tvoří
turbulentní mezní vrstvu. Uvažujme nejjednodušší případ - tenkou desku paralelní s proudem
tekutiny, viz obr. 4.1. Tlak je v celém objemu tekutiny konstantní. Tekutina na desce lpí,
protože vlivem viskozity se zabrzdí nejbližší vrstvy tekutiny u povrchu desky. Rychlost
tekutiny s odlehlostí od stěny narůstá až na hodnotu rychlosti nenarušeného proudu v ¥ .
Tloušťka "zabržděné" tekutiny, tj. tloušťka mezní vrstvy je u náběžné hrany nulová a na
odtokové hraně je maximální. V mezní vrstvě a oblasti kolem desky nejsou proudnice
paralelní přímky, ale tvoří mírně se rozbíhající svazek. Složka rychlosti kolmá k desce je
mnohem menší a lze ji zanedbat.
obr. 4.1 Mezní vrstva při obtékání desky [22]
V přední části je mezní vrstva laminární, v zadní turbulentní, mezi nimi přechodová oblast.
Okamžitá hranice turbulentní mezní vrstvy – plná nepravidelná křivka - se s časem mění.
Střední tloušťka turbulentní mezní vrstvy je zakreslena čárkovaně. Kritérium pro stanovení
přechodu laminární mezní vrstvy na turbulentní je opět kritické Reynoldsovo číslo, jehož
hodnota se mění se stupněm turbulence proudu. Zpravidla se udává Re k =
v ¥ xk
= 5.105 ,
n
kde x k je vzdálenost od náběžné hrany, ve které laminární mezní vrstva přechází do
turbulentní. Je vidět, že stanovení typu proudění není zcela jednoduché a jednoznačné a
záleží na zkušenostech řešitele.
32
Turbulentní proudění skutečných kapalin
V případě jednorozměrného proudění v potrubí je přechod k turbulenci dán
experimentálně určeným kritickým Reynoldsovým číslem Re , definováným vztahem
v sd
, kde v s je střední rychlost v potrubí, d jeho průměr a n kinematická viskozita.
n
Re =
Kritická hodnota Rekrit pro potrubí kruhového průřezu je 2320. Při Re £ Rekrit se v potrubí
vyvine uspořádané laminární proudění, pohyb se děje ve vrstvách a částice tekutiny se
nepohybují napříč průřezem. Je-li
Re ³ Rekrit , proudění je turbulentní. Při vyšších
Reynoldsových číslech částice tekutiny konají neuspořádaný pohyb všemi možnými směry.
Tento pohyb je nepravidelný, náhodný a připomíná pohyb molekul plynu, ale na rozdíl od
molekul se částice tekutiny mohou rozpadat a ztrácet tak svou identitu. Pohyb částic kolmo
ke stěně zvyšuje tok hybnosti ke stěně a proto je pokles tlaku ve směru proudění mnohem
větší než u laminárního proudění. Následkem promíchávání tekutiny jsou rozdíly rychlosti na
různých místech průřezu mnohem menší než u laminárního proudění mimo oblast poblíž
stěny.
Proudění lze vizualizovat různými metodami a pozorovat odlišnosti laminárního a
turbulentního proudění, viz obr. 4.2. U turbulentní mezní vrstvy lze definovat turbulentní
(koherentní) vírové struktury charakteristické právě pro turbulentní proudění.
obr. 4.2 Srovnání laminární a turbulentní mezní vrstvy [5]
Vlastnosti turbulentního proudění jsou:
·
náhodný pohyb částic tekutiny, tedy objemů obsahujících velké množství molekul,
přičemž pohyb částic se skládá z uspořádaného středního pohybu a z náhodných
fluktuací, z čehož vyplývá analogie mezi chováním molekuly (Brownův pohyb
molekul) a chováním částice tekutiny. Vlivem fluktuací se může dostat molekula
z oblasti větší makroskopické rychlosti do oblasti menší makroskopické rychlosti a při
nárazu na jinou molekulu se zpomalí, přičemž molekulu, na niž narazila, zrychlí a
odevzdá jí část své hybnosti. Opačně je tomu, přechází-li molekula z oblasti menší
rychlosti do oblasti větší rychlosti, kdy se její hybnost při nárazu zvětší. Tím dochází
ke sdílení hybnosti mezi oblastmi tekutiny s rozličnou rychlostí, což se projevuje
rostoucím odporem proti proudění jako vnitřní tření tekutiny.
33
Turbulentní proudění skutečných kapalin
·
tečné napětí, vznikající u turbulentního proudění, není určeno pouze vnitřním třením
v tekutině a rychlostním gradientem jako tomu je u laminárního proudění (Newtonův
zákon t = h
dv
), ale změnou hybnosti makroskopických částeček, jako následek
dy
jejich pronikání mezi sousední vrstvy. Tento neuspořádaný pohyb vyvolá tzv.
přídavné turbulentní napětí.
·
turbulentní viskozita, o níž nelze mluvit jako o fyzikální konstantě tekutiny, jako
tomu je u molekulové viskozity laminárního proudění, ale jako o složité funkční
závislosti stavu proudící tekutiny a poloze uvažovaného bodu, tedy sdílení hybnosti
fluktuacemi a odlehlosti od stěny. Proto rychlostní profil u turbulentního proudění ve
srovnání s laminárním je více plochý (nemá parabolický charakter).
·
difuzivní charakter turbulence, kdy gradienty rychlosti vyvolané turbulentními
fluktuacemi rychlostí jsou zdrojem vazkých napětí a disipace energie. Zvyšuje se tak
vnitřní energie tekutiny na úkor kinetické energie turbulence. Turbulence proto
potřebuje trvalý přísun energie ke krytí těchto ztrát, jinak rychle zaniká.
4.2.
Turbulentní proudění
Proudění se obecně nazývá turbulentní, jestliže jeho proměnné vykazují chaotické
fluktuace jak v prostoru, tak v čase, viz obr. 4.3.
obr. 4.3 Plně vyvinuté turbulentní proudění – rychlost jako funkce času
První práce z oblasti turbulentního pohybu tekutin, jednoho z fascinujících jevů
přírody, se datují od r.1883 a jejich autorem je Osborn REYNOLDS, viz [1]. Rovnice popisující
takové proudění jsou známy již desítky let. Bohužel problém turbulence z hlediska fyziky
není stále vyřešen. Ačkoliv byl v současné době udělán významný pokrok, zvláště v oblasti
nelineárních dynamických systémů nebo teorie chaosu, úplné řešení turbulence nelze
v blízké budoucnosti očekávat. Avšak zájem o turbulenci není inspirován pouze přáním
34
Turbulentní proudění skutečných kapalin
porozumět její podstatě, ale nutností předpovídat turbulentní proudění v mnoha technických
aplikacích. Navzdory náhodnosti turbulence detailní studie ukazují, že turbulentní proudění
sestává z prostorových struktur, které se obvykle nazývají „eddies“, (turbulentní víry), viz
obr. 4.4.
a) Mezní vrstva turbulence
b) Turbulentní vírová cesta
generovaná v blízkosti stěny [2]
generovaná obtékáním překážky [2]
obr. 4.4 Příklady typických turbulentních struktur
Je snahou charakterizovat turbulenci pomocí těchto struktur, aby bylo možno vysvětlit
dynamiku turbulence při vzniku, vývoji a zániku vírů „eddies“ jako funkci času. Je zřejmé, že
tento výzkum závisí na možnostech získat informace o prostorových strukturách turbulence a
jejich vývoji v čase.
Jak bylo prezentováno dříve, rovnice proudění tekutin jsou dobře známy. Rychlý rozvoj
výpočetní techniky v posledních letech umožňuje řešit tyto rovnice přístupem, který se
nazývá numerická simulace, což je jeden z nástrojů studia základních aspektů turbulence.
Její hlavní výhodou je, že dává detailní informace o trojdimenzionálních strukturách, které
nelze získat měřením v laboratoři.
4.2.1. Teorie turbulence
Jak bylo řečeno v úvodu, turbulentní proudění obsahuje prostorové struktury, nazývané
„eddies“, tj. turbulentní víry různých velikostí (viz obr. 4.4). Velké víry obsahující energii se
rozpadají na menší. Tento kaskádní proces je ukončen disipací energie nejmenších vírů na
teplo, viz obr. 4.5.
Turbulentní víry jsou charakterizovány délkovým měřítkem l [m] (tj. geometrií oblasti
resp. charakteristickým rozměrem) a rychlostním měřítkem u [m.s-1]. V následujících
kapitolách budou tato měřítka označována jako makroměřítka.
Na druhé straně tekutina, ve které se objevuje turbulence, je charakterizována svými
molekulárními vlastnostmi, jako je kinematická viskozita n [m2.s-1]. Hlavním důsledkem
35
Turbulentní proudění skutečných kapalin
kinematické viskozity je vyhlazení rychlostních gradientů pomocí molekulární difúze.
Výše zavedené parametry dovolují zavést bezrozměrnou veličinu, známou jako
Reynoldsovo číslo:
Re l =
u.l
n
( 4.2.1)
které lze upravit následovně
l2
u.l l n Tn
Re l =
=
=
l
n l
Tt
u
( 4.2.2)
kde Tt označuje časové měřítko přenosu turbulentních vírů o makroměřítku l a Tn označuje
časové měřítko molekulární difúze. Proudění lze charakterizovat na základě hodnoty
Reynoldsova čísla následovně [6], jestliže
·
Tn áTt tj. Re l á1
laminární
proudění,
procesy
molekulární
difúze
převažují
a
turbulentní víry zanikají.
·
Tn ñTt tj. Re l ñ1
turbulentní proudění, turbulentní víry přetrvávají. Nerovnost je
splněna pro poměrně malé hodnoty parametrů proudění, takže lze udělat závěr, že
většina proudění je turbulentní.
·
Tn ññTt tj. Rel ññ1 plně vyvinutá turbulence, což znamená, že viskózní děje, které
ovlivňují časové měřítko Tn mohou být zanedbány vzhledem k dynamice vírů, které se
objevují nad hodnotou Tt . Turbulentní víry v plně vyvinutém turbulentním proudění jsou
téměř neviskózní.
·
Tn » Tt
,
tj. Re l » 1 přechodový stav, laminární stacionární proudění se mění na
turbulentní, nestacionární, pokud je překročeno kritické Reynoldsovo číslo Rel. Proudění
se zpočátku stává periodické. Tato kvalitativní změna v chování proudění se nazývá
bifurkace, Při zvyšování Reynoldsova čísla se vytvářejí další nestability, až se proudění
stane plně turbulentní.
Důležitým důsledkem Reynoldsovy podobnosti je disipace, resp. rychlost disipace,
při které turbulentní víry ztrácejí svou kinetickou energii a mění ji na teplo. Disipace (pro
jednotku hmotnosti) e [m2s-3] je definována vztahem
e =
u3
l
( 4.2.3)
Výsledkem rozměrové analýzy je délkové mikroměřítko h těchto disipačních oblastí, které
se nazývá Kolmogorovovo měřítko
36
Turbulentní proudění skutečných kapalin
1
æn 3 ö 4
h = çç ÷÷
èe ø
( 4.2.4)
Lze snadno odvodit vztah
3
3
h
= Rel 4 Þ h = Rel 4 .l
l
( 4.2.5)
12
10
1
10
0.1
0.01
l , h [m]
l , h [m]
8
generace
makroměřítko
0.001
0.00001
6
0.000001
1.E+02
mikroměřítko
přenos
0.0001
disipace
1.E+03
1.E+04
1.E+05
1.E+06
1.E+07
Re [1]
4
2
0
1.00E+03
1.00E+04
1.00E+05
Re [1]
1.00E+06
v
1.00E+07
obr. 4.5 Mikroměřítko a makroměřítko turbulentních vírů v dekadických a logaritmických
souřadnicích
Turbulentní proudění sestává ze spojité oblasti vírových struktur, jejichž délková měřítka
(rozměry) leží mezi l a h, viz obr. 4.5. Rozlišují se dva typy vírových struktur, struktury
v blízkosti stěny (vlásenkové víry - hairpin, pukání - bursts, proužky - streaks) a struktury
uprostřed proudu.
obr. 4.6 Turbulentní víry v proudu tekutiny, jejich štěpení, přeměna v teplo
37
Turbulentní proudění skutečných kapalin
4.2.2. Metody matematického modelování turbulentního proudění
MODELY TURBULENCE
Prostorová filtrace
LES
Časové středování
Modely
turbulentní
viskozity
DES
Přímá simulace
Reynoldsův
napěťový
model - RSM
DNS
Modelování turbulence je stále ve stádiu výzkumu a vývoje, který se neustále mění
s pokrokem v matematickém, fyzikálním a technickém odvětví. Při numerické simulaci
turbulentního proudění existují tři
teoreticky odlišné přístupy,
které
vyplývají ze
zjednodušujících modifikací výchozích rovnic popisujících proudění ([4]).
·
Metoda přímé simulace (DNS-Direct Numerical Simulation)
·
Metoda velkých vírů (LES-Large Eddy Simulation)
·
Metoda časového středování (RANS-Reynolds Averaged Navier-Stokes equations)
u
u
u
DNS
LES
RANS
obr. 4.7 Metody modelování turbulence
Metoda přímé numerické simulace je dána velkými nároky na kapacitu počítače z
důvodu velmi jemné sítě. Počet uzlových bodů sítě pro DNS lze odhadnout řádově z
38
Turbulentní proudění skutečných kapalin
Kolmogorovova mikroměřítka turbulence (rozměr nejmenších turbulentních vírů) Np » Rel9 / 4.
Počet uzlových bodů sítě tedy prudce narůstá s Reynoldsovým číslem, což vede k technické
nereálnosti výpočtů při stávající výpočetní technice.
1E+16
1E+15
1E+14
1E+13
Np [1]
1E+12
1E+11
1E+10
1E+09
1E+08
10000000
1000000
1.E+02
1.E+03
1.E+04
1.E+05
1.E+06
1.E+07
Re [1]
obr. 4.8 Odhad počtu uzlů sítě pro metodu DNS
Metoda velkých vírů
(LES-Large Eddy Simulation) je založena na modelování
velkých vírů, jako prostorových časově závislých útvarů, které lze zachytit sítí. Turbulentní
víry o malých měřítcích se málo se podílejí na transportních jevech, ale jejich prostřednictvím
dochází k disipaci kinetické turbulentní energie v důsledku viskozity na teplo. Tyto malé víry
jsou parametrizovány tzv. subgridními modely a odstraněny pomocí filtrace turbulentního
pole. Volbou šířky pásma filtru, většinou odpovídajícího rozměru buněk sítě, je možné
dosáhnout takový počet buněk sítě, který lze se současnou výpočetní technikou řešit.
Pro většinu inženýrských úloh turbulentního proudění zůstávají nejpoužívanějším
nástrojem statistické modely turbulence, které jsou založeny na metodě časového
(Reynoldsova) středování (RANS-Reynolds Averaged Navier-Stokes equations) veličin
turbulentního proudění a na následující proceduře časového středování bilančních rovnic.
4.2.3. Reynoldsova rovnice
Turbulentní proudění se vyznačuje náhodným charakterem, ale je statisticky stabilní.
V zásadě je možno takové proudění řešit i pomocí Navierových - Stokesových rovnic při
použití statistické metody časového středování ([12]), které činí úlohu technicky
39
Turbulentní proudění skutečných kapalin
zvládnutelnou. Dle O. Reynoldse (1895) okamžité hodnoty veličin popisujících turbulentní
proudění lze tedy rozložit na část časově středovanou a fluktuační složku (viz obr. 4.9),
přičemž platí
z = z +z¢
kde
z =
1T
z dt
T ò0
resp. z =
z ¢=0
1
åz i
N i
obr. 4.9 Fluktuace a časove středovaná část
Platí Reynoldsova pravidla (podrobně viz. 14.2).
z = z , z + z ¢ = z + z ¢ = z , z z ¢= 0, z + y = z + y ,
zy = z .y + z ¢y ¢,
( 4.2.6)
¶z ¶z
=
¶x ¶x
kde z ¢y ¢ je korelační moment fluktuačních složek.
Po dosazení součtu časově středované a fluktuační hodnoty do rovnice kontinuity se
dostane:
¶ (u j + u ¢j )
¶ xj
¶u j
¶ xj
+
¶u ¢j
¶ xj
= 0
( 4.2.7)
=0
Po časovém středování platí
¶uj
¶ xj
+
¶ u ¢j
¶ xj
=0
( 4.2.8)
Rovnice kontinuity pro středovanou hodnotu má tvar:
¶uj
¶ xj
=0
( 4.2.9)
Rovnice kontinuity pro fluktuační složku se dostane odečtením rovnice ( 4.2.9) od rovnice (
4.2.8):
¶u ¢j
¶ xj
=0
( 4.2.10)
Obdobně lze dosadit i do Navierovy - Stokesovy rovnice za předpokladu, že f i = f i , r = r :
40
Turbulentní proudění skutečných kapalin
¶ (u i + u ¢i ) ¶ (u i + u i¢ )(u j + u ¢j )
¶ 2 (u i + u i¢ )
1 ¶ (p + p ¢ )
=+n
+ fi
+
¶t
¶ xj
r ¶ xi
¶x 2j
( 4.2.11)
Výsledkem je Reynoldsova rovnice podobná formálně Navierově-Stokesově rovnici pro
středované veličiny s tím, že obsahuje navíc člen na levé straně:
¶ u i ¶ ui u j
¶ 2ui
1 ¶p
¶
¢
¢
+
+
+n .
+ fi
ui u j = ¶t
¶ xj
¶ xj
r ¶ xi
¶ x 2j
( 4.2.12)
Rovnice pro fluktuační složky rychlosti je
¶ (u i¢ ) ¶u i u ¢j ¶u ¢i u j ¶ (u i¢u ¢j )
¶ 2 (u i¢ )
1 ¶ (p ¢)
=+n
+ fi
+
+
+
¶t
¶x j
¶x j
¶ xj
r ¶ xi
¶x 2j
( 4.2.13)
- r u ¢i u ¢j jsou tzv. Reynoldsova (turbulentní) napětí, která existují jen při turbulentním
proudění. Turbulentní vír může směšovat tekutinu o různých rychlostech v objemu tvaru
krychle, viz obr. 4.10a).
b)
a)
e)
-ru1¢u3¢
c)
- r u 3¢ u1¢
d)
obr. 4.10 Deformační účinky Reynoldsových napětí [12]
Jestliže tato tekutina o odlišné rychlosti obtéká jednu plochu a protilehlou plochu krychle
neovlivňuje, pak se krychle deformuje z důvodu odchylek rychlostí na těchto dvou plochách,
viz. obr. 4.10 b). Rychlost, kterou je tekutina o odlišné rychlosti transportována přes plochu
krychle, je tok hybnosti. Účinek tohoto toku na krychli je shodný s účinkem síly na plochu
krychle, která se deformuje. Tedy turbulentní tok hybnosti působí jako napětí. V tomto
případě, tekutina pohybující se svislou fluktuační rychlostí u3´ se míchá s tekutinou
41
Turbulentní proudění skutečných kapalin
pohybující se vodorovnou fluktuační rychlostí u1´. Výsledkem je složka Reynoldsových napětí
- r u1¢u 3¢ . Dokonce, jestliže se uvažuje pouze jedna plocha krychle, tekutina pohybující se
v jednom ze tří souřadnicových směrů může být v ní míchána a výsledkem jsou deformace
u1¢u1¢, u1¢u ¢2 , u1¢u 3¢ . Dále předpokládejme, že rychlost u1´ je turbulentně míchána na horní ploše
původní krychle, obr. 4.10 c), s předpokládanou rychlostí u3´. Výsledná deformace souvisí
s Reynoldsovým napětím - r u3¢ u1¢ , obr. 4.10 d). Deformace je identická, jen má jinou
orientaci. Na obr. 4.10 e) jsou zobrazeny tyto dvě deformace bez ohledu na rotaci. Závěrem
lze říci, že u 3¢ u1¢ = u1¢u 3¢ . Stejný počet kombinací vznikne pro plochy kolmé k dalším dvěma
souřadnicovým směrům, získá se devět složek Reynoldsových napětí, podobně jako u
viskózních napětí, a také platí symetrie.
é u1¢u1¢ u1¢u ¢2
ê ¢ ¢
êu 2 u1 u 2¢ u ¢2
êëu 3¢ u1¢ u 3¢ u ¢2
u1¢u 3¢ ù éu1¢u1¢
ú ê
u 2¢ u 3¢ ú = êu1¢u 2¢
u 3¢ u 3¢ úû êëu1¢u 3¢
u1¢u 2¢
u 2¢ u 2¢
u 2¢ u 3¢
u1¢u 3¢ ù
ú
u ¢2 u 3¢ ú
u 3¢ u 3¢ úû
Turbulentní tok hybnosti působí tedy jako napětí a je nazván Reynoldsovo napětí, pro které
lze odvodit také transportní rovnice. Vyjde se z ( 4.2.13), ve které se formálně zamění sčítací
index j za k (j se vyskytuje v proměnné ui¢u¢j ).
¶ (u ¢i ) ¶u i u k¢ ¶u ¢i u k ¶ (u ¢i u ¢k )
¶ 2 (u ¢i )
1 ¶ (p ¢)
+
+
+
=+n
¶t
¶x k
¶x k
¶ xk
r ¶ xi
¶x k2
( 4.2.14)
Tato rovnice se násobí u ¢j , dále rovnice ( 4.2.14) se napíše pro složku j a násobí u i¢ . Takto
získané rovnice se sečtou, časově středují a dostanou se transportní rovnice pro
Reynolsdova napětí:
¶ u ¢i u ¢j
¶ u ¢i u ¢j
¶ é
+ uk
=(u ¢i u ¢j uk¢ ) + p¢ (d kj u ¢i + d ik u ¢j ) -n ¶ (u ¢i u ¢j )ùú ê
¶t
¶x k
¶x
r
¶x k
14k4ë444444424444444
443û
difúzní transport
¶u j
é
é ¶u i¢ ¶u ¢j ù
¶u i ù p¢ é ¶u i¢ ¶u ¢j ù
- êu i¢u k¢
+ u ¢j u ¢k
+
ú - 2n ê
ú+ ê
ú
¶x k
¶x k û r êë ¶x j ¶x i úû
¶x k ¶x k û
ë
ë
14444244443 1442443 14
4244
3
produkce
redistribuce
( 4.2.15)
disipace
Reynoldsova napětí tvoří tenzor o devíti členech, přitom nezávislých je šest, proto je i rovnic
šest, což tvoří rozsáhlý systém diferenciálních rovnic obtížně řešitelných. Proto bude
pozornost věnována teoriím, zabývajícím se jednodušším vyjádřením Reynoldsových napětí
v rovnici ( 4.2.12) (tzv. modely turbulence).
42
Turbulentní proudění skutečných kapalin
4.2.4. Boussinesquova hypotéza o vírové (turbulentní) viskozitě
Základem matematických modelů turbulence, hlavně těch jednodušších, je popis
lokálního
stavu
turbulence
vírovou
(turbulentní)
viskozitou
vyjádřenou
pomocí
rychlostního měřítka u a délkového měřítka l
m t » l .u
( 4.2.16)
Úkolem jednotlivých modelů turbulence je pak vyjádřit turbulentní napětí a toky tepla nebo
jiných skalárních veličin pomocí zvoleného měřítka a určit rozložení tohoto parametru
v proudovém poli. Většina modelů přitom využívá Boussinesqovy hypotézy o vírové
(turbulentní) viskozitě. Tato hypotéza předpokládá, že podobně jako při laminárním proudění,
kdy platí v zjednodušeném dvourozměrném proudění pro smykové napětí Newtonův vztah,
jsou turbulentní napětí a turbulentní toky úměrné gradientu střední rychlosti, teploty,
koncentrace apod., tj .
laminární proudění
Boussinesquova
turbulentní proudění
molekulová viskozita
hypotéza (analogie)
vírová turbulentní viskozita
t =m
du
dy
t t = - r u ¢v ¢ = m t
¶u
¶y
Obecně
æ ¶u ¶ uj
- r ui¢u ¢j = m t çç i +
è ¶ x j ¶ xi
æ ¶u
¶uj
- u ¢i u ¢j = n t ç i +
ç¶ x
¶ xi
j
è
ö 2
÷ - rkd ij resp.
÷ 3
ø
ö 2
÷ - kd ij
÷ 3
ø
( 4.2.17)
kde k je turbulentní kinetická energie
k=
1
u ¢j u ¢j
2
Podrobněji
k=
( 4.2.18)
(
)
1
1
u ¢j u ¢j = u1¢ 2 + u 2¢ 2 + u 3¢ 2 . Na rozdíl od laminárního proudění
å
2 j =1..3
2
turbulentní viskozita není fyzikální vlastností kapaliny, ale proudění. Je silně závislá na míře
turbulence a může se výrazně lišit v rámci proudového pole. Pohybová rovnice ( 4.2.12) se
upraví použitím Boussinesquovy hypotézy takto:
¶ ui ¶ ui u j
¶ 2u i
¶ 2u i
1 ¶p
+
=+n
+n t
+ fi
¶t
¶ xj
r ¶ xi
¶ x 2j
¶ x 2j
( 4.2.19)
Rovnice pro přenos tepla resp. jiného skaláru se upraví způsobem shodným
s odvozením rovnice pro přenos hybnosti ( 4.2.12) pro středované hodnoty. Po úpravě se
v těchto rovnicích objeví člen - u i¢z ¢ reprezentující turbulentní tok tepla resp. jiné skalární
43
Turbulentní proudění skutečných kapalin
veličiny a podobně jako ( 4.2.17) se vyjádří
- ui¢z ¢ = Gt
¶z
¶ xj
( 4.2.20)
kde Gt turbulentní difusivita. Většina modelů předpokládá, že turbulentní difuzivita je úměrná
kinematické turbulentní viskozitě podle vztahu Gt =
nt
, kde st je turbulentní Prandtlovo
st
resp. Schmidtovo číslo pro přenos tepla resp.koncentrace.
Z hlediska modelování turbulentní viskozity v proudovém poli lze rozdělit modely
turbulence do tří skupin, a to nularovnicové (algebraické), jednorovnicové a dvourovnicové
modely, nazvané podle počtu doplňujících diferenciálních rovnic.
44
Statistické modely turbulence
5. Statistické modely turbulence
Základní problém výpočtu turbulentního smykového proudění spočívá v přítomnosti
Reynoldsova napětí v rovnicích popisujících střední pohyb tekutiny, takže systém
pohybových rovnic není uzavřen jako v případě laminárního proudění. Soubor přídavných
rovnic a empirických vztahů, které společně s pohybovými rovnicemi tvoří řešitelný systém
rovnic, se nazývá modelem turbulence. Modely turbulence lze rozdělit do několika skupin dle
obr. 5.1.
Laminární
proudění
Přímá metoda
DNS
Matematické
modely
proudění
Přímá metoda
Turbulentní
proudění
Metoda velkých
vírů - LES
Metoda časového
středování RANS
Metoda
Reynoldsových napětí
Boussinesquova
hypotéza
Nularovnicový
model
k-e
Boussinesquova
Jednorovnicový
model
hypotéza
RNG k- e
Dvourovnicový
model
k-e available
modeblel
k-w
obr. 5.1 Schéma metod řešení proudění.
45
Statistické modely turbulence
5.1.
Model směšovací délky (nularovnicový model) - Prandtl
První model popisující rozložení turbulentní viskozity n t navrhl Prandtl. Turbulentní
viskozita je vyjádřena v závislosti na střední hodnotě rychlosti opět ve zjednodušené
dvourozměrné podobě vztahem
m t = rn t = rl m
2
¶u
¶y
( 5.1.1)
kde lm je směšovací délka, která se stanoví na základě empirických vztahů l m = ky , kde k
je von Kármánova konstanta (0.41). Tento model je vhodný pro modelování proudění ve
smykové vrstvě. Nedostatkem je, že se předpokládá lokální rovnováha, tj. produkce
turbulentní kinetické energie (TKE) je rovna rychlosti disipace TKE. Nepostihuje transport
turbulence.
5.2.
Jednorovnicový model
Aby bylo možné postihnout transport turbulentních parametrů, je nutné řešit pro tyto
parametry diferenciální
transportní rovnici. Nejjednodušší modely používají transportní
rovnici pro rychlostní měřítko turbulentního pohybu
k , kde k =
(
)
1 2
1
u1¢ + u 2¢2 + u3¢ 2 = u ¢j 2
2
2
je kinetická (středovaná) energie turbulentního pohybu vztažená na jednotku hmotnosti. Pro
k lze odvodit exaktní rovnici z Navierových - Stokesových rovnic ve tvaru:
p ¢ öù
¶ k ¶ uj k
¶ é æ u ¢l u l¢
¶ 2k
¶ ul
¶ u ¢l ¶ u ¢l
+
=+ d jl ÷÷ú + n
- u l¢u ¢j
-n
êu ¢j çç
2
¶t
¶ xj
¶ xj ë è 2
r øû
¶x j
¶ xj
¶ xj ¶ xj
{
123
123
1424
3 142
43
1
4
4
4
4
2
4
4
4
4
3
I
II
kde
III
IV
V
( 5.2.1)
VI
I – rychlost změny
II – konvektivní transport
III – turbulentní difúze v důsledku fluktuací rychlosti a tlaku
IV – molekulová difúze
V – produkce v důsledku smykového tření
VI – vazká disipace
Na pravé straně se objevují členy vyjadřující turbulentní difúzi v důsledku fluktuací rychlosti a
tlaku, dále produkci k v důsledku interakce Reynoldsových napětí a gradientu střední
rychlosti a disipaci e v důsledku přeměny energie na energii tepelnou. Působí-li v oblasti
rovněž Archimedovy síly, je na pravé straně rovněž člen odpovídající produkci (destrukci)
kinetické energie v důsledku vztlakových sil. V takto odvozené rovnici se vyskytují neznámé
korelace v difúzním a disipačním členu. Aby se získala uzavřená soustava rovnic, je nutné
46
Statistické modely turbulence
tyto členy modelovat pomocí vztahu
æ u ¢u ¢ p¢ ö n ¶ k
,
- u ¢j çç l l + ÷÷ = t
r ø sk ¶ xj
è 2
e =n
¶u l¢¶ul¢
k 3/2
=
C
D
l
¶x 2j
( 5.2.2)
kde s k a CD jsou empirické konstanty. Po dosazení za tyto členy do rovnice pro k a
s využitím vztahů ( 4.2.17) a ( 4.2.20) pro n t a Gt má rovnice pro k tvar:
¶ k ¶ ujk
¶
+
=
¶t
¶ xj
¶ xj
3/ 2
æ nt ¶ k ö
æ¶u
ö
ç .
÷ + n t ç j + ¶ u l ÷ ¶ u l - CD k
çs ¶ x ÷
ç¶ x
¶ x ÷ø ¶ x j 1
l3
j ø
è k
1è44l 424j4
43 424
P
( 5.2.3)
e
a podle Kolmogorova-Prandtlova vztahu je
( 5.2.4)
n t = Cn k l
je empirická konstanta. Délkové měřítko l charakterizující turbulentní pohyb je
kde Cn
definováno pomocí empirických vztahů podobně jako v případě modelů směšovací délky.
Jednorovnicové modely postihují transport turbulence a jsou vhodné hlavně v případech, kdy
lze reálně popsat rozložení délkového měřítka l , nejsou však vhodné pro modelování
složitějších případů proudění, kdy nelze s dostatečnou přesností definovat jeho rozložení
pomocí empirického vztahu. Zde je nutno definovat další transportní rovnici a přejít na
dvourovnicový model turbulence.
5.3.
Dvourovnicový k-e model
Dvourovnicový k-e model určuje turbulentní viskozitu pomocí dvou transportních rovnic
pro k a e . Model využívá Boussinesqovy hypotézy o vírové viskozitě a vztahuje mt ke k ,
e a Cm
m t = Cn
k2
e
( 5.3.1)
Rozložení k je dáno transportní rovnicí. Exaktní tvar transportní rovnice pro e je možné
odvodit opět z Navierových - Stokesových rovnic, tato rovnice obsahuje komplexní korelace,
které je opět nutné aproximovat. Výsledný tvar rovnice pro rychlost disipace používaný v tzv.
modelu k - e , je uváděn v této podobě:
¶ k ¶ ujk
¶
+
=
¶t
¶ xj
¶ xj
3/ 2
æ nt ¶ k ö
æ¶u
ö
ç .
÷ + n t ç j + ¶ u l ÷ ¶ u l - CD k
çs ¶ x ÷
ç¶ x
¶ x ÷ø ¶ x j 1
l3
j ø
è k
1è44l 424j4
43 424
P
47
e
( 5.3.2)
Statistické modely turbulence
¶ e ¶ u je
¶
+
=
¶t ¶ xj
¶ xj
5.4.
æ nt ¶e
ç .
çs ¶ x
j
è e
ö
æ ¶u
÷ + C1en t ç j + ¶u l
÷
ç ¶x
ø
è l ¶x j
ö ¶u l
e2
÷
- C 2e
÷ ¶x
k
ø j
( 5.3.3)
RNG k-e model
Tento model je odvozen z klasického k - e modelu při využití matematického
postupu nazvaného metoda renormalizačních grup (RNG). Renomalizační procedura
aplikovaná na turbulenci spočívá v postupné eliminaci malých vírů, přitom se přetransformují
pohybové rovnice (Navierovy - Stokesovy rovnice) tak, že se modifikuje turbulentní viskozita,
síly a nelineární členy. Předpokládá-li se, že tyto víry souvisí s disipací e , pak turbulentní
viskozita m t resp. n t =
mt
je závislá na měřítku turbulentních vírů a RNG metoda konstruuje
r
tuto viskozitu pomocí iteračního odstraňování úzkých pásem vlnových čísel. Podrobně je tato
metoda uvedena v [13] s tím, že pro iterační proces se využívá relace
dm eff
A el 3
= 1 2,
dl
m (l )
kde
( 5.4.1)
m efektivní = m turbulentn í + m molekulová
Integrací této rovnice přes délkové měřítko l pro počáteční podmínku n = n mol a pro měřítko
l = l d = L / Re 3 / 4 , což je Kolmogorovo disipační měřítko odpovídající malým turbulentním
vírům, se dostane
é
3 A1e 4 4 ù
(l - l d )ú
m eff (l ) = m mol 3 ê1 +
3
ë 4 m mol
û
(l ³ ld )
( 5.4.2)
Tato rovnice je interpolačním vzorcem pro výpočet m eff (l ) mezi molekulovou viskozitou a
viskozitou disipačních vírů s limitou lññ l d odpovídající vysokým Reynoldsovým číslům. Pro
vysoké Reynoldsovo číslo se dá dokázat, že rovnice ( 5.4.2) má tvar
m eff » m t = (0.094l ) Ñu
( 5.4.3)
2
Metoda RNG je asi o jednu desetinu pomalejší než klasický dvourovnicový model, ale
v oblastech zavíření (kde se kapalina zpomalí a je tam nižší Re číslo) je přesnější.
5.5.
Reynoldsův napěťový model (RSM)
Reynoldsův napěťový model zahrnuje výpočet jednotlivých Reynoldsových napětí
prostřednictvím diferenciálních transportních rovnic ve tvaru (viz ( 4.2.15))
48
Statistické modely turbulence
¶ u ¢i u ¢j
¶ u ¢i u ¢j
p¢
¶ é
¶
+ uk
=(u¢i u¢j )ùú ê(u ¢i u ¢j u ¢k ) + (d kj u i¢ + d ik u ¢j ) - n
¶t
¶x k
¶x
r
¶x k
14k4ë444444424444444
443û
difúzní transport
¶u j
é
é ¶u ¢ ¶u ¢j ù
¶u ù p¢ é ¶u ¢ ¶u ¢ ù
- êu i¢u ¢k
+ u ¢j u ¢k i ú + ê i + j ú - 2n ê i
¶x
¶x k û r ëê ¶x j ¶x i ûú
¶x k ¶x k úû
ë4
1ë4444k24444
3 1442443 14
244
3
produkce
redistribuce
( 5.5.1)
disipace
Za učelem uzavření soustavy rovnic Fluent aproximuje některé členy z rovnice ( 5.5.1), viz
[13]. Vypočtená Reynoldsova napětí jsou pak dosazována do rovnice pro přenos hybnosti (
4.2.12). Fluent tedy řeší:
·
šest transportních rovnic ve tvaru ( 5.5.1) pro Reynoldsova napětí
·
tři transportní rovnice pro středované složky rychlosti a rovnici kontinuity
·
transportní rovnici pro disipaci e
·
transportní rovnici pro turbulentní energii k v blízkosti stěny v případě explicitní
okrajové podmínky pro napětí jako součet rovnic ( 5.5.1)
·
5.6.
Modelování proudění v blízkosti stěny, stěnové funkce
Modelování proudění u stěny ovlivňuje přesnost numerického řešení v celé oblasti.
V blízkosti stěny se řešené veličiny rychle mění, výrazně se zde uplatňuje přenos hybnosti a
skalárních veličin. Turbulence je těsně u stěny potlačena, ve vnější části mezní vrstvy však
dochází k výrazné produkci turbulentní kinetické energie v důsledku Reynoldsových napětí a
gradientu střední rychlosti. Četné experimenty prokázaly, že oblast u stěny, tzv. mezní
vrstva, může být rozdělena na více části. Bezprostředně u stěny se nachází viskózní
(laminární) podvrstva, proudění je zde téměř laminární a molekulární viskozita má
dominantní vliv na přenos hybnosti, tepla a hmotnosti. Vnější část mezní vrstvy se označuje
jako plně turbulentní vrstva, dominantní úlohu zde hraje turbulence. Mezi laminární
podvrstvou a plně turbulentní vrstvou se vyskytuje přechodová vrstva, kde se stejnou
měrou uplatňují účinky molekulární viskozity i turbulence. Rozdělení mezní vrstvy je
znázorněno na obr. 5.2.
49
Statistické modely turbulence
8
8
7
u-exp
u-lam
u-turb
7
6
6
experiment
5
4
u [m/s]
5
3
vnější vrstva
2
4
plně vyvinutá turbulence
1
0
0.0
3
y [m]
0.5
u=u*(y+yo)/y
)
u=(u*/к). ln((y+yo)/y
)
1.0
2
1
viskózní podvrstva
přechodová vrstva
vyvinutá
0.010
0
0.001
plně vyvinutá turbulence
vnější vrstva
0.100
1.000
y [m]
obr. 5.2 Rozdělení oblasti v blízkosti stěny – v lineárních a logaritmických souřadnicích [26]
5.6.1. Teorie stěnových funkcí dle Laundera a Spaldinga
Stěnové funkce vycházející z teorie Laundera a Spaldinga jsou široce používané
hlavně pro průmyslové aplikace. V turbulentním proudění se mezní vrstva skládá z viskózní
podvrstvy a tzv. oblasti logaritmického zákona pro středovanou rychlost v turbulentní oblasti
ve zjednodušeném dvourozměrném případě:
+
u =
(
1
ln E.y +
k
)
( 5.6.1)
Bezrozměrné veličiny v této rovnici jsou definovány takto:
+
u =
u
ut
y+ =
rut y
m
ut =
tw
r
( 5.6.2)
kde
k
= von Kármánova konstanta (=0.42)
E
= empir. konstanta (=9.81)
u
= střední rychlost proudění v bodě P
ut
= třecí rychost
y
= vzdálenost bodu P od stěny ve směru normály
m
= dynamická viskozita tekutiny
Třecí rychlost ut
je určena pomocí smykového napětí definovaného podobně jako
Reynoldsovo napětí (veličiny mají index s ).
50
Statistické modely turbulence
ut2 º
tw
2
2
= u ¢v s¢ + w ¢v s¢
r
( 5.6.3)
kde totální vertikální toky hybnosti měřené v blízkosti stěny jsou určeny vztahy
2
2
t xy = - r u ¢v s¢ a t zy = - r w ¢v s¢ Þ t w = t xy
+ t zy
Reynoldsova napětí se vypočítají užitím Boussinesqovy aproximace. Pro k - e turbulentní
modely v případě lokální rovnováhy produkce turbulentní energie a rychlosti disipace
definované vztahem
- r u ¢v ¢
¶u
¶u
= tw
= re
¶y
¶y
( 5.6.4)
lze odvodit v logaritmické vrstvě vztah pro turbulentní kinetickou energii
kp =
ut2
Cm
( 5.6.5)
a disipaci
ep =
ut3
ky
( 5.6.6)
V rovnicích ( 5.6.1) až ( 5.6.6) vystupuje třecí (smyková) rychlost jako konstantní
veličina. Tyto standardní stěnové funkce jsou také definovány ve Fluentu a jsou vyhovující
pro celou řadu aplikací. Jestliže budou výše definované funkce použity pro obecnější
proudění, může právě smyková rychlost zapříčinit problémy. Navíc pro vyčíslení turbulentní
kinetické energie a disipace se předpokládá lokální rovnováha, což je velmi silné omezení
v proudění komplexního charakteru.
5.6.2. Modelování proudění v blízkosti stěny ve Fluentu
MODELOVÁNÍ
PROUDĚNÍ TEKUTINY
Stěnová funkce
Dvouvrstvý model
Standardní
Nerovnovážná
stěnová funkce
stěnová
51
Statistické modely turbulence
Proudění v blízkosti stěny lze modelovat dvěma způsoby:
·
užití stěnové funkce („wall function“), pomocí níž se překlene oblast laminární
podvrstvy a přechodové vrstvy, kde se uplatňuje molekulární i turbulentní viskozita, tj.
oblast mezi stěnou a oblastí plně vyvinutého turbulentního proudění.
·
modelování proudění u stěny („near-wall modelling“) včetně vazké podvrstvy
v souvislosti s jemností sítě. Podstata obou přístupů je znázorněna na obr. 5.3.
plně vyvinutá turbulence
proud tekutiny
P
Dy
viskózní +přechodová vrstva
P
stěna
užití logaritmické
stěnové funkce
metoda modelování
v blízkosti stěny
obr. 5.3 Přístup k modelování proudění u stěny ve Fluentu
Stěnové funkce představují soubor poloempirických vztahů a funkcí, které umožňují
„propojení“ řešené proměnné v buňce v blízkosti stěny s korespondující hodnotou na stěně.
Stěnové funkce zahrnují logaritmický zákon pro střední rychlost a teplotu a vztahy pro
turbulentní veličiny v blízkosti stěny.
FLUENT nabízí dvě základní varianty stěnové funkce:
·
standardní stěnové funkce
·
nerovnovážné stěnové funkce
V případě proudění s velkým Reynoldsovým číslem uplatnění stěnových funkcí podstatně
snižuje nároky na výpočet a poskytuje tak ekonomické a přitom dostatečně přesné řešení
pro většinu inženýrských problémů. Tento přístup je však nevhodný v případě proudění
s nízkým Reynoldsovým číslem. V těchto případech je nutné zvolit druhý přístup, který
umožňuje modelování v bezprostřední blízkosti stěny ovlivněné viskozitou proudící tekutiny a
52
Statistické modely turbulence
ve FLUENTU je definován jako tzv.
·
metoda modelování v blízkosti stěny - dvouvrstvý model
Při využití konečných objemů byly odvozeny obecnější stěnové funkce vhodné pro
komplexní proudění a jsou dané vztahy
u* =
(
1
ln E.y *
k
)
( 5.6.7)
kde
*
u =
uP k P1/ 2C m1/ 4 r
tw
,
y* =
rk P1/ 2C m1/ 4 Dy P
m
( 5.6.8)
Tato logaritmická formule pro střední rychlost je platná pro y * ñ 30 ¸ 60 . Ve Fluentu je
logaritmický předpis doporučován pro y * ñ11.225 . Když je hustota sítě u stěny taková, že
y * á11.225 v buňkách přiléhajících ke stěně, pak se aplikuje laminární vztah, tedy
(
ì1
*
ï ln E .y
u = ík
ïy *
î
*
)
( 5.6.9)
kde u * a y* jsou definovány v ( 5.6.8). Bezrozměrná tloušťka podvrstvy yv* je jednoduše
určena jako průsečík lineárního a logaritmického profilu a přitom přibližně platí
y v* =
rk p1 / 2Cm1/ 4 y v
» 11.23
m
Reynoldsova
( 5.6.10)
analogie mezi přenosem momentu a energie umožňuje definovat
podobný logaritmický zákon pro střední teplotu. Zákon stěny pro teplotu zahrnuje dvě
odlišné závislosti:
·
lineární závislost pro teplotně vodivou podvrstvu, kde vodivost má rozhodující
význam
·
logaritmický zákon pro turbulentní oblast, kde účinky turbulence převládají nad
vodivostí
Tloušťka teplotně vodivé podvrstvy je odlišná od tloušťky vazké podvrstvy a závisí na druhu
tekutiny. Zákon stěny pro střední teplotu je definován takto:
T
(T
=
w
(
ìPr .y y á y T
- Tp )rc pC m1/ 4 k p1/ 2 ï
=í
æ1
q ¢¢
ïPrt .ç k ln Ey
è
î
(
kde P je počítáno pomocí vztahu
53
)
) + P ö÷(y ñ y )
ø
T
( 5.6.11)
Statistické modely turbulence
1/ 2
P=
p /4 æ Aö
ç ÷
sin(p / 4 ) è k ø
æ Pr
öæ Pr ö
çç
- 1÷÷ç t ÷
è Prt
øè Pr ø
1/ 4
( 5.6.12)
a kde
k
= von Karmanova konstanta(=0.42)
E
= empir. konstanta stěny (=9.793)
Tp
= teplota v buňce přilehlé ke stěně
Tw
= teplota na stěně
Pr
= molekulové Prandtlovo číslo m c p / l
Prt
= turbulentní Prandtlovo číslo (0.85 na stěně)
A
= Van Driestova konstanta (=26)
´´
q
= teplotní tok
Bezrozměrná tloušťka tepelné podvrstvy yT v rovnici ( 5.6.11) se spočítá jako hodnota y, ve
které se protíná lineární a logaritmická závislost, při daném molekulovém Prandtlově čísle.
Postup při uplatnění zákona stěny pro teplotu je takový, že Prandtlovo molekulové číslo se
spočítá ze zadaných fyzikálních vlastností tekutiny, pak se počítá tloušťka tepelné podvrstvy
yT jako průsečík lineárního a logaritmického profilu a tato hodnota se ukládá. Během
iterování je aplikován buď lineární nebo logaritmický profil v závislosti na hodnotě y pro
výpočet teploty na stěně Tw nebo teplotního toku q´´ (viz rovnice pro ( 5.6.11)).
5.6.3. Nerovnovážná stěnová funkce
Tam, kde je proudění u stěny vystaveno účinkům velkého tlakového gradientu a kde
se nedá předpokládat splnění podmínky lokální rovnováhy, je možné dosáhnout přesnějších
výsledků pomocí nerovnovážné stěnové funkce. Podstata nerovnovážné stěnové funkce
spočívá v následujících krocích:
·
logaritmický zákon pro střední rychlost podle Laundera a Spaldinga je upřesňován
v závislosti na účinku tlakového gradientu
·
bilance turbulentní kinetické energie a disipace v buňce sousedící se stěnou je
počítána ve dvou vrstvách, tj laminární a turbulentní.
Stěnová funkce pro střední teplotu se neliší od standartní stěnové funkce pro teplotu.
Logaritmický profil střední rychlosti je přepočítán pomocí tlakového gradientu podle vztahu
u~ 1
= ln(Ey p )
ut k
æy
1 dp é y v
lnçç
u~ = u ê
1/ 2
2 dx ë rkk
è yv
( 5.6.13)
ö y - y v y v2 ù
÷÷ +
+ ú
1/ 2
mû
ø rkk
54
( 5.6.14)
Statistické modely turbulence
a yv je fyzikální tloušťka vazké podvrstvy definovaná výrazem
yv =
my *
rCm1/ 4 k p1/ 2
( 5.6.15)
kde y*= 11.225.
5.6.4. Použití stěnových funkcí a omezení
Standardní stěnové funkce
umožňují dostatečně
přesné
řešení pro
velká
Reynoldsova čísla. Nerovnovážná stěnová funkce dále rozšiřuje možnost aplikace stěnové
funkce na případy, kdy je proudění vystaveno účinkům tlakového gradientu a nerovnováhy.
Tento přístup není vhodný, jestliže se proudění příliš liší od ideálních předpokladů,
na nichž je metoda stěnové funkce založena. Jedná se například o tyto případy:
·
proudění s nízkým Reynoldsovým číslem, velký vliv stěny (např. proudění úzkou
štěrbinou, proudění velmi vazkých tekutin, proudění s malou rychlostí)
·
silný tlakový gradient vedoucí k odtržení mezní vrstvy
·
významné působení objemových sil (např. odstředivé síly, proudění v blízkosti
rotujícího disku, Archimedovy síly)
·
trojrozměrné proudění v blízkosti stěny (Ekmanova spirála, silně zakřivená 3D mezní
vrstva).
Odpovídá-li charakter proudění některému z výše uvedených případů a je-li nutné tyto jevy
zahrnout do řešení problému, pak musíme přistoupit k podrobnému modelování proudění u
stěny (near-wall modelling). K tomuto účelu Fluent nabízí dvouvrstvový model.
5.6.5. Dvouvrstvový model (near-wall modelling)
V tomto modelu je celá oblast rozdělena na část, ve které se projevuje vliv viskozity, a
na plně turbulentní oblast. Hranice mezi oběma oblastmi je definována pomocí turbulentního
Reynoldsova čísla
Re y º
r kn
m
( 5.6.16)
kde n je normálová vzdálenost středu buňky od stěny, ve FLUENTU je n interpretováno
jako vzdálenost od nejbližší stěny:
r r
n º min
r - rw
r
(5.6.17)
rw e Gw
r
r
kde r je polohový vektor v bodě pole a rw je polohový vektor na stěně. Gw je sjednocení
všech stěn na hranici. Tato interpretace umožňuje jednotnou definici n
komplexního tvaru, která zahrnuje mnoho stěn.
55
v oblasti
Statistické modely turbulence
V plně turbulentní oblasti (Rey>200) se používají turbulentní modely popsané v kapitole 4.1 4.3. V oblasti ovlivněné viskozitou v blízkosti stěny (Rey< 200), je aplikován jednorovnicový
model podle Wolfsteina, kde rovnice pro přenos hybnosti a rovnice pro k je definována
známým způsobem, turbulentní viskozita je počítána ze vztahu
mt = rC m k 1 / 2l m
(5.6.18)
a rychlost disipace
e=
k3/2
le
(5.6.19)
Délková měřítka lm a le , která se v těchto rovnicích vyskytují, jsou definována následovně:
é
æ Re y
l m = Cl y ê1 - expçç è Am
ëê
öù
÷ú
÷
øûú
(5.6.20)
é
æ Re y
l e = Cl y ê1 - expçç è Ae
ë
öù
÷÷ú
øû
(5.6.21)
Jestliže v celé oblasti je Rey <200, pak e není řešeno pomocí transportní rovnice, ale
z algebraického vztahu (4.4.22 ). Konstanty Am, Ae , vyskytující se v těchto vztazích, jsou
definovány podle autorů Chena a Patela:
(5.6.22)
Am = 70, Ae =2 cl
Jestliže se používá RSM model, jako stěnové funkce mohou být využity standardní a
nerovnovážné funkce, dvouvrstvový model stěnových funkcí není vhodný.
5.6.6. Vliv kvality sítě na volbu stěnové funkce pro různé modely
turbulence
Vzdálenost středů buněk sousedících se stěnou od těchto stěn je určující pro to, je-li
správný přístup logaritmické stěnové funkce nebo je třeba volit jiný.
·
logaritmický předpis je platný pro y * ñ 30 ¸ 60
·
dvouvrstvý předpis je platný pro y * á 4 ¸ 5 , v ideálním případě nejméně 10 buněk má
být v laminární podvrstvě
·
Spalart-Allmaras model využívá logaritmické stěnové funkce za předpokladu velmi
jemné sítě ( y * = 1 ) nebo sítě, pro kterou je y * ³ 30 .
·
Large Eddy Simulation model používá logaritmickou stěnovou funkci pro velmi
jemnou síť (řádu y * = 1 )
56
Statistické modely turbulence
5.6.7. Vliv drsnosti na stěnovou funkci
Drsnost stěn se zohledňuje ve vztahu pro logaritmickou stěnovou funkci následovně:
up =
u* æ r u* Dy ö
lnç E
÷ - u * DB
k çè
m ÷ø
( 5.6.23)
kde u p je střední hodnota rychlosti v bodě P nejblíže stěny a Dy je vzdálenost bodu P od
stěny.
DB =
Funkce
DB
závisí na
typu
a
velikosti
drsnosti a
je určena
vztahem
rK s u *
1
ln 1 + CKs K s+ , kde K s+ =
je bezrozměrná drsnost a K s skutečná fyzická
k
m
(
)
drsnost v metrech, CKs = 0.5 pro pravidelnou drsnost, vyšší hodnoty (CKs = 0.5 ¸ 1.0) jsou
doporučovány pro nepravidelnou drsnost.
Drsností lze nahradit pravidelné objekty vyskytující se v určité oblasti (soustava
trubek v chladiči, lesní porost při sledování proudění v atmosféře atd.), ale je třeba zvážit, jak
hustá bude síť u stěny, jestli drsnost nebude rozměrem převyšovat velikost buněk u stěny.
Pak by se musely překážky modelovat jako obtékané objekty.
Drsnost ovlivní tedy profil rychlosti v blízkosti stěny (tj. rychlost se sníží), ale i všechny
další veličiny, jako je profil teploty, koncentrace příměsi, apod.
100
U(z0=0.001)
U(z0=0.01)
U(z0=0.1)
U(z0=0.5)
U(z0=1)
U(z0=2)
U(moc)
90
80
70
z[m]
60
50
40
30
20
10
0
-4
-2
0
2
4
6
8
U [m/s]
obr. 5.4 Grafy rychlosti podle různých aerodynamických drsností při proudění v atmosféře
57
Matematický model turbulence pro stlačitelné neizotermní proudění
6. Matematický model turbulence pro stlačitelné
neizotermní proudění
Pro jednoduchost je uveden standardní k - e dvourovnicový model turbulence pro
proudění stlačitelné tekutiny s okrajovými podmínkami zahrnujícími turbulentní veličiny [15],
[11].
k-e dvourovnicový model turbulence
6.1.
Rovnice lze odvodit postupem uvedeným v kap.4.2.3 a 4.2.4, a jejich tvar je
následující:
rovnice kontinuity platná pro středované veličiny
¶ r ¶ (r u j )
+
=0
¶t
¶ xj
( 6.1.1)
rovnice pro přenos hybnosti
¶ ui
¶ (r u i ) ¶ (r u i u j )
¶p
¶ æç
+ mt )
+
=+
(
m
¶x j
¶t
¶ xj
¶ x i ¶ x j çè
ö
÷+
÷
ø
( 6.1.2)
+ r d i 3 g + r fc e ij 3 u j + r f i
123
1424
3
vztlakové síly
Coriolisov y síly
V případě dvourovnicového k-e modelu jsou tyto rovnice doplněny rovnicí pro přenos
kinetické turbulentní energie k a rychlosti disipace e.
¶ (rk ) ¶ (r u j k ) ¶
+
=
¶t
¶ xj
¶ xj
æ mt ¶ k
ç .
çs ¶ x
j
è k
ö
ö
æ
÷ + mt ç ¶ u j + ¶ u l ÷ ¶ u l - g j mt ¶r - re
÷
ç ¶ x ¶ x ÷¶ x
rs h ¶x j
j ø
ø 14
è 4l4244
4
3
43j 14424
P
( 6.1.3)
G
¶ (re ) ¶ (r u j e )
¶ æç m t ¶ e ö÷
e2
+
=
.
+ rc1e (P + c 3e G ) - rc 2e
¶t
¶ xj
¶ x j çè s e ¶ x j ÷ø
k
( 6.1.4)
kde P a G představují produkci turbulentní kinetické energie v důsledku napětí a vztlakových
sil
é¶ u j ¶ u l ù ¶ u j
m
¶r
P = mt ê
+
, G = -g j × t ×
ú
rs h ¶ x j
êë ¶ x l ¶ x j úû ¶ x l
( 6.1.5)
kde C1e=1,44 , C2e=1,92, C3e=1, s k = 1, s e =1.3 jsou konstanty určené empiricky a
sh =
mt
c p je Prandtlovo turbulentní číslo.
lt
Pro doplnění jsou Reynoldsova napětí u¢i u¢j definována vztahem
- r u i¢u ¢j = m t
¶ ui
¶ xj
( 6.1.6)
58
Matematický model turbulence pro stlačitelné neizotermní proudění
kde turbulentní viskozita mt se předpokládá jako funkce délkového a rychlostního měřítka
dle Kolmogorov-Prandtlovy hypotézy:
mt » l .v = rC m
k2
e
( 6.1.7)
Okrajové podmínky pro k-e turbulentní model
6.2.
6.2.1. Hmotnostní průtok
Rychlostní podmínka se používá k definování okrajové podmínky na průtočné hranici
do oblasti při proudění nestlačitelného proudění. U stlačitelného proudění se předpokládá
nekonstantní hustota, která je závislá na stavových veličinách tlaku a teplotě a ovlivňuje
objemový průtok a tím rychlost, což může vést k nereálným výsledkům. V tomto případě se
zadává hmotnostní průtok Qm = rQ = ruS .
6.2.2. Turbulentní veličiny
Velký
význam
v souvislosti
se vstupní
okrajovou
podmínkou
má
nastavení
turbulentních parametrů v podobě hodnot turbulentní kinetické energie a rychlosti
disipace. Přesnější je samozřejmě vyjádření těchto veličin profilem získaným z empirických
dat nebo z empirických formulí. Pokud není profil přesně znám, lze zadat konstantní hodnotu
odhadnutou na základě zkušenosti. Tyto turbulentní veličiny mohou být určeny případně
pomocí veličin snadněji určitelných jako je intenzita turbulence, poměr turbulentní a
molekulové viskozity, hydraulického průměru a délkového měřítka turbulence. Velikost
turbulentních fluktuací se obvykle popisuje intenzitou turbulence. Za předpokladu izotropní
turbulence ( u1/ 2 =
u 2/ 2 =
u3/ 2 ) se vyjadřuje relativní intenzita turbulence jako poměr
efektivní hodnoty fluktuační složky rychlosti ke střední rychlosti ve stejném místě
proudu obvykle vyjádřený v procentech. Zpravidla se měří pouze jedna směrová složka:
I=
u1/ 2
u1
( 6.2.1)
Běžné turbulentní proudění je anizotropní (nesourodé v souřadnicových směrech), ale
anizotropie bývá malá. Největší rozdíly jsou mezi podélnou a příčnou složkou pohybů.
Obecně
/
I=
/
uj uj
,
3uu
( 6.2.2)
Rozdíl mezi fluktuacemi rychlostí v příčném směru u2/ a u3/ je zpravidla velmi malý. Hodnota
intenzity je přibližně:
59
Matematický model turbulence pro stlačitelné neizotermní proudění
I [%]
aerodynamický tunel
0.05%
turbulentní proudění generované mříží
1-5%
úplavy
2-10%
proudění v mezní vrstvě a při průtoku
5-20%
trubicí
zatopený proud
20%
recirkulační proudění s malou rychlostí u
100%
Turbulentní měřítko l je omezeno velikostí oblasti, protože turbulentní víry nemohou
být větší než je rozměr oblasti. Přibližná hodnota turbulentního měřítka se určí ze vztahu
l = 0.07L kde L je charakteristický rozměr oblasti, případně hydraulický průměr. Intenzita
turbulence a hydraulický průměr jsou dostupné veličiny, které je možno zadat jako okrajové
podmínky, ostatní se pak přepočítají dle následujících vztahů
intenzita turbulence
uj uj
I=
3uu
turbulentní měřítko
l = 0.07L
poměr turbulentní viskozity
/
mt
m
n =
turbulentní kinetická energie
k=
rychlost disipace
/
3
u.I.l
2
3
(u¢)2 nebo k = 3 (uI )2
2
2
3
4
m
3
2
k
k 2 æ mt
= rCm
ç
e =C
l
m çè m
ö
÷÷
ø
-1
Samozřejmě turbulentní energii a rychlost disipace lze definovat také přímo. Podle složitosti
matematického modelu se definují také další veličiny související s přenosem tepla případně
další skalární veličiny. Hodnota turbulentní intenzity v případě LES se definuje pomocí
náhodné fluktuace rychlostního pole na vstupu.
6.2.3. Tlak na vstupu
Tlaková podmínka na vstupu se používá, pokud je znám celkový (totální) tlak nebo
statický tlak a průtok nebo rychlost. Tato podmínka je vhodná i pro proudění s uvažováním
vztlakových sil.
Na vstupu se definuje celkový (totální) relativní tlak (vztažený k operačnímu tlaku)
vztahem odvozeným z Bernoulliho rovnice, přitom hustota je konstantní nebo je funkcí
teploty resp. hmotnostních zlomků příměsí:
60
Matematický model turbulence pro stlačitelné neizotermní proudění
ptot = pstat +
1
2
ru
2
( 6.2.3)
Pro stlačitelné proudění
k
ptot
kde
é k - 1 2 ù (k -1)
Ma ú
= pstat ê1 +
2
ë
û
( 6.2.4)
ptot
celkový (totální) tlak
pstat
statický tlak
Ma
Machovo číslo Ma =
r
měrná plynová konstanta r =
k
poměr měrných tepel k = c p
c
rychlost zvuku v tekutině
u
=
c
u
(krTs )0.5
R
, M je molekulová váha
M
cV
Rozdíl při vyhodnocení celkového tlaku při proudění vzduchu do Machova čísla Ma = 0.1 je
patrný z obr. 6.1
1200
0.14
0.12
p dyn
p stat
p celk
p tot stlac
Ma
p [Pa]
800
0.1
0.08
600
0.06
400
Machovo číslo [1]
1000
0.04
200
0.02
0
0
0
5
10
15
20
25
30
35
40
u [ms-1]
obr. 6.1 Zvislost tlaků na rychlosti (vzduch)
V případě, že je proudění ovlivněno vztlakovými silami, je tlakové pole zvětšeno o
hydrostatický tlak:
61
Matematický model turbulence pro stlačitelné neizotermní proudění
p¢ = r ref gx i + p
( 6.2.5)
Tedy počítá se odchylka od hydrostatického tlaku a také se vyhodnocuje. Je ale nutné zadat
rozumnou hodnotu referenční hustoty rref při referenční teplotě. Vstupní hodnoty tlaku
celkového i statického jsou o hydrostatický tlak zvětšeny.
Při zadání tlakové podmínky je nutné definovat směr proudění pomocí složek
rychlosti případně pomocí proudění v normálovém směru k hranici.
Statický tlak na vstupu musí být specifikován v případě supersonického proudění.
Turbulentní veličiny jsou určovány shodně jako v případě okrajové podmínky hmotnostního
průtoku.
6.2.4. Tlak na výstupu
Tlaková okrajová podmínka na výstupu se zadává v podobě statického tlaku. Statický
tlak se definuje jen v případě subsonického proudění. Pokud je proudění supersonické, tak
se tlak i ostatní veličiny extrapolují z proudění uvnitř oblasti. Pokud se objevuje během
výpočtu zpětné proudění, je tato podmínka vhodnější než outflow, protože dosahuje lepší
konvergence. Pro zpětné proudění je ale nutné určit reálné okrajové podmínky ostatních
počítaných veličin, což je teplota a turbulentní veličiny, případně další skalární veličiny.
6.2.5. Outflow
Outflow je nevhodné pro stačitelné proudění, nestlačitelné nestacionární proudění
s měnící se hustotou a v případě zadaného tlaku na vstupu.
62
Řešení přenosu tepla (konvekce, kondukce)
7. Řešení přenosu tepla (konvekce, kondukce)
7.1.
Úvod do problematiky modelování přenosu tepla
Teplo, jakožto forma energie, může být přenášeno 3 základními způsoby [23]:
·
kondukcí (vedením),
·
konvekcí (prouděním)
·
radiací (sáláním).
Jedním ze základních pojmů z oblasti sdílení tepla je pojem „teplotní pole“. Jedná se
o druh pole, které udává rozložení teplot v určitém časovém okamžiku ve všech bodech
sledovaného prostoru. Mezi další základní pojmy patří například „teplotní gradient“, „tepelný
tok“ nebo „hustota tepelného toku“.
Teplotní gradient ( grad T ) představuje vektorovou veličinu, která udává změnu
teploty ve směru normály k izotermickému povrchu. Lze jej matematicky vyjádřit vztahem:
grad T = lim
n ®0
D T ¶T
=
Dn ¶n
.
( 7.1.1 )
kde DT představuje „rozdíl relativních či absolutních teplot“ mezi zkoumanými izotermickými
plochami a Dn je „rozdílem vzdáleností izotermických vrstev“ ve směru normály k těmto
vrstvám (plochám).
Tepelný tok – tepelný výkon Q& představuje množství tepelné energie přenesené za
jednotku času, a je definován vztahem:
dQ
Q& =
dt
( 7.1.2 )
Hustota tepelného toku q vyjadřuje tepelný tok vztažený k jednotce plochy, a je
definována vztahem:
q=
dQ&
dS
( 7.1.3 )
7.1.1. Kondukce – vedení tepla
Sdílení tepla kondukcí (vedením) vzniká v důsledku pohybu strukturních částic hmoty
(molekul). V čisté formě existuje pouze u látek pevného skupenství. V tekutinách existuje
pouze za předpokladu, že je vliv pohybu tekutiny na přenos tepla zanedbatelný (tj. že se
nejedná o přenos tepla konvekcí). Ve všech ostatních případech přispívá k přenosu tepla i
konvekce nebo dokonce konvekce spolu s radiací. Při většině inženýrských úloh do
fyzikálního děje vstupují všechny tři způsoby sdílení tepla.
Nejčastěji je tento fyzikální děj demonstrován na jednoduchém příkladu vedení tepla
homogenní nekonečnou stěnou o určité tloušťce tvořenou konkrétním materiálem. V oblasti
63
Řešení přenosu tepla (konvekce, kondukce)
sdílení tepla kondukcí (vedením) se ještě rozlišují dva dílčí pojmy. Prvním z nich je pojem
„vedení tepla“ stěnou, který vyjadřuje situaci, kdy jsou pro výpočet zadány obě povrchové
teploty stěny, viz obr. 7.1. Druhým pojmem je pak „prostup tepla“ stěnou, o kterém
hovoříme, pokud jsou pro výpočet zadány teploty obou médií (tekutin), jež stěnu z obou stran
obklopují a které se mohou také pohybovat.
VEDENÍ TEPLA
STĚNOU
Ts1
l
Ts 2
Dx
obr. 7.1 Princip „vedení tepla“
Základním zákonem vedení tepla je Fourierův zákon, který udává vztah mezi
hustotou tepelného toku q a teplotním gradientem grad T :
q = -l grad T
.
( 7.1.4 )
kde l je tepelná vodivost, která závisí na druhu materiálu a mění se s teplotou. Záporné
znaménko na pravé straně rovnice vyjadřuje skutečnost, že hustota tepelného toku a teplotní
gradient mají jako vektory opačný smysl (teplo se šíří ve směru klesající teploty).
Nejjednodušším případem vedení tepla je stacionární vedení tepla homogenní
neomezenou izotropní rovinnou stěnou (viz obr. 7.1), tj. stěnou, která je materiálově
stejnorodá a teplo se v ní šíří ve všech směrech stejným způsobem. Pro výpočet hustoty
tepelného toku v tomto případě platí základní vztah:
q = -l ×
TS 1 - TS 2
Dx
.
( 7.1.5 )
kde TS 1 a TS 2 představuje teploty na stěnách a Dx je tloušťkou stěny. Obdobný princip lze
aplikovat u tzv. stěny složené, která se skládá z několika vrstev o různém materiálu. Tento
postup se provede pro každou dílčí vrstvu zvlášť.
Nejjednodušším případem prostupu tepla je stacionární prostup tepla homogenní
neomezenou izotropní rovinnou stěnou (viz obr. 7.1). Podmínkou však je, aby se tekutina
64
Řešení přenosu tepla (konvekce, kondukce)
obklopující stěnu z obou stran (viz výše uvedený výklad) se výrazněji nepohybovala a
nedocházelo tak ke sdílení tepla prouděním. Pro výpočet hustoty tepelného toku v tomto
případě platí základní vztah:
q=
1
× (T2 - T1 )
1 Dx 1
+
+
a1 l a 2
.
( 7.1.6 )
kde a1 a a 2 představuje součinitele přestupu tepla na rozhraní stěn a tekutiny, T1 a T2
představují teploty obou stěn obklopujících tekutin a Dx je tloušťkou stěny. Na rozdíl od
případu vedení tepla nelze stejný postup uplatnit u stěn složených.
7.1.2. Konvekce
Sdílení tepla konvekcí (prouděním) je na rozdíl od kondukce spjato s pohybem tekutiny
okolo nebo podél stěny, jíž teplo prochází. Opět stejně jako v případě kondukce, neexistuje
v praxi nikdy ve zcela čisté podobě, protože uvnitř proudící tekutiny stejně jako na rozhraní
tekutiny a pevné stěny je vždy doprovázeno kondukcí. Výpočet přenosu tepla je však oproti
případu kondukce ztížen faktem, že je třeba do něj zakomponovat kromě rovnic z oblasti
sdílení tepla ještě rovnice z oblasti hydrodynamiky. Tepelnou konvekci tedy v souhrnu
popisuje následující matematický model přenosu tepla.
Ts>Tvz
u
T
qK
Ts
obr. 7.2 Princip „konvekce tepla“ a profily rychlosti a teploty v blízkosti stěny obtékané
vzduchem [27]
Jestliže je tekutina zahřívána a její hustota se mění v závislosti na teplotě, může
v důsledku rozdílu hustot a působení tíhového zrychlení dojít k proudění, které je
označováno jako smíšená (příp. přirozená) konvekce. Význam vztlakových sil při smíšené
konvekci lze posoudit pomocí poměru Grashofova a druhé mocniny Reynoldsova čísla:
Gr
Drgh
=
2
Re
ru 2
( 7.1.7)
65
Řešení přenosu tepla (konvekce, kondukce)
Gr
ññ1 - vztlakové síly významně ovlivňují proudění v oblasti.
Re 2
Gr
áá1 - vztlakové síly se zanedbávají.
Re 2
V případě pouze přirozené konvekce je možné určit podle hodnoty tzv. Rayleighova čísla
definovaného vztahem Ra =
a=
1 ¶r
gb DTl 3 r
, kde b = je součinitel teplotní roztažnosti a
ma
r ¶T
l
je součinitel teplotní vodivosti, zda proudění je laminární nebo turbulentní.
r.c p
Ra<108
laminární proudění
8
10 <Ra<10
7.2.
10
přechod do turbulence
Matematický model přenosu tepla
7.2.1. Přenos tepla v tekutinách
Výše uvedené principy kondikce a konvekce se v reálných problémech vyskytují
současně a řeší se v celém komplexu ([32], [33]), z čehož plne složitost matematického
modelu přenosu tepla v tekutinách, který vychází opět ze zákona zachování hmoty,
hybnosti a navíc energie. Rovnice energie je formálně podobná předchozím rovnicím a za
předpokladu turbulentního proudění (všechny veličiny jsou uvažovány turbulentní) je
definována takto:
¶ (t u )
¶
(r E ) + ¶ (r u j E ) = ¶ p + r u j f j + j l j + ¶
¶t
¶t
¶x j
¶ xl
¶ xj
kde
E =U +
æ
ö
ç leff ¶ T ÷
ç
¶ x j ÷ø
è
( 7.2.1)
1 2
u j celková energie, která je součtem vnitřní a kinetické energie
2
leff = l + lt , l je součinitel molekulové teplotní vodivosti, lt je součinitel trubulentní
teplotní vodivosti.
Zavede-li se pojem entalpie, definované vztahem
h =U +
pak E = h -
p
r
p 1 2
+ u j . Entalpie h je definovaná pro ideální plyny jako
r 2
T
h=
ò c dT
p
[Jkg-1]
Tref
a pro nestlačitelné médium (nestlačitelné plyny a kapaliny) jako
66
( 7.2.2)
Řešení přenosu tepla (konvekce, kondukce)
T
h=
òc
p
dT +
Tref
p
r
[Jkg-1]
( 7.2.3)
Ve výše uvedených rovnicích pro entalpie je výpočet definován pro referenční teplotu (např.
Tref = 298,15 K ), kterou lze měnit podle situace.
7.2.2. Přenos tepla ve vodivých stěnách
Přenos tepla ve vodivých stěnách je definován obdobnými veličinami jako
v tekutinách. Neuvažují se fluktuace teploty, neboť v pevných látkách nemají smysl. Rovnice
entalpie je:
rw
kde
¶hw
¶
=
¶t
¶x j
æ ¶T ö
÷ & //
ç lw
ç ¶x ÷ + q
j ø
è
( 7.2.4)
rw
hustota stěny,
hw
entalpie stěny, cw(T – Tref),
lw
tepelná vodivost stěny,
T
teplota stěny,
q& //
jednotkový zdroj tepla ve stěně, který je dán jako objemový vývin tepla disipací
[Wm-3] nebo např. generací tepla v elektronických komponentech
Pokud jsou řešeny úlohy, kde ještě dochází k pohybu či rotaci stěn, pak tyto efekty
jsou zahrnuty v řešení rovnice energie:
rw
¶hw
¶
+ rw
(v i,w hw ) = ¶ æçç lw ¶T ö÷÷ + q& ///
¶t
¶xi
¶xi è ¶xi ø
( 7.2.5)
Konvekce tepla je z důvodu pohybu stěny rychlostí v i ,w zahrnuta v rovnici energie pro oblasti
ohraničující proudění. Na pohybující se vodivé stěně je nutné zadat následující parametry:
·
v i,w rychlost pohybující se stěny
·
drsnost stěny (pouze při řešení turbulentního proudění)
·
specifické teplo (neustálené proudění, pohybující se stěny)
Zadání teplotní vodivosti umožňuje řešit úlohy, kde pevná vodivá oblast je tvořena
oddělenými stěnami z různých materiálů a různých vlastností. Hustota a specifické teplo
stěny jsou důležité při řešení časově závislých úloh a při řešení ustáleného stavu pouze
tehdy, když se stěna pohybuje. Typickými příklady jsou řešení dopravníkových pásů,
pohybujících se ocelových válcovaných pásů v pecích a úlohy s rotačními strojními
součástmi.
Všechny fyzikální vlastnosti mohou být podle charakteru úlohy konstantní nebo závislé
67
Řešení přenosu tepla (konvekce, kondukce)
na teplotě případně na tlaku. Nejvýznamnější veličinou je hustota.
7.3.
Hustota závislá na teplotě a tlaku, Boussinesqova
aproximace
Definování hustoty a případně vztlakových sil v matematickém modelu proudění závisí na
tom, zda proudění je modelováno jako stlačitelné nebo nestlačitelné
7.3.1. Vyjádření hustoty pro stlačitelné médium
Pro stlačitelné médium (ve Fluentu jen plyny) se využívá stavové rovnice pro ideální
plyn definující hustotu ve tvaru
r=
pop + p
rT
( 7.3.1)
7.3.2. Vyjádření hustoty pro nestlačitelné médium
Pro nestlačitelné médium (ve Fluentu plyny i kapaliny) se může vybrat jedna
z následujících metod:
r = konst .
·
konstantní hustota
·
stavová rovnice pro nestlačitelný plyn, pokud jsou změny tlaku velmi malé a plyn je
tedy nestlačitelný, pak
r=
pop
rT
Operační tlak pop musí být definován dosti přesně.
·
závislost hustoty na teplotě (u přirozené konvekce) je vyjádřena jako:
r n +1 - r n
(T - Tn ) pro Tn áT áTn +1
Tn +1 - Tn
po částech lineární funkce
r (T ) = r n +
polynomická funkce
r (T ) = A1 + A2T + A3T + ...
2
Na obr. 7.3 je vyjádřena závislost hustoty vzduchu pro výše uvedené varianty, přitom
polynomická funkce je volena pro jednoduchost jako přímka a je získána metodou
nejmenších čtverců
68
Řešení přenosu tepla (konvekce, kondukce)
1.6
1.5
-3
r[kg.m ]
1.4
1.3
1.2
ro konst
ro ideal (prel=0kPa)
ro ideal (prel=5kPa)
ro ideal (prel=10kPa)
ro ideal (prel=20kPa
ro(A1,A2)
1.1
1
0.9
0
5
10
15
20
25
30
35
40
45
50
t[0C]
obr. 7.3 Grafické vyjádření hustoty v závislosti na teplotě dle předchozích vztahů a pro různé
tlaky
7.3.3. Boussinesquova aproximace
Je-li proudění modelováno jako nestlačitelné, lze definovat závislost hustoty na
teplotě nebo použít Boussinesqův model, kde hustota vystupuje jako konstantní (referenční)
hodnota ve všech řešených rovnicích kromě vztlakového členu v pohybové rovnici,
¶ (r u i ) ¶ (r u i u j )
¶p
¶ æç ¶ u i
+
=+
m
¶t
¶ xj
¶ x i ¶ x j çè ¶ x j
ö
÷ + rd i 3 g + r fc e ij 3 u j + r fi
÷
123
1424
3
ø vztlakové síly Coriolisov
y síly
který se formálně nahradí změnou teploty takto:
(r - rref ) × g i
= - r ref b (T - Tref ) × g i
( 7.3.2)
kde r ref a Tref jsou referenční hustota a teplota a b je součinitel teplotní roztažnosti. Platí
tedy, že r = rref (1 - b DT ) . Tato aproximace je dostatečně přesná pokud jsou změny
hustoty malé.
*
Nechť p = rrefd i 3 g i xi + p (tj. pro i=1,2 je p=p*), kde první člen na pravé představuje
hydrostatický tlak definovaný pomocí referenční hustoty a p
*
je přírůstek tlaku, platí za
předpokladu hydrostatické rovnováhy
¶p
= rref × d i 3 g i
¶ xi
( 7.3.3)
69
Řešení přenosu tepla (konvekce, kondukce)
Po dosazení do vztlakového
členu pomocí rozdílu hustot
(r - r ref ) × g i
a nahrazením
hustoty r konstantní referenční hustotou rref se dostane upravená pohybová rovnice [13] pro
korekce tlaku p*:
æ
é
¶ u j ù æ 2 ¶ u l ö ö÷ ¶ p*
ç mt × ê ¶ u i +
÷÷ +
ú - çç m t
÷ ¶ xi
ç
x
x
x
3
¶
¶
¶
ê
ú
j
i
l
è
ø
ë
û
ø
è
+ d i 3 (r - r ref )g i + Fi
¶
(rref u i ) + ¶ (rref u i u j ) = ¶
¶t
¶ xj
¶ xj
( 7.3.4)
Tato rovnice se obdrží formálně následující úpravou tak, že k pravé straně se přičte a odečte
člen rrefgi, kde r ref g i =
¶ pref
, přičemž pref je hydrostatický tlak odpovídající referenční
¶ xi
hustotě rref . Pak pravá strana má tvar:
....... =
¶
¶ xj
æ
ç mt
ç
è
é ¶ u i ¶ u j ù æ 2 ¶ u l ö ö ¶ p ¶ pref ¶ pref
÷÷ ÷ ×ê
+
+
+ r g i + Fi
ú - çç mt
÷ ¶ xi
¶
¶
3
¶
¶ xi
¶ xi
x
x
x
ú è
i û
l øø
ëê j
..... =
¶
¶ xj
æ
ö
é
ù
ç mt × ê ¶ u i + ¶ u j ú - æç 2 m ¶ u l ö÷ ÷ - æç ¶ p - ¶ pref
t
÷÷ ç
ç
ç
¶ xi
ëê ¶ x j ¶ xi ûú è 3 ¶ xl ø ø è ¶ x i
è
ö
÷÷ + (r - r ref ) × g i + Fi
ø
..... =
¶
¶ xj
æ
ö
é
ù
ç mt × ê ¶ u i + ¶ u j ú - æç 2 m ¶ u l ö÷ ÷ - æç ¶ p - ¶ pref
t
ç
÷
ç
ç
¶ xi
êë ¶ x j ¶ xi úû è 3 ¶ xl ø ÷ø è ¶ x i
è
ö
÷÷ - r ref b (T - T ref ) × g i + Fi
ø
*
¶ p ¶ pref ¶ p
kde
.
=
¶ xi ¶ xi
¶ xi
1.4
1.3
-3
r[kg.m ]
1.2
ro (Tref=0)
ro (Tref=10)
ro (Tref=20)
ro (Tref=30)
ro (Tref=40)
ro (Tref=50)
ro (Tref=60)
ro ideal (prel=0kPa)
1.1
1
0.9
0
5
10
15
20
25
30
35
40
45
50
55
60
0
t[ C]
obr. 7.4 Závislost hustoty na referenční teplotě, srovnání s hustotou počítanou stavovou
rovnicí.
70
Řešení přenosu tepla (konvekce, kondukce)
Zanedbání změny hustoty v členu na levé straně rovnice (storage člen – viz [12]) a její
zachování ve vztlakovém členu na pravé straně rovnice (buyoancy člen – viz [12]) se nazývá
Boussinesqova aproximace. Na obr. 7.4 je porovnána závislost hustoty vzduchu pro různé
referenční teploty a hustoty vzduchu počítané z rovnice pro ideální plyn a pro atmosférický
tlak.
7.4.
Okrajové podmínky pro neizotermní proudění
Řešená oblast, sestávající z proudící tekutiny a případně vodivých stěn a okolního
prostředí, je omezena hranicemi, na kterých se definují okrajové podmínky.
Při izotermickém proudění je hranicí myšlena tenká stěna obklopující tekutinu
s okrajovými podmínkami danýmí pouze hodnotami rychlosti. V případě neizotermního
proudění závisí okrajové podmínky na každém konkrétním případu, tj. jestli je nutno použít
úplný model, částečně zjednodušený model nebo zjednodušený model:
·
úplný model - řeší se rozložení teploty v proudícím médiu, stěně trubky (vodivá
oblast) i v okolí (vzduch), okrajové podmínky jsou definované vnějšími tepelnými
podmínkami okolí
·
částečně zjednodušený model - řeší se rozložení teploty v proudícím médiu a stěně
trubky (vodivá oblast), je nutné definovat teplotu nebo hustotu toku tepla na vnější
stěně trubky,
·
zjednodušený model – řeší rozložení teploty v proudicím médiu s hranicí definovanou
stěnou nulové tloušťky (lze navíc respektovat tepelný odpor stěny zadané tloušťky)
s přesně definovanými neizotermickými vlastnostmi a okrajovými podmínkami.
Jednotlivé přístupy jsou zobrazeny na obr. 7.5.
hranice
úplný model
částečně zjednodušený
zjednodušený model
tekutina+stěna potrubí+okolí
model
tekutina+tenká stěna
tekutina+stěna potrubí
obr. 7.5 Proudění tekutiny potrubím
tekutina,
71
stěna potrubí,
okolí
Řešení přenosu tepla (konvekce, kondukce)
Hranice je pak nejvzdálenější plocha nulové tloušťky s definovanými okrajovými
podmínkami. Na hranici je možné nastavit okrajové podmínky pro přestup tepla (teplotu nebo
hustotu tepelného toku případně radiaci), rychlost (u pohybujících se a rotujících stěn),
smykové napětí, drsnost, podmínky pro příměs, chemické reakce, podmínky pro multifázové
proudění a volnou hladinu.
Pokud se řeší přenos tepla, na hranici se nastavují následující teplotní podmínky:
·
konstantní hustota toku tepla
·
teplota
·
konvektivní přenos tepla
·
externí radiační přenos tepla
·
kombinace externí radiace a externího konvekčního přenosu tepla
Rozhraní
může
být
zjednodušené
považována
za
jednostrannou
nebo
dvoustrannou stěnovou zónu.
dvoustranná stěnová zóna
two-sided wall
jednostranná stěnová zóna
one-sided wall
obr. 7.6 Hranice jako jednostranná a dvoustranná stěnová zóna
Je-li stěnová zóna dvoustranná (two-sided wall, tj. stěna tvoří rozhraní mezi dvěma oblastmi,
jako je rozhraní tekutina/stěna pro problém přenosu tepla), je možné modelovat vodivost
uvnitř hraničních stěn a vnitřních stěn (dvoustranných), viz obr. 7.6. Dále existuje také
možnost si vybrat, zda jsou nebo nejsou podmínky na dvoustranné stěně propojeny
(coupled).
Na obr. 7.7 je možno sledovat rozložení teploty v oblasti proudícího média a ve stěně,
kdy na vnější ploše stěny je dána konstantní podmínka pro teplotu. Šíření tepla je ovlivněno
materiálem, tj. vodivostí stěny a vody, a samozřejmě prouděním.
72
Řešení přenosu tepla (konvekce, kondukce)
obr. 7.7 Detail konce oblasti s rozložením teploty pro částečně zjednodušený přístup řešení.
7.4.1. Okrajové podmínky na tenké stěně
Podle předpokladu má stěna nulovou tloušťku. Je-li stěna nenulové tloušťky, lze
nastavit parametry pro výpočet tepelného odporu pro tenkou stěnu a modelovat tak tenkou
vrstvu materiálu mezi dvěma zónami. Např. lze modelovat účinek kousku plátu mezi dvěma
zónami tekutiny, potahování pevné látky, nebo kontaktního odporu mezi dvěma pevnými
oblastmi. FLUENT pak řeší 1D rovnici vodivosti, aby spočítal tepelný odpor definovaný
stěnou a generaci tepla ve stěně.
Aby se mohly tyto účinky zahrnout do
výpočtu
přenosu
tepla,
je
TENKÁ STĚNA
třeba
vnitřní plocha
(outer wall)
specifikovat typ materiálu, tloušťku stěny
vnější plocha
(inner wall)
a generaci tepla ve stěně. Tedy vybere
se
materiál, specifikuje
se
tloušťka
stěny. Tepelný odpor stěny je
Dx
, kde
l
Tw
resp. qw
l je vodivost materiálu stěny a Dx je
tloušťka stěny. Teplotní podmínka resp.
podmínka hustoty tepenlého toku bude
specifikována na vnější straně stěny, jak
je patrno na obr. 7.8. Dle konvence užité
buňky kapaliny
nebo pevné
vodivé stěny
ve Fluentu bude nazvána vnitřní plocha
(inner wall). Tw je konstantní teplota
lw
Dx
obr. 7.8 Okrajová podmínka na tenké stěně [15]
stěny. Je třeba poznamenat, že pro tenkou stěnu je možné definovat pouze konstantní
73
Řešení přenosu tepla (konvekce, kondukce)
tepelnou vodivost. Je-li třeba užít nekonstantní tepelnou vodivost pro nenulovou tloušťku, je
nutno stěnu definovat konkrétní geometrií a vysíťovat.
7.4.2. Okrajové podmínky na tenké dvoustranné stěně
Jestliže stěna má na každé straně kapalinu nebo pevnou stěnu, nazývá se tato stěna
dvoustranná (two-sided wall), a je schématicky zobrazená na obr. 7.9.
TENKÁ STĚNA
wall
shadow wall
Tw1
qw1
buňky
kapaliny nebo
pevné vodivé
stěny
Tw2
qw2
lw1
lw2
buňky
kapaliny nebo
pevné vodivé
stěny
obr. 7.9 Okrajová podmínka na stěně se dvěma povrchy [15]
Když je vložena z Gambitu síť s tímto typem stěny do Fluentu, vytvoří se automaticky
„shadow“ zóna tak, že každá strana stěny je stěnová zóna. V panelu WALL se ukáže jako
Shadow Face Zone. Pak lze definovat odlišné tepelné podmínky na každé zóně označené
WALL a SHADOW WALL, nebo propojit (coupled) obě zóny:
·
Při propojení zón je třeba vybrat Coupled option v Thermal Conditions (tento
parametr se objeví ve WALL panelu, když stěna je dvoustranná). Žádné doplňující
tepelné okrajové podmínky nejsou požadovány, protože přestup tepla bude řešen
přímo z rovnic pro sousedící buňky. Lze ale definovat typ materiálu, tloušťku stěny a
generaci tepla pro výpočty tepelného odporu, jak bylo uvedeno výše. Parametry
odporu tepla budou automaticky nastavené na „shadow“ stěnové zóně.
·
Při odlišných (nepropojených) stěnových zónách mohou být definovány odlišné
tepelné podmínky na každé z nich. Je třeba vybrat Temperature nebo Heat Flux
(Convection a Radiation nejsou možné pro dvoustrannou stěnu). Obě nepropojené
stěny mohou mít odlišnou drsnost a jsou navzájem izolovány. Pokud je třeba
specifikovat nenulové tloušťky stěn pro nepropojené zóny, tepelné podmínky budou
definovány na vnější ploše nenulových stěn, jak je patrno z obr. 7.9, kde Tw 1 a Tw 2 je
74
Řešení přenosu tepla (konvekce, kondukce)
teplota ( qw 1 a qw 2 je hustota tepelného toku) definovaná na jedné a druhé stěně.
lW 1 a lW 2 jsou tepelné vodivosti na nepropojených nenulových stěnách. Mezera
mezi stěnami není částí modelu, je pouze z ilustrativních důvodu zahrnuta do
obrázku.
7.4.3. Výpočet teploty a hustoty tepelného toku na stěně
Výpočet tepelného toku při zadané teplotě
Je-li na stěně definována okrajová podmínka teplotou, pak lze počítat hustotu
tepelného toku z tekutiny do stěny a pro turbulentní proudění je
q = hf (Tw - Tf ) + qrad
kde
hf
( 7.4.1)
lokální součinitel přestupu tepla v oblasti, kde je tekutina, je určen
logaritmickou závislostí v blízkosti stěny,
Tw
teplota povrchu stěny,
Tf
lokální teplota tekutiny v buňce blízko stěny,
q
konvekční hustota tepelného toku ze stěny,
qrad
radiační hustota tepelného toku.
Jestliže stěna hraničí s buňkou vodivé stěny, pak Fluent počítá přenos tepla k okraji stěny
jako
q=
kde
lcw
(Tw - Tcw ) + qrad
Δn
( 7.4.2)
lcw
tepelná vodivost vodivých stěn,
Tcw
lokální teplota vodivých stěn,
Dn
vzdálenost mezi povrchem stěny a středem buňky vodivé stěny.
Příklad
Řešte proudění v trubce s přestupem tepla stěnou, uvažujte vedení tepla stěnou a přestup
tepla mezi stěnou a proudícím médiem
75
Řešení přenosu tepla (konvekce, kondukce)
Tw=400 K, ocel
q=0 Wm-2
-2
q=0 Wm
T=300K
v=0.5 m/s
voda
Oblast proudění
s
d
obr. 7.10 Schéma řešené oblasti
Geometrie oblasti:
Vnitřní průměr trubky d [mm]
25
Tloušťka stěny trubky s [mm]
3
Délka trubky l [m]
0.9
Fyzikální vlastnosti kapaliny a stěny:
Hustota ρ [kgm-3]
Měrná tepelná kapacita cp [J(kgK)-1]
-1
Tepelná vodivost l [W(mK) ]
-1
Dynamicá viskozita h [kg(ms) ]
voda
ocel
998.2
8030
4182
502.48
0.6
16.27
0.001003
Okrajové podmínky:
Teplota vody na vstupu T [K]
300
Rychlost vody na vstupu u [ms-1]
0.5
Vnější teplota stěny Tw [K]
400
Teplota vnější plochy tenké stěnyTw = Tb - inner [K]
400
Hustota tepelného toku q [Wm-2]
0
Matematický model:
Reynoldsovo číslo Re [1]
25 000
Re =25 000 Þ proudění je turbulentní
Výsledky:
76
Řešení přenosu tepla (konvekce, kondukce)
obr. 7.11 Detail sítě konce trubky a rozložení teploty
390
vstup
vystup
370
350
0
T [ C]
x=0.45
330
310
290
0
0.005
0.01
0.015
0.02
0.025
0.03
r [m]
obr. 7.12 Teplotní profil v radiálním směru na vstupu, výstupu a uprostřed ve vzdálenosti
x = 0.45 m
77
Řešení přenosu tepla (konvekce, kondukce)
400
250000
395
240000
T [K]
390
230000
teplota-outer
385
teplota-inner
220000
380
hustota toku
210000
375
200000
370
190000
365
180000
360
170000
355
160000
350
150000
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
x [m]
obr. 7.13 Teplotní profil v podélném směru na vnitřní ploše (outer), vnější poše (inner) a
hustota tepelného toku na stěně
Výpočet teploty při zadané hustotě tepelného toku
Pro výpočet teploty povrchu stěny přiléhající k živé buňce se využije rovnice ( 7.4.1) a
uživatelem zadaná hustota tepelného toku.
Tw =
q - q rad
+ Tf
hf
( 7.4.3)
Hustota tepelného toku je kladná pro směr ze stěnových buněk do buněk živých a naopak.
Jestliže stěnová buňka hraničí s buňkou vodivé stěny, pak se teplota povrchu stěny počítá
dle následujícího vztahu.
Tw =
(q - q rad ) × Dn + T
lcw
( 7.4.4)
cw
Adiabatické stěny jsou stěny bez prostupu tepla, mohou být nadefinovány zadáním
nulové hustoty tepelného toku jako okrajové podmínky na stěně.
Okrajová podmínka na stěně definovaná konvektivním přenosem tepla
Software požaduje zadání externího součinitele přestupu tepla a také zadání externí
teploty uživatelem, což je ve většině případů teplota vnějšího okolí. Slovy externí součinitel
přestupu tepla je součinitel přestupu tepla ze stěny do okolí. Tepelný tok je pak počítán dle
následujícího vztahu
78
Řešení přenosu tepla (konvekce, kondukce)
q = hf (Tw - Tf ) + q rad = h ext (Text - Tw )
kde
hext =
l
Dx
( 7.4.5)
externí součinitel přestupu tepla, Dx tloušťka stěny
Text
externí teplota (teplota okolí),
q rad
hustota radiačního toku tepla
Další výpočty souvisejí s radiací a nebudou v této učebnici vysvětlovány.
79
Modelování proudění příměsí
8. Modelování proudění příměsí
8.1.
Transportní rovnice pro přenos příměsí
Typickým příkladem pro modelování příměsi je šíření znečišťujících látek z komína
nebo jiných zdrojů polutantů. K výpočtu je třeba znát parametry zdroje, proudění vzduchu
v atmosféře, povrch terénu, viz obr. 8.1
Faktory ovzduší
· směr a rychlost větru
· vertikální teplotní gradient
· turbulence
· ostatní meteorologické
prvky
Parametry zdroje:
· množství emisí
· geometrie
· poloha, výška
Faktory prostředí
· konfigurace terénu
· zástavba
· drsnost povrchu
· zdroje tepla, teplotní toky
obr. 8.1 Faktory, ovlivňující šíření znečištění v atmosféře.
Matematický model je dán rovnici kontinuity, pohybovou rovnicí, rovnicemi pro přenos
turbulentních veličin a navíc počítá lokální hmotnostní zlomky příměsí Yi ¢ . Hmotnostní
zlomek příměsi je definován vztahem
Yi ¢ =
kde
mi ¢ r i ¢Vi ¢ r i ¢
=
=
ai¢
m
rV
r
( 8.1.1)
mi ¢ [kg]
hmotnost příměsi i ¢
m [kg]
celková hmotnost směsi
Yi ¢ [1]
hmotnostní zlomek příměsi i ¢ ve směsi
a i ¢ [1]
objemový zlomek příměsi i ¢ ve směsi
Další veličina, která se užívá se spojení s šířením příměsi je molární koncentrace Ci ¢
[kmol.m-3]. Koncentrace definovaná vztahem Mi ¢Ci ¢ je uváděna v jednotkách [kgkmol1
.kmolm-3=kgm-3]. Označení ppm běžné při vyhodnocování koncentrací definuje miliontinu
dané hodnoty (analogie procenta, může se vztahovat k hmotnostnímu nebo objemovému
zlomku).
80
Modelování proudění příměsí
Bilanční rovnice přenosu příměsi (tj. hmotnostního zlomku) v konzervativní formě má
tvar [24], [28]
¶
¶(rYi ¢ ) ¶ (r u jYi ¢ )
=+
J j ,i ¢ + Ri ¢ + Si ¢
¶xi
¶x j
¶t
( 8.1.2)
kde u i je časově středovaná složka rychlosti proudění a na pravé straně je Ri ¢ rychlost
produkce příměsí i ¢ vlivem chemické reakce a Si ¢ rychlost tvorby přírůstku z distribuované
příměsi [1]. Výše uvedená rovnice platí pro N - 1 příměsí, kde N je úplný počet komponent
prezentovaných v matematickém modelu. Distribuce příměsí může být realizována za
různých podmínek, obecně lze rozlišovat distribuci za laminárního nebo turbulentního
proudění [24]. J j ,i ¢ představuje difúzní tok i ¢ -té komponenty směsi, který se liší pro laminární
a turbulentní proudění.
Difúzní tok pro laminární proudění
V předchozí rovnici J i ¢ představuje difúzní tok i ¢ -té složky jednotkou plochy, který je
definovaný vztahem
J i ¢ = - r Di ¢,m
¶Yi ¢
¶x j
[kgm-2s-1]
( 8.1.3)
kde Di ¢,m je difúzní koeficient i ¢ -té příměsi ve směsi.
Difúzní tok pro turbulentní proudění
Při turbulentním proudění pro vyjádření difúzního toku jednotkou plochy i ¢ -té složky se
uplatňuje vztah
æ
m ö ¶Y
J i ¢ = -çç r Di ¢,m + t ÷÷ i ¢
Sct ø ¶x j
è
[kgm-2s-1]
( 8.1.4)
kde Sct je Schmidtovo turbulentní číslo, jehož hodnota je 0,7.
8.2.
Definice zdroje příměsi
Parametry zdroje zahrnují data o typu, velikosti a umístění zdroje v oblasti, složení a
množství vypouštěných znečisťujících látek. Je nutné definovat polohu zdroje v kartézském
souřadném systému, jeho objemový průtok, teplotu, rychlost a její směr na výstupu ze zdroje
a koncentraci (hmotnostní zlomek) jednotlivých příměsí na výstupu ze zdroje. Vydatnost
zdroje může být konstantní nebo se může měnit v závislosti na čase. V případě havarijních
úniků je obtížné definovat množství a koncentraci látky, která se bude v okamžiku havárie
v daném místě nacházet. Je nejprve nutné identifikovat a porozumět dějům, které se
81
Modelování proudění příměsí
odehrávají v samotném zdroji, jako je například výbuch, odpařování, současný únik plynné a
kapalné fáze, atd. K výpočtu potřebných parametrů pro různé typy zdrojů lze použít
empirické vztahy. Zadání zdroje je umožněno definováním objemových zdrojových členů
hmotnosti, hybnosti, energie, turbulence a dalších skalárních veličin. Je možné také
definovat tyto zdroje jako nekonstantní, závislé na čase případně dalších veličinách [8], [24].
V praktických úlohách se nejčastěji vyskytují tyto typy zdrojů: bodový zdroj (komín),
plošný zdroj a objemový zdroj. Tyto typy se liší geometrickými vlastnostmi a definicí
příměsi vystupující ze zdroje. Např. komín představuje válcovou vertikální konstrukci, ze
které jsou vypouštěny příměsi plynného i pevného charakteru. Tento typ zdroje je vhodný,
pokud je předmětem zájmu modelování obtékáni komína jako překážky v řešené oblasti. Na
příkladu komína budou definovány všechny varianty definice zdroje, jako je obvyklý inlet,
umístěný na ploše definující výstup z komína nebo na hranici oblasti a objemový zdroj
umístěný uvnitř oblasti.
Inlet2
obr. 8.2 Zdroje příměsi:
Inlet1
Source
u
plošný zdroj typu Inlet jako výstup z komína nebo otvor
v horní ploše oblasti
objemový zdroj typu Source uvnitř oblastk
8.2.1. Zdroj-inlet
Pokud není nutné modelovat proudění v samotném komíně, případně vliv přestupu
tepla přes stěny komína, pak se konstruuje jeho geometrie již v pre-procesoru, horní průřez o
ploše S se pak definuje jako vstup inlet, komín se geometricky netvoří vůbec. Na vstupu se
zadají okrajové podmínky pro plošný zdroj ve výšce komína, které jsou: hmotnostní zlomek
příměsi, teplota na výstupu z komína, rychlost na výstupu z komína a turbulentní parametry
(např. hydraulický průměr, intenzita turbulence).
8.2.2. Objemový zdroj
Pokud je průměr komína příliš malý ve vztahu k rozměrům oblasti nelze jej definovat
82
Modelování proudění příměsí
pomocí podmínky inlet, je lepší jej nahradit objemovým zdrojem. Může obsahovat pouze
jednu buňku nebo skupinu buněk, které však musí tvořit oddělenou zónu. Tu je možné
vytvořit už v preprocesoru při přípravě sítě nebo dodatečně přímo v prostředí FLUENTu
pomocí nástrojů adaptace sítě a tvorby registru. Protože zdrojové členy jsou definovány na
jednotku objemu, je nutné nejprve určit velikost objemu buněk v dané zóně pomocí integrace
(Report-Volume Integral). Pak se objemový zdroj hmotnosti
(pro jednu i více příměsí)
definuje vztahem [24]
mass source =
Qm
V
[kgm-3s-1]
( 8.2.1)
V rovnici kontinuity se definovaný zdroj hmotnosti objeví jako zdrojový člen Sm . Obdobným
způsobem je možné definovat i zdrojový člen v pohybové rovnici jako momentum source X,
Y, Z
momentum source =
ui Qm
V
[Nm-3=ms-1kgs-1m-3]
( 8.2.2)
kde ui je složka rychlosti ve směru x, y, z. V pohybové rovnici se zdrojový člen objeví
v hmotnostní síle F . Tento člen definuje hybnost příměsi. Pokud je roven nule, pak příměs,
která je vložena do oblasti pouze jako zdrojový člen v rovnici kontinuity, je unášena okolním
proudem.
Analogicky lze definovat zdrojový člen v rovnici energie S h , v transportní rovnici pro
kinetickou turbulentní energii S k a rychlost disipace S e . V případě zdrojů, jejichž parametry
se mění v závislosti na čase, je možné závislosti popsat pomocí uživatelských funkcí UDF.
8.3.
Fyzikální vlastnosti plynů a jejich směsí
Nejdříve je nutno definovat fyzikální vlastnosti jednotlivých plynů tvořících směs. Pokud
se řeší izotermní proudění (míchání látek o konstantní teplotě), pak tyto vlastnosti mohou být
konstantní. V případě neizotermního proudění je možno uvažovat jednotlivé vlastnosti jako
funkce teploty. Teprve potom se určí celková vlastnost směsi.
8.3.1. Hustota
Pro izotermní proudění plynů s konstantní hustotou, nebo neizotermní s hustotou
vyjádřenou jako funkce teploty (ne jako ideální plyn), je hustota definována
r=
kde
1
Y
åi ¢ r i ¢
i¢
( 8.3.1)
Yi ¢ je hmotnostní zlomek příměsi i ¢ ve směsi
83
Modelování proudění příměsí
ri ¢ je hustota příměsi i ¢ ve směsi
Hustotu pro neizotermní proudění je také možno zadat polynomem nebo po částech
lineární funkcí
r i ¢ (T ) = A1 + A2T + A3T 2 + ... pro Tmin < T < Tmax
r i ¢ (T ) = r i ¢,n +
(8.3.2)
r i ¢,n +1 _ r i ¢,n
(T - Tn )
Tn +1 - Tn
Je také možno definovat vlastní závislost pomocí UDF (User Defined Function),
Pokud je hustota plynu základního proudění a všech příměsí zadána stavovou
rovnicí
ri ¢ =
(p
op
+ p)
(8.3.3)
rT
pak hustota směsi je dána v daném objemu součtem hustot příměsí modifikovaných podílem
hmotnostního zlomku a molekulové hmotnosti (což odpovídá objemovému zlomku)
r=
pop + p
Y
RT åi ¢ i ¢
Mi ¢
( 8.3.4)
kde M i ¢ je molekulová váha příměsi i ¢ ve směsi.
Shrnutí
hustota i ¢ -té příměsi
hustota směsí
r i¢ = konst
r i ¢ (T ) = A1 + A2 × T + A3 × T 2 + ...
r i ¢ (T ) = r i ¢,n +
ri ¢
(p
=
op
r=
r i ¢,n +1 _ r i ¢,n
(T - Tn )
Tn +1 - Tn
+ p)
r=
rT
1
Y
åi ¢ ri ¢
i¢
pop + p
Y
RT åi ¢ i ¢
Mi ¢
8.3.2. Viskozita
Kinematická viskozita jednotlivých plynných látek může být konstantní, tj. nezávisí na
teplotě, nebo je vyjádřena funkcí teploty (např. jako polynom n -tého řádu, po částech
lineární funkce nebo jiná funkční závislosti na teplotě),
viskozita i ¢ -té příměsi
viskozita směsí
mi ¢ = konst
m = åi ¢ Yi ¢ m i¢
84
Modelování proudění příměsí
m i ¢ (T ) = A1 + A2 × T + A3 × T 2 + ...
m i ¢ (T ) = m i ¢,n +
m i ¢,n +1 _ m i ¢,n
(T - Tn )
Tn +1 - Tn
X i ¢ × mi ¢
m =å
i¢
Ideální plyn
å
j¢
Xi¢
é æ m ö1/ 2 æ M ö1/ 4 ù
ê1 + ç i ¢ ÷ × ç j ¢ ÷ ú
çM ÷ ú
ê çè m j ¢ ÷ø
è i¢ ø û
ë
é æ M ¢ öù
ê8 × çç 1 + i ÷÷ú
êë è M j ¢ øúû
2
1/ 2
kde X i ¢ je molový zlomek příměsi i ¢ (počet molů příměsi v jednom molu směsi).
8.3.3. Měrná tepelná kapacita
Podobně předchozím vlastnostem měrná tepelná kapacita plynů je vyjádřena např.
funkcí teploty jako polynom n -tého řádu
c p,i ¢ (T ) = A1 + A2 × T + A3 × T 2 + ... pro
Tmin < T < Tmax
(8.3.5)
a tepelná kapacita směsí je dána vztahem
c p = å Yi ¢c p,i ¢
(8.3.6)
i¢
8.3.4. Tepelná vodivost
Tepelná vodivost jednotlivých plynných látek je konstantní nebo je vyjádřena funkcí
teploty jako polynom n -tého řádu
li ¢ (T ) = A1 + A2 × T + A3 × T 2 + ... pro
Tmin < T < Tmax
(8.3.7)
pak při definici hustoty pro ideální plyn je tepelná vodivost směsí je dána vztahem
X i ¢ × li ¢
l=å
i¢
å
X i¢
j¢
é æ l ö1 / 2 æ M ö1 / 4 ù
ê1 + ç i ¢ ÷ × ç j ¢ ÷ ú
çM ÷ ú
ê çè l j ¢ ÷ø
è i¢ ø û
ë
é æ
M öù
ê8 × çç1 + i ¢ ÷÷ú
êë è M j ¢ øúû
2
(8.3.8)
1/ 2
8.3.5. Standardní slučovací entalpie a entropie
V případech, kdy se řeší proudění s uvažováním chemické reakce, je nutné definovat
standardní slučovací entalpii (slučovací teplo) h0j ¢ pro každou příměs j ¢ . Tato vlastnost
85
Modelování proudění příměsí
slouží k definování entalpie směsi
T
é
ù
H = å m j ¢ êh 0j¢ + ò c p, j ¢dT ú
êë
úû
Tref , j ¢
(8.3.9)
kde Tref , j ¢ je referenční teplota, za které je definována h0j ¢ .
Jestliže se uvažuje vratná chemická reakce, je nutné definovat standardní slučovací
entropii s 0j ¢ pro každou příměs j ¢ . Entropie směsi je definována
T
c p, j ¢
é
ù
S = å m j ¢ ês 0j ¢ + ò
dT ú
Tref , j ¢ T
j¢
ë
û
(8.3.10)
Hodnoty standardní slučovací entalpie a entropie lze vyhledat v tabulkách.
8.4.
Okrajové podmínky pro
příměsi
na vstupu,
výstupu a stěně
Na vstupu je definován hmotnostní zlomek příměsi. Je třeba poznamenat, že je
definován hmotnostní zlomek pouze pro N-1 příměsí, přitom poslední je doplňkem do 1,
neboť součet všech hmotnostních zlomků musí být roven 1.
Pro podmínku Outfow je nutné definovat hmotnostní zlomek příměsi pro případ zpětného
proudění.
Na stěně je předdefinována podmínka nulového gradientu příměsi, tj. příměs neprochází
stěnou. Výjimkou je příměs reagující s povrchem stěny. Je možno definovat také hmotnostní
zlomek.
Příklad
Řešte šíření příměsi lehkého a těžkého plynu v oblasti.
Zdroj 2
y
u
Oblast proudění
Zdroj 1
x
obr. 8.3 Schéma řešené oblasti
86
z
Modelování proudění příměsí
Geometrie oblasti
délka oblasti x [m]
3.5
výška oblasti y [m]
0.5
šířka oblasti z [m]
1.5
Fyzikální vlastnosti:
vzduch
Zdroj 1 CO2
Zdroj 2 SO2
Hustota ρ [kgm ]
1.225
1.7878
2.77
Měrná tepelná kapacita cp [J(kgK)-1]
1006.43
840.37
622.28
Tepelná vodivost l [W(mK)-1]
0.0242
0.0145
0.0104
Dynamicá viskozita h [kg(ms)-1]
1.7894e-05
1.37e-05
1.2e-05
vzduch
Zdroj 1 SO2
Zdroj 2 CO2
rychlost u [ms ]
1
0.5
0.5
Intenzita turbulence I [%]
2
5
5
Hydraulický průměr d h [m]
0.4
0.0794
0.0794
Hmotnostní zlomek CO2 h [1]
0
0
1
Hmotnostní zlomek SO2 h [1]
0
1
0
-3
Okrajové podmínky:
-1
Plyny jsou o rozdílné hustotě, proto se bude sledovat také vliv vztlakové síly na šíření plynů.
Matematický model:
Reynoldsovo číslo
Re =11 784 Þ proudění je turbulentní
Výsledky:
statický tlak
amplituda rychlosti
obr. 8.4 Detail sítě konce trubky a rozložení teploty
87
Modelování proudění příměsí
YCO2=0.01
YSO2=0.05
obr. 8.5 Izoplochy hmotnostního zlomku v podélném a šesti příčných řezech
Vztlakové síly aktivované
Vztlakové síly neaktivované
obr. 8.6 Prostorové izoplochy hmotnostního zlomku YCO2=0.01 a YSO2=0.05
88
Vícefázové modely
9. Vícefázové modely
9.1. Vícefázové modely obecně
Velké množství aplikací v přírodě a v průmyslových technologiích se týká směsí fází.
Fázemi se předpokládá:
·
plyn
·
kapalina
·
pevná látka
Ale koncepce fází je aplikována v širším smyslu. V multifázovém proudění je fáze definována
jako identifikovatelná třída materiálu, který má
částečnou inertní odezvu na interakci
s prouděním a potenciálem pole, ve kterém se vyskytuje. Např. pevné částice různých
velikostí téhož materiálu mohou být považovány za různé fáze, protože každé seskupení
částic o téže velikosti má podobné dynamické vlastnosti v proudovém poli. Do vícefázového
proudění zahrnujeme řešení následujících problémů: kavitace, aerace, nukleární reaktory,
bezpečnost nukleárních reaktorů, vlnění vodní hladiny, extrakce, emulsifikace,
separace, homogenizace, promíchávání, hydraulická doprava, sedimentace, flotace,
cyklony, pneumatická doprava, fluidizační pole reaktoru, sila, spalování, odpařování
aj. [30], [31].
Prakticky se lze setkat s různými variantami vícefázových systémů, které jsou
upřesněny v následujícím přehledu:
plyn – kapalina nebo kapalina - kapalina
·
proudění bublin plynu nebo velkých kapek kapaliny ve spojitém prostředí
·
proudění kapek ve spojité fázi plynu
·
pomalé proudění velkých bublin
·
proudění s volnou hladinou, s jasně definovanou hladinou
plyn – pevná látka
·
proudění pevných částic v plynu
·
pneumatická doprava
·
fluidizační pole
kapalina – pevná částice
·
proudění kalu
·
sedimentace
třífázové proudění
·
kombinace výše uvedených variant
89
Vícefázové modely
proudění s částicemi,
proudění s volnou
hydraulická a pneumatická
bublinami, kapkami
hladinou
doprava
obr. 9.1 Přiklady vícefázového proudění [15]
Přístupy k modelování vícefázového proudění
V současné době existují dva přístupy pro modelování vícefázového proudění
·
Euler - Lagrangeův přístup
·
Euler - Eulerův přístup
9.1.1. Euler-Lagrangeův přístup
Při Euler - Lagrangeově přístupu je tekutá fáze uvažována jako kontinuum a je řešena
Navierovými - Stokesovými rovnicemi, zatímco dispergovaná fáze (částice) je řešena
stopováním velkého počtu částic, bublin nebo kapek v proudovém poli. Tato dispergovaná
fáze může vyměňovat moment, hmotu a energii se spojitou fází. Základním předpokladem je,
že v tomto modelu dispergovaná fáze zaujímá malý objemový zlomek, ačkoliv pro hmotnost,
resp. hmotnostní průtok to nemusí platit ( Qm,částic ³ Qm,tekutiny ). Trajektorie (dráhy) částic nebo
kapek jsou počítány individuálně v předdefinovaných intervalech během výpočtu spojité fáze.
Toto umožňuje modelovat proudění částic ve sprejích, sušičkách, spalování uhelných a
kapalných paliv a částicemi ovlivněné proudění. Je nevhodný pro modelování směsi kapalina
- kapalina, fluidizačního lože a dalších aplikací, kde objemový zlomek druhé fáze není
zanedbatelný.
9.1.2. Euler-Eulerův přístup
Při Euler-Eulerově přístupu jsou různé fáze řešeny matematicky jako vzájemně se
prostupující kontinua. Protože objem jedné fáze není překryt objemem druhé fáze, je
zaveden pojem objemového zlomku fáze. Tyto objemové zlomky se předpokládají jako
funkce spojité v čase a prostoru a jejich součet je roven 1. Rovnice jsou definovány pro
každou fázi.
90
Vícefázové modely
9.2. Vícefázové matematické modely
Vícefázové modely [15] umožňují modelování většího počtu oddělených, ale vzájemně se
ovlivňujících fází. Fáze mohou být kapalné, plynné a pevné v různých kombinacích (plyntekutina, tekutina-tekutina, plyn-pevná látka, kapalina-pevné částice, trojfázové proudění
(kombinace předchozích)). V literatuře se diskutují tři multifázové modely, jsou to:
·
VOF model
·
Model směsi
·
Eulerův model
Každý z těchto modelů má své typické aplikace a také své klady i zápory. Proto
v následujících odstavcích je uveden krátký přehled těchto modelů, jejich kompletní popis lze
najít v literatuře [15].
VOF MODEL je vhodný pro stratifikované (vrstvené) proudění a proudění s volnou
hladinou. Tímto modelem se může řešit proudění dvou a více nesmísitelných kapalin
řešením pohybové rovnice a sledováním objemového zlomku každé kapaliny v oblasti.
Typické aplikace zahrnují předpověď odtržení proudu, pohyb velkých bublin v kapalině,
pohyb kapaliny za hrází a ustálené nebo neustálené sledování jakýchkoliv rozhraní kapalina
- plyn.
MODEL SMĚSI (MIXTURE MODEL) je zjednodušený vícefázový model, který lze použít
k modelování vícefázového toku, kde se jednotlivé fáze posouvají různou rychlostí.
Předpokládá se ale lokální rovnováha na krátkém prostorovém délkovém měřítku. Vazba
mezi fázemi musí být silná. Toho se může využít také k modelování homogenního
vícefázového proudění s velmi silnou vazbou a fázemi pohybujícími se stejnou rychlostí.
Model směsi může modelovat n-fází (tekutina nebo částice) řešením pohybové rovnice,
rovnice kontinuity a rovnice energie pro směsi, rovnice objemového zlomku pro druhou fázi
(dispergovanou) a algebraického výrazu pro relativní rychlosti. Typická aplikace zahrnuje
sedimentace, cyklónové separátory, částice s nízkým zatížením a bublinkovité proudění, kde
objemový zlomek plynu je nízký.
EULERŮV MODEL dovoluje modelování vícenásobných oddělených interaktivních fází.
Fáze mohou být tekutina, plyn a pevné látky v nějaké kombinaci. U Eulerova multifázového
modelu je počet dalších fází limitován pouze požadavky na paměť a konvergenci řešení. To
znamená, že lze modelovat libovolný počet dalších fází, pokud je k dispozici dostatečná
paměť počítače.
91
Vícefázové modely
V dalším textu je uvedena krátká rozvaha, kdy a za jakých podmínek je vhodnější použit
model směsi a kdy Eulerův model.
·
Pokud se očekává velký výskyt dispergované fáze, pak je preferován model směsi
kvůli menší náročnosti na výpočet. Je-li dispergovaná fáze koncentrována jen v určité
oblasti, je lépší použít Eulerův model.
·
Je-li nutno aplikovat odpor proti pohybu v systému (budˇ přímo ve Fluentu nebo
pomocí uživatelem definovaná funkce UDF), pak Eulerův model může obvykle
poskytovat přesnější výsledky než model směsí. Jestliže nejsou známy interfázové
odporové koeficienty nebo jejich použitelnost na zkoumaný systém, pak model směsi
může být lepší volbou.
·
Pokud se řeší jednodušší problém, který požaduje menší výpočtovou náročnost, je
model směsi lepší volbou, jelikož řeší menší počet rovnic než model Eulerův. Když
přesnost je důležitější než výpočtová náročnost, pak Eulerův model zaručí kvalitnější
výsledek. Musí se však pamatovat na to, že komplexnost Eulerova modelu může být
příčinou menší stability řešení než u modelu směsi.
Pro snazší identifikaci multifázového proudění se zavádějí pojmy jako je:
·
hmotnostní konzistence je definována parametrem b =
a d rd
(mass density ratio),
a c rc
kde d je index dispergované fáze, c je index nosné fáze, a je objemový zlomek, r je
hustota.
·
podíl hustot g =
rd
(material density ratio), který je větší než 1000 pro fáze plyn rc
pevná látka, kolem 1pro fáze kapalina - pevná látka a menší než 0.001 pro fáze plyn kapalina.
1
·
b
L æ p 1+ k ö3
průměrná vzdálenost mezi částicemi
=ç
÷ , kde k = .
g
dd è 6 k ø
Informace o těchto parametrech je důležitá pro určení chování dispergované fáze. Např. pro
směs plyn - částice s b » 1 je prostor mezi částicemi
jako izolovaná.
Interakce mezi fázemi se pak dělí na tři kategorie:
92
L
» 8 . Částice se tedy bude chovat
dd
Vícefázové modely
·
Pro malé hodnoty b je vztah mezi fázemi jedním směrem, tj. spojitá fáze ovlivňuje
částice přes odpor a turbulenci, ale částice neovlivňují spojitou fázi. Model diskrétní fáze,
směsi a Eulerův se vypořádá s tímto typem velmi dobře. Protože Eulerův model je velmi
náročný, je možno použít model diskrétní fáze nebo model směsi.
·
Pro střední hodnoty b je vztah mezi fázemi dvousměrný, tj. spojitá fáze ovlivňuje
částice přes odpor a turbulenci a částice ovlivňují spojitou fázi pomocí středního momentu a
turbulence. Modely diskétní fáze, směsi a Eulerův jsou vhodné pro řešení, ale pro výběr
modelu je nutné uvážit další parametry, jako je Stokesovo číslo.
·
Pro vyšší hodnoty b je vztah mezi fázemi dvousměrný a navíc se uvažuje tlak a
viskózní napětí z důvodu existence částic (čtyřsměrný vztah). Pro tento typ je vhodný
Eulerův model.
Stokesovo číslo je definováno jako vztah mezi časovou odezvou částic a časovou
odezvou systému:
St =
kde t d =
td
ts
( 9.2.1)
r d d d2
a ts je založeno na charakteristické délce Ls a charakteristické rychlosti Vs
18 mc
a platí t s =
Ls
. Pro St áá1 budou částice unášeny proudem a všechny tři modely budou
Vs
použitelné. Je tedy vhodné vybrat rychlejší a levnější variantu modelu, tj. model směsi
(mixture model). Pro St ñ1 se částice budou pohybovat nezávisle na proudu a je možno
použít model diskrétní fáze a Eulerův model. Pro St » 1 je možno vybrat kterýkoliv ze tří
modelů.
9.2.1. Vícefázový model směsi (mixture model)
Mixture model je zjednodušený vícefázový model, který může být použit pro proudění, kde
fáze se pohybují odlišnými rychlostmi, ale předpokládá se rovnováha při krátkém
prostorovém měřítku. Vazba mezi fázemi je velmi silná. Tímto modelem lze také modelovat
homogenní multifázové proudění s velmi silnou vazbou mezi fázemi pohybujícími se stejnou
rychlostí. Model je navíc schopen modelovat i proudění neNewtonských kapalin.
Model může řešit proudění n fází kapaliných nebo částic řešením pohybové rovnice a
rovnice kontinuity pro směs, rovnice pro objemový zlomek druhých fází a algebraického
vztahu pro relativní rychlosti. Typickými aplikacemi jsou sedimantece, cyklóny, proudění
s částicemi o malém zatížení a proudění bublin o malém objemovém zlomku.
93
Vícefázové modely
Model umožňuje řešit prolínání fází. Za tím účelem jsou definovány objemové zlomky fáze
q a p pro daný konečný objem označené jako a q a a p , jejichž hodnota je mezi nulou a
jedničkou a závisí na velikosti objemu, který zabírá daná fáze. Fáze se mohou pohybovat
různými rychlostmi, přitom se aplikuje koncepce relativních (slip) rychlostí.
Rovnice kontinuity pro směs je dána vztahem
¶r m ¶ (r mum, j )
+
=0
¶t
¶x j
( 9.2.2)
kde um, j jsou složky rychlosti zprůměrované podle hmotnosti
n
u m, j =
åa
k =1
k
r k uk , j
( 9.2.3)
rm
a r m je hustota směsi
n
rm = åa k rk
( 9.2.4)
k =1
kde a k je objemový zlomek fáze k .
Rovnice zachování hybnosti pro směs je získána sečtením rovnic zachování hybnosti
pro jednotlivé fáze
¶ (rmum,i ) ¶
¶p
¶
+
+
(
rmum,i um, j ) = ¶t
¶x j
¶xi ¶x j
¶
+ rfi +
¶x j
æ æ ¶um,i ¶um, j
ç mm ç
+
ç ç ¶x j
¶xi
è è
ö
æ n
çç å a k rk udr ,km,i udr ,,k , j ÷÷
ø
è k =1
ö
¶u
÷ - md ij 2 m,l
÷
3 ¶x l
ø
ö
÷+
÷
ø
( 9.2.5)
kde n je počet fází, fi jsou složky vnějších hmotnostních sil, m m je dynamická viskozita
směsi
mm =
n
åa
k= 1
k
mk
( 9.2.6)
a udr ,k ,i je složka unášivé rychlosti
udr ,k , i = uk ,i - um,i
( 9.2.7)
Relativní (slip) rychlost je definována jako rychlost sekundární fáze p k rychlosti sekundární
fáze q
u p,q ,i = u p,i - uq, i
( 9.2.8)
Jestliže hmotnostní zlomek fáze k je dán vztahem
ck =
a k rk
rm
( 9.2.9)
94
Vícefázové modely
pak unášivá rychlost a relativní (slip) rychlost jsou ve vztahu
n
udr , p,i = u p,q ,i - å ck uk ,q ,i
( 9.2.10)
k =1
Upřesnění unášivé rychlosti je závislé na definování odporových sil částic atd.
Rovnice objemového zlomku sekundární fáze je dána rovnicí:
¶ (a p r p ) ¶ (a p r pum, j )
¶ (a p r pudr , p, j ) n
+
=+ å Qm,qp - Qm, pq
¶t
¶x j
¶x j
q =1
( 9.2.11)
Příklad
Zobrazte objemové zlomky vody při proudění vody a vzduchu.
y
u
Oblast proudění
z
Vstup-vzduch
x
obr. 9.2 Schéma řešené oblasti
Geometrie oblasti
délka oblasti x [m]
3.5
výška oblasti y [m]
0.5
šířka oblasti z [m]
1.5
Fyzikální vlastnosti:
Hustota ρ [kgm-3]
-1
Dynamicá viskozita h [kg(ms) ]
Voda
Zdroj vzduch
998
1.225
0.001003
1.7894e-05
Okrajové podmínky:
Vstup
intenzita turbulence I [%]
mixture
2
hydraulický průměr d h [m]
mixture
0.4
rychlost u [ms-1]
voda
1
95
Vícefázové modely
vzduch
1
voda
-
vzduch
0
intenzita turbulence I [%]
mixture
5
hydraulický průměr d h [m]
mixture
0.0794
rychlost u [ms-1]
voda
0.8
vzduch
0.8
voda
-
vzduch
1
objemový zlomek a [1]
Vstup-vzduch
objemový zlomek a [1]
Je možno také sledovat také vliv vztlakové síly na šíření vzduchu.
Matematický model:
Reynoldsovo číslo
Re =400 000 Þ proudění je turbulentní
Výsledky:
obr. 9.3 Izoplochy objemového zlomku vzduchu v podélném a třech příčných řezech, je
uvažována vztlaková síla
obr. 9.4 Izoplochy objemového zlomku vzduchu v podélném a třech příčných řezech, není
uvažována vztlaková síla
96
Vícefázové modely
9.2.2. Trajektorie pevných částic v plynu
Modelování pohybu částic je jedním z typů vícefázového proudění, kdy částice jsou
unášeny proudem v Lagrangeově pojetí. Tato druhá fáze obsahuje kulové částice (také
kapky nebo bubliny). Částice nejsou modelovány jednotlivě, ale jsou definovány vzorky,
charakterizující dostatečně proudění. Je možno řešit následující varianty:
·
výpočet trajektorie diskrétní fáze užitím Lagrangeovy formulace, přitom se zahrnuje
diskrétní fáze pevných částic, hydrodynamický odpor a gravitační síly pro stacionární
i nestacionární případ
·
vliv turbulentních vírů na rozložení částic
·
ohřívání a ochlazování částic
·
hoření částic
·
vznik kapek
Celkový průtok částic je modelován tak, že se sleduje malý počet částic pohybujících
se ve spojité fázi. Částice se pohybuje z bodu, ze kterého je vypuštěna (injection), až do
opuštění oblasti.
Pohyb
částic
(prachových částic,
kapek, bublin)
je
ovlivňován
hydrodynamickým odporem a gravitací. Pohyb částic a zároveň přestup tepla a hmoty mezi
nimi je definován soustavou obyčejných diferenciálních rovnic dle času, obsahující rovnice
pro pohyb, rychlost, teplotu a hmotnost částic. Tyto rovnice jsou integrovány relativně
jednoduchými integračními metodami, aby se vypočetlo chování částic napříč oblastí.
Částice i spojitá fáze se mohou vzájemně ovlivňovat.
Rovnice pohybu částic
Rovnováha sil při použití Lagrangeova přístupu je dána vztahem
duP
g
= FD (u - uP ) + x (ρP - ρ ) + FX
dt
ρP
kde
( 9.2.12)
Fx je vnější objemová síla
FD (u - up ) je síla hydrodynamického odporu vztažená na jednotku hmotnosti částice.
FD =
18 μ CD Re
ρP DP2 24
u – rychlost kapaliny
u p – rychlost částic, řešená integrací podle času
m - molekulární viskozita kapaliny
r - hustota kapaliny
97
Vícefázové modely
r p – hustota částic
Dp – průměr částic
ρDP uP - u
μ
Re – Reynoldsovo číslo
Re =
CD - koeficient hydrodyn. odporu
CD = a1 +
a2
a
+ 32
Re Re
kde a1,2,3 jsou konstanty pro určitý rozsah Re, definované Morsim a Alexandtem
.
Rozptyl částic v důsledku turbulence
Jestliže proudění je turbulentní, pak se počítá trajektorie částice pomocí střední
rychlosti u . Fluent používá stochastické metody (RANDOM WALK MODEL) k určení
fluktuační složky rychlosti, pak u = u + u (t ) . Výpočet trajektorií umožňuje sledovat vliv
turbulence na jejich rozptyl.
Okrajové podmínky na stěně pro DPM
Pokud se částice setkají se stěnou, mohou nastat následující varianty:
·
odraz částice od stěny s uvažováním pružné kolize (reflected)
·
částice projdou hranicí a již se v řešené oblasti nevyskytují (escape)
·
částice se přichytí na stěně (trapped)
·
částice projde vnitřní zónou jako je radiátor nebo porézní prostředí (pass)
·
částice může klouzat podél stěny (slide)
Definování zdroje diskrétní fáze a počátečních podmínek
Pro definici zdroje částic se zadává počáteční poloha zdroje částic, rychlost, velikost
částic případně jejich teplota (hmotnostní průtok). V nabídce materiálů je nutno definovat
fyzikální vlastnosti materiálu a v nabídce fyzikálních modelů zvážit vliv eroze, Brownova
pohybu a dalších jevů.
V počátečních podmínkách se uvažuje rozložení částic podle velikosti. Je možno
uvažovat částice o konstantním průměru nebo lineární změně průměru. Velmi uživané je
Rosin-Rammlerovo rozdělení, kdy jednotlivé prachové částice nejsou o stejných průměrech
a podle jejich velikostí je lze statisticky rozdělit do určitých tříd. Velikost uhelných prachových
částic může např. vyhovovat následujícímu rozdělení, viz tab. 9.1.
98
Vícefázové modely
PRŮMĚR
HMOTNOSTNÍ ZLOMEK
D
YD
0 - 70
0.05
70 - 100
0.1
100 - 120
0.35
120 - 150
0.3
150 - 180
0.15
180 - 200
0.05
tab. 9.1 Rozdělení částic podle velikostí do tříd
Ve Fluentu je rozdělení velikosti částic definováno pomocí Rosinovy - Rammlerovy
exponenciální závislosti:
æ æ D ön ö
YD = exp× ç - ç ÷ ÷
ç èDø ÷
è
ø
D - průměr částice
YD - hmotnostní zlomek částic
D - střední hodnota (MEAN DIAMETER)
n - koeficient distribuce (SPREAD PARAMETER).
Nejprve se určí hmotnostní zlomek příměsi s průměrem větším než D
PRŮMĚR
HMOTNOSTNÍ ZLOMEK
D(mm)
S PRŮMĚREM > D , YD
0
1
70
0.95
100
0.85
120
0.50
150
0.20
180
0.05
200
0.00
1
0.9
0.8
0.7
YD [1]
0.6
0.5
0.4
0.3
0.2
0.1
0
0
50
100
150
200
D [mm]
Logaritmická křivka se určí metodou nejmenších čtverců nebo zjednodušeným postupem
popsaným dále.
D
lze určit z hodnot
D
za zjednodušujícího předpokladu, že
YD = exp(- 1) = 0.368 . Z grafu hodnotě YD = 0.368 odpovídá přibližně průměr D @134 mm
99
Vícefázové modely
a tato hodnota je shodná se středním průměrem D =134 mm. Po dosazení hodnot YD , D a
D do Rosin-Rammlerova vztahu se získá koeficient distribuce n =4.52.
n=
ln ( - ln MD )
D
ln ( )
D
Dále se zadává maximální a minimální průměr částic a hmotnostní průtok.
Specifikace zdroje částic
Vstup částic je definován jako zdroj (injection) v následujících možných variantách
·
single – pro každou částici je nutné definovat počáteční podmínky
·
group – počáteční podmínky jsou definovány pro skupinu částic
·
kužel
·
plocha
·
atd.
obr. 9.5 Definice single, group a spray zdrojů částic
Pro single zdroj je nutno definovat souřadnice zdroje, souřadnice rychlosti, průměr částic,
hmotnostní průtok.
Pro group zdroj se definuje souřadnice prvního a posledního paprsku, ostatní budou mít
souřadnice rozložené lineárně mezi těmito krajními. Ostatní zadání je podobné předchozí
úloze. Následující obrázek znázorňuje varianty group zdroje, viz obr. 9.6.
1) Bodový zdroj umístěn na počátku oblasti–
2) Bodový zdroj umístěn na počátku oblasti-
prachové částice vystupují ve směru
prachové částice vystupují pod úhlem 90°.
rovnoběžném
j1
j1
jn
jn
obr. 9.6 Varianty group zdroje pro směr vodorovný, podobně je zdroj usořádaný pro směr
svislý.
100
Vícefázové modely
Další speciální možnost zadání geometrie zdroje jsou v manuálu [15].
Příklad
Využijte geometrii z předchozího příkladu, který řeší proudění vzduchu a na vstup vložte
zdroj uhelných prachových částic. Zdroj je typu group - úsečka. Pohyb hmotných částic je
ovlivněn vztlakovou silou. Vyhodnotťe trajektorie prachových částic.
y
u
Oblast proudění
z
Vstup - injection
x
obr. 9.7 Schéma řešené oblasti
Geometrie oblasti:
délka oblasti x [m]
3.5
výška oblasti y [m]
0.5
šířka oblasti z [m]
1.5
Fyzikální vlastnosti:
Hustota ρ [kgm-3]
-1
Dynamická viskozita h [kg(ms) ]
Vzduch
Částice
1.225
2000
1.7894e-05
Okrajové podmínky:
Vstup - vzduch
Vstup-injection
intenzita turbulence I [%]
2
hydraulický průměr d h [m]
0.4
rychlost u [ms-1]
1
částice (DPM)
escape
typ
group
počet proudů částic
10
materiál
uhlík (carbon)
rozložení velikosti částic
Rosin-Rammler
101
Vícefázové modely
1. proud
10. proud
x-souřadnice [m]
0
0
y-souřadnice [m]
0.25
0.25
z-souřadnice [m]
0.5
1
-1
0
0
-1
y-ová rychlost [ms ]
0
0
z-ová rychlost [ms-1]
0
0
celkový hmot. průtok [ms-1]
0.0005145
min. průměr [m]
0.000001
max. průměr [m]
0.0001
střední průměr [m]
0.00001
spread parametr [1]
3.5
x-ová rychlost [ms ]
Je možno také sledovat také vliv vztlakové síly na šíření částic.
Matematický model:
Reynoldsovo číslo
Re =23 529 Þ proudění je turbulentní
Výsledky:
obr. 9.8 Trajektorie částic bez uvařování vztlakových sil a s uvažováním vztlakových sil
obr. 9.9 Trajektorie prachových částic s přihlédnutím k turbulentním fluktuacím
102
Časově závislé řešení
10. Časově závislé řešení
Pokud má úloha časově závislý charakter proudění generovaný především vytvořením
sekundárního proudění (vírové cesty) při obtékání těles, viz obr. 10.1, nebo jako odezva na
časově závislé okrajové podmínky, viz obr. 10.2, pak je nutné podstoupit složitější a časově
náročnější řešení, které je funkcí času. Samozřejmě časová závislost se promítá do jevů
přenosu tepla a chemických reakcí.
0
0.005
0.01
0.015
0.02
0.025
0
1
2
3
t [s]
4
5
6
u-konst
u-period
7
obr. 10.1 Vznik vírové cesty při obtékání válce [22]
u [ms ]
-1
obr. 10.2 Průběh rychlosti jako odezva na konstantní a sinusovou rychlost na vstupu zleva
do oblasti (periodická okrajová podmínka vyvolá periodické proudění v celé oblasti) [25]
103
Časově závislé řešení
10.1.
Diskretizace časově závislé rovnice
V případě časově závislého proudění se předpokládá výchozí bilanční rovnice
(pro jednoduchost jednorozměrný tvar) pro obecnou proměnnou veličinu z ve tvaru
¶z ¶
+
(uz ) = ¶
¶t ¶x
¶x
é ¶z ù
êa z ¶x ú + Sz
ë
û
(10.1.1)
v integrální formě je ve tvaru
¶z
ò ¶t
V
é ¶z ù
dV + ò (uz )dA = ò êa z
ú dA + ò Sz dV
x
¶
ë
û
A
A
V
( 10.1.1)
Výchozí rovnice musí být diskretizována v čase i prostoru. Prostorová diskretizace
pro časově závislé rovnice je shodná se stacionární úlohou. Časová diskretizace
zahrnuje integraci každého členu diferenciální rovnice s časovým krokem Dt .
Integrace časových výrazů je jednoduchá, jak bude uvedeno dále.
Výše uvedená rovnice se zapíše v obecné podobě
¶z
= F (z )
¶t
( 10.1.2)
kde funkce F obsahuje prostorovou diskretizaci. Pokud se na časovou derivaci
použije diferenční aproximace prvního řádu vpřed, pak je diskretizovaná rovnice
daná jako
z n +1 - z n
= F (z )
Dt
( 10.1.3)
a případně diskretizace druhého řádu přesnosti je
3z n +1 - 4z n + z n -1
= F (z )
Dt
kde
z
( 10.1.4)
obecná skalární veličina
n + 1 hodnota v následujícím čase t + Dt
n
hodnota v čase t
n - 1 hodnota v předešlém čase t - Dt
Časová diskretizace výchozí rovnice ( 10.1.1) předpokládá implicitní přístup, tedy
konvektivní, difúzní a zdrojový člen jsou vyhodnoceny v čase t + Dt .
¶z
¶z n +1
n +1
n +1 n + 1
dV
u
dA
+
=
z
a
òV ¶t
òA
òA z ¶x dA + Vò Sz dV
( 10.1.5)
V iteračním schématu jsou všechny rovnice řešeny iteračně pro daný časový krok, až je
dosažena konvergence. Tedy řešení v každém časovém kroku požaduje určitý počet
vnějších iterací, dokud nezkonverguje v každém časovém kroku (odpovídá zkonvergování
104
Časově závislé řešení
stacionární úlohy v každém časovém kroku).
t=t+nDt
Řešení pohybové rovnice
Řešení tlakové korekce
Korekce rychlosti a
vnější iterace
tlakového toku
(20)
Řešení skalárních veličin
(T, k, e)
Konvergence
ne
ano
Další časový krok
n=n+1
obr. 10.3 Schéma řešení při použití segregace řešiče.
Volba časového kroku je problematická. Pokud je časová závislost způsobena
známou okrajovou podmínkou, pak je možné přibližně odhadnout časový krok. V opačném
případě je časová závislost způsobena např. odtrhávajícími se víry za ostrou hranou, pak je
nutné velikost časového kroku testovat na počátku výpočtu a splnit následující požadavky:
·
ideální doporučený počet vnějších iterací v každém časovém kroku je 10-20
·
větší počet iterací znamená velký časový krok
·
menší počet iterací znamená malý časový krok
·
počátek výpočtu realizovat pro relativně malý časový krok Dt a v průběhu výpočtu
postupně zvyšovat.
105
Časově závislé řešení
Vyhodnocování výpočtu s časově závislým krokem je možné při automatickém
ukládání datových souborů příkazem FILE-WRITE- AUTO-SAVE . Jde o pravidelné ukládání
výsledků řešení po určitém počtu časových kroků během výpočtu. Jiná možnost je ukládat
hodnoty vybraných proměnných v určitém místě oblasti během průběhu časově závislého
řešení a tak sledovat jejich změny v čase a posoudit, zda se např. řešení blíží ustálenému
stavu při sledování rozběhu systému do ustáleného stavu. Nejdříve se vytvoří tzv.
monitorovací body v nabídce SURFACE-POINT zadáním přesných souřadnic sledovaného
bodu nebo odhadem myší. V příkazu SOLVE-MONITORS-SURFACE INTEGRALS je pak
možno vybrat daný bod a vyhodnocovanou proměnou. Záznam v závislosti na čase lze
zaznamenávat v souboru a zároveň kreslit do grafu na monitor.
Samozřejmě optimální vyhodnocení je animací, vytvořenou přímo softwarově během
výpočtu.
10.2.
Okrajové podmínky
Časově závislé okrajové podmínky se mohou zadávat dvěma způsoby:
·
pomocí souboru (tabulky) pro definici profilu
·
UDF (User Defined Function) - funkce se definuje pomocí C jazyka, uloží, zkompiluje,
přiřadí v okrajových podmínkách pomocí souboru (tabulky) pro definici profilu
10.2.1.
Tabulka pro časovou okrajovou podmínku
Tabulka se vytvoří v textovém editoru s příponou TXT. Format takové tabulky je následující:
profile-name n_field n_data periodic?
field-name-1 field-name-2 field-name-3 .... field-name-n_field
v-1-1 v-2-1 ... ... ... ... v-n_field-1
v-1-2 v-2-2 ... ... ... ... v-n_field-2
...
v-1-n_data v-2-n_data ... ... ... ... v-n_field-n_data
kde
profile-name
sournný název všech proměnných
n-field
počet proměnných
n-data
počet dat charakterizujících funkční závislost (počet
řádků v tabulce)
periodic?
roven 1 pro periodickou podmínku,
roven 0 pro neperiodickou podmínku
106
Časově závislé řešení
field-name-1
je výhradně použito pro vektor času, jehož hodnoty
musí narůstat
field-name-i
vektory dalších proměnných, závislých na čase
v-1-1 … v-n_field-n_data položky v matici, jejíž sloupce odpovídají časové
závislosti vektorů proměnných
Příklad
Definujte soubor pro tabelovanou závislost rychlosti na čase
Tabulka vstupních hodnot
Zadání tabulky vstupních hodnot
pro Fluent
time
u
sampletabprofile 2 3 0
1
10
time u
2
20
1 10
3
30
2 20
3 30
Příklad
Definujte soubor pro tabelovanou periodickou závislost rychlosti na čase.
Proměnná periodicity se nastaví 1.
n_data definuje počet dat charakterizujících jednu
periodu.
Tabulka vstupních hodnot
Zadání tabulky vstupních hodnot
pro Fluent
time
u
sampletabprofile 2 4 1
0
10
time u
1
20
0 10
2
30
1 20
3
10
2 30
3 10
Všechny veličiny se musejí zadat v jednotkách soustavy SI (při čtení profilu se neprovádí
konverze dat a je nutno použít jen malá písmena pro označení proměnných. Profil se přečte
z textového menu následujícími příkazy:
FILE-READ TRANSIENT TABLE
Je možno požít zkratky (f-rtt). Zadá se jméno souboru i s příponou, na obrazovce se objeví
informace o přečtení souboru. Profil se pak zadá do okrajové podmínky příkazy
DEFINE-BOUNDARY CONDITIONS
107
Časově závislé řešení
10.2.2.
UDF pro okrajovou podmínku
Okrajové podmínky závislé na čase se mohou definovat procedurou v C-jazyku.
Proměnné mají přesně definované označení, které je nutné najít v manuálu [15], tam jsou
také jednoduché příklady.
Příklad
Nadefinujte
x-ovou
složku
rychlosti
na
vstupu
pomocí
sinové
funkce
času
u x (t )= u 0 + A sin(wt ) :
/**********************************************************************
unsteady.c
UDF for specifying a transient velocity profile boundary condition
***********************************************************************/
#include "udf.h"
DEFINE_PROFILE(unsteady_velocity, thread, position)
{
face_t f;
real t = CURRENT_TIME;
begin_f_loop(f, thread)
{
F_PROFILE(f, thread, position) = 10. + sin(7.*t);
}
end_f_loop(f, thread)
}
Soubor se vytvoří jako soubor *.txt
způsobem
pomocí
příkazů
a uloží s příponou C. Zkompiluje se interaktivním
DEFINE-UDF-ITERPRETED-COMPILE.
Pak
se
připojí
v okrajových podmínkách pro danou vstupní hranici.
10.3.
Příklad vyhodnocení časově závislé úlohy
Časově závislá úloha ve srovnání s časově nezávislou (stacionární úlohou) je
mnohem složitější, neboť v každém časovém kroku dochází ke změně proudového pole a
tedy všech sledovaných veličin.
108
Časově závislé řešení
Nejdokonalejší obraz řešení umožňuje animace např. vektoru rychlostí, tlaku a
dalších veličin. To je ale velmi náročné z hlediska hardwarového. Navíc prezentování
výsledku vyžaduje počítač, tedy nehodí se do textových zpráv. Je-li opravdu nutné
prezentovat časovou závislost v textu, je možno vytvořit sérii obrázků tak, že se uloží datové
soubory v předem definovaných časových krocích, pak se vytvoří grafické prezentace a vloží
jako obrázky do textového souboru, viz obr. 10.2.
Z důvodu časové a hardwarové náročnosti se využívá jednodušších prostředků
k vyhodnocení, jako je graf závislosti určité veličiny v předem definovaném bodě, nebo
vyhodnocení střední hodnoty veličiny na ploše.
Příklad
Řešte proudění v oblasti, kde na vstupu je dána rychlost, která se periodicky mění dle
funkční závislosti rychlosti na čase. Vyhodnotťe rychlost a tlak ve vybraných bodech.
y
bod-vstup
bod-vir
bod-schod
u
Oblast proudění
z
x
obr. 10.4 Schéma řešené oblasti
Geometrie oblasti:
délka oblasti x [m]
3.5
výška oblasti y [m]
0.5
šířka oblasti z [m]
1.5
Fyzikální vlastnosti:
Vzduch
Hustota ρ [kgm-3]
1.225
Dynamická viskozita h [kg(ms)-1]
1.7894e-05
Okrajové podmínky:
Vstup
rychlost u [ms-1]
2 + sin(7t )
109
Body
vyhodnocení
Časově závislé řešení
střední rychlost us [ms-1]
2
intenzita turbulence I [%]
2
hydraulický průměr d h [m]
0.4
Matematický model:
Re =235 294 Þ proudění je turbulentní
Reynoldsovo číslo
Výsledky:
Je nutno nadefinovat body, ve kterých se provede grafický i textový záznam průběhu
rychlosti a tlaku. Příkazy
Definice bodu:
SURFACE-POINT
zadání souřadnic a název
BOD-VSTUP
BOD-SCHOD
BOD-VSTUP
Definice záznamu:
SOLVE-MONITOR-SURFACE
jméno, plot, write,time step, define
DEFINE-AREA WEIGHTED AVERAGE-FLOW TIME-PRESSURE
výběr bodu BOD-VSTUP
SOLVE-MONITOR-SURFACE
jméno, plot, write,time step, define
DEFINE-AREA WEIGHTED AVERAGE-FLOW TIME-VELOCITY
výběr bodu BOD-SCHOD
SOLVE-MONITOR-SURFACE
jméno, plot, write,time step, define
DEFINE-AREA WEIGHTED AVERAGE-FLOW TIME-VELOCITY
výběr bodu BOD-VIR
Určí se perioda závislosti vstupní rychlosti
w = 2pf = 7 Þ f =
7
w
=
= 1.115 Þ T = 0.897
2p 2.3.14
Spustí se časově závislý výpočet, odhadne se časový krok (menší než desetina periody)
Dt = 0.02 s
Při výpočtu se kontroluje, zda je počet vnitřních iterací menší než 20, jinak se časový krok
bude korigovat. Výsledek výpočtu se zapisuje do souborů BOD-VSTUP.OUT, BODSCHOD.OUT A BOD-VIR.OUT. Soubory jsou textové a přečtou se do EXCELu a připraví se
grafy. Na záznamu reziduálů je vidět periodičnost děje kromě prvních několika iterací, které
jsou ovlivněny tím, že výpočet začíná od počáteční aproximace, která je dána nulovými
hodnotami proměnných, viz obr. 10.5
110
Časově závislé řešení
obr. 10.5 Reziduály periodického řešení
obr. 10.6 Okamžité hodnoty tlaku
Hodnoty rychlosti na vstupu, v bodě BOD-SCHOD a tlaku v bodě BOD-VSTUP jsou
vyhodnoceny a zobrazeny v EXCELu, viz obr. 10.7.
111
Časově závislé řešení
0
2.5
-0.005
2
1.5
u [m/s]
p-vstup
-0.015
p [Pa]
-0.01
r-vir
r-schod
1
-0.02
0.5
-0.025
vliv poč. podmínek
0
-0.03
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
t [s]
obr. 10.7 Vyhodnocení rychlosti a tlaku v bodě v závislosti na čase
Na obr. 10.7 je patrný periodický průběh rychlostí a tlaku, jejichž perioda je shodná
s periodou vstupní rychlosti určené dříve. Perioda je na počátku deformovaná tím, že
výpočet začíná s nulovými počátečními podmínkami uvnitř proudového pole. Po přibližně 1
s je amplituda všech zobrazovaných funkcí již konstantní.
112
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
11. Měření a výpočet turbulence při obtékání válce
v aerodynamickém tunelu
Úloha obtékání válce [22] je jednou ze základních úloh řešení problematiky
obtékaných těles umístěných v proudovém poli. Válec jako takový může reprezentovat
například komín, pilíř, lano, anténu, anténní věž, stožár či jinou vysokou a štíhlou konstrukci
válcového tvaru. Při obtékání takovýchto symetrických těles dochází při určité rychlosti
proudění k odtržení proudící tekutiny od povrchu válce. Toto odtržení je doprovázeno
vznikem vírových struktur, které mají významný vliv jak na samotné obtékané těleso, tak také
například na tělesa nacházející se v jeho úplavu. Vzhledem k tomu, že tato úloha je již
dostatečně zpracována v literatuře, je dobře realizovatelná v laboratorních podmínkách
katedry hydromechaniky a hydraulických zařízení a také ji samozřejmě lze namodelovat
numericky, byla zvolena jako první testovací příklad pro modernizovaný aerodynamický
tunel. Úkolem bylo srovnání fyzikálního experimentu tj. měření jednokanálovým žárovým
anemometrem CTA a matematického modelu vytvořeného v programu GAMBIT a řešeného
v programu FLUENT.
11.1.
Teoretický rozbor úlohy obtékání válce
Při řešení úlohy obtékání válce lze vyhodnotit kromě základních fyzikálních veličin
jako rychlost a tlak a jejich statistického zpracování také Reynoldsovo číslo, Strouhalovo
číslo (frekvenci největších odtrhávajících se vírových struktur), odporové koeficienty, místo
odtržení mezní vrstvy, případně délku úplavu.
Reynoldsovo číslo
Působení proudového pole skutečné (vazké) tekutiny na obtékané těleso je závislé
na hodnotě Reynoldsova čísla. Základní rozdělení charakteru proudění okolo válce při
různých Reynoldsových číslech stanovil experimentálně Roshko. Rozdělil proudění okolo
válce v závislosti na Reynoldsově čísle na následující oblasti:
40 < Re < 150
-
stabilní oblast
150 < Re < 300
-
přechodová oblast
300 < Re < 200 000
-
nestabilní oblast
Podrobnější rozdělení je dosud vzhledem k charakteru turbulence problematické. Další
zkoumané parametry jsou uvedeny v literatuře [22]
11.2.
Fyzikální experiment
Po nastudování metodiky měření žárového anemometru a zvládnutí měřicího zařízení
MiniCTA výrobce DANTEC, viz [28], je možno začít provádět experimentální měření. Pro
113
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
všechna měření byla použita jednodrátková sonda typu 55P14 výrobce DANTEC viz obr.
11.1. Pro vyhodnocení proudového pole byly zvoleny následující parametry proudového
pole:
-
profil střední rychlosti umean
-
profil směrodatných odchylek rychlosti urms
-
profil intenzity turbulence urms umean
-
výkonové spektrum
obr. 11.1 Sonda v měřicím prostoru za
obr. 11.2 Traverzovací zařízení
válcem
Před měřením samotného obtékání válce je nutno provést kalibrací ([28]) žhaveného
drátku sondy 55P14 v rozsahu rychlostí 1 až 20 ms-1. Po vytvoření kalibrační křivky byla
sonda z kalibračního zařízení přemístěna do měřícího prostoru aerodynamického tunelu a
připravena na měření zvolených veličin. V tomto okamžiku je ještě nutné zvolit vzorkovací
frekvenci a dobu měření v jednom bodě. Vzhledem k tomu, že frekvence turbulentního
proudění se pohybuje v rozmezí do fmax = 10 kHz, což je Nyquistova frekvence, je zvolena dle
Shannonova-Kotelnikova kriteria fvz =
1
= 2fmax vzorkovací frekvence pro měření 20 kHz
Dt
[29]. Tato frekvence byla zvolena jako vzorkovací a nastavena v obslužném programu
MiniCTA. Dalším parametrem je volba počtu vzorků. Ten byl zvolen 32 768 v jednom bodě (z
důvodů časových a zejména hardwarových). Výsledná doba měření datového záznamu
v jednom bodě proudového pole Tcelková je tedy:
vzorkovací perioda
T =
1
1
=
= 5 × 10- 5 ,
f 20000
celková doba měření Tcelková = T × n = 5 × 10 -5 × 32768 = 1,64s
kde n = počet vzorků. Na základě zkušebních měření bylo zvoleno, že rychlostní profil bude
obsahovat 160 bodů a vzdálenost mezi jednotlivými body bude 0,5 mm, tj, rychlostní profil
114
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
bude proměřen v šíři y = 80 mm. Pro tyto rozměry byla v programu MiniCTA vytvořena
měřící síť a dokončeny tak všechny přípravy na samotné měření.
x/D=1,2
x/D=5
x/D=2
y
obr. 11.3 Vyobrazení míst pro měření profilů rychlosti a intenzity turbulence
Při měření je posouvání sondy mezi jednotlivými měřenými body prováděno
manuálně pomocí traverzovacího zařízení, viz obr. 11.2. Takto byly proměřeny rychlostní
profily ve třech vzdálenostech za měřeným objektem a to x/D = 1,25, x/D = 2,5 a x/D = 5, viz
obr. 11.3 , kde x je vzdálenost za válcem a D je průměr válce. Výsledky měření profilu
střední rychlosti tj. umean a profilu intenzity turbulence tj. urms / umean jsou uvedeny na obr.
11.4 a obr. 11.5, kde je příklad naměřeného profilu ve vzdálenosti x/D = 2,5 za válcem.
Profil rychlosti za válcem
profil x/D = 1,25
profil x/D = 2,5
profil x/D = 5
u m ean/u vstup [1]
1,2
1
0,8
0,6
0,4
0,2
0
-2
-1,5
-1
-0,5
0
0,5
y/D [1]
obr. 11.4 Profil střední rychlosti
115
1
1,5
2
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
Profil intenzity za válcem
profil x/D = 2,5
profil x/D = 1,25
profil x/D = 5
u rm s/ u m ean [1]
1
0,8
0,6
0,4
0,2
0
-2
-1,5
-1
-0,5
0
0,5
1
1,5
2
y/D [1]
obr. 11.5 Profil intenzity turbulence
Rychlost je zobrazena jako podíl střední naměřené rychlosti v bodě a vstupní
rychlosti. Intenzita turbulence je zobrazena jako podíl směrodatné odchylky v měřeném bodě
a střední rychlosti v tomtéž bodě. Tato volba vyhodnocení byla zvolena z důvodu
maximálního zobecnění naměřených dat.
Naměřený profil střední rychlosti je dle předpokladů symetrický. V místě za válcem
(v oblasti -1 až +1) je zřetelný pokles rychlosti. V tomto místě jsou také patrné mírné
odchylky jednotlivých naměřených bodů od trendu křivky. Tyto odchylky si lze vysvětlit
zvýšením zavíření a vzniku velkých vírových struktur v úplavu za válcem. Jejich odstranění a
vyhlazení naměřených profilů by bylo možné pomocí prodloužení doby měření v jednom
bodě. Zde však vyvstal problém s hardwarovým omezením PC určeného pro měření, na
kterém nebylo možno vyhodnotit potřebné množství dat z jednoho bodu tak, aby byl profil
vyhlazen.
Z vykreslení profilu intenzity turbulence je patrný nárůst intenzity turbulence v úplavu
za válcem. Ze srovnání profilů a to jak rychlosti, tak intenzity turbulence je dále patrné
rozšiřování úplavu. Tento jev je dán disipativním charakterem turbulence. Tzn. velké vírové
struktury vytvořené bezprostředně za obtékaným tělesem se postupně rozpadají a předávají
svou energii dalším částicím unášeným v proudu, až do jejich úplného rozpadu. Právě tento
proces je spojen s rozšiřováním úplavu.
Kromě profilů proudového pole lze pomocí žárového anemometru vyhodnotit
frekvenci odtrhávajících se vírových struktur za obtékaným tělesem a to pomocí vyhodnocení
výkonového spektra, které je v tomto případě získáno ze záznamu hodnot okamžité rychlosti
závislé na čase užitím metody FFT [29], [15]. Tato časová řada byla měřena a po té
vyhodnocena v bodě X = [25mm,10mm] od osy válce ve směru proudění, viz schéma na obr.
11.6.
116
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
9
8
7
6
signal [V]
5
4
3
2
1
0
Spektrální výkonová hustota
-1
Výkonové spektrum
-2
0
P [m 2.s -2]
1
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
t [s]
0.1
0.01
0.001
0.1
1
10
100
1000
f [Hz ]
obr. 11.6 Schématicky znázorněný postup k získání výkonového spektra
V grafu na obr. 11.6 je patrná špička na hodnotě 105 Hz. Tato špička ve výkonovém
spektru indikuje zvýšené množství energie pulzující na této frekvenci, tzn. při této frekvenci
dochází ke zvýšenému množství přenosu energie v úplavu za obtékaným tělesem, tj.
zvýšenému výskytu vírů v této oblasti, což indikuje vírovou cestu. V okamžiku shody této
frekvence s vlastní frekvencí obtékaného tělesa může dojít k nebezpečnému rozkmitání
tohoto obtékaného tělesa až do možných nevratných změn na tělese. Z grafu tedy vyplývá,
že víry s největší energií se odtrhávají s frekvencí 105 Hz. Z vyčíslení Reynodsova čísla:
Re =
v × d 10 × 20 × 10 -3
=
= 11177
1,7894 × 10 - 5
u
( 11.2.1)
lze předpokládat vývin vírové cesty za obtékaným válcem. Tento předpoklad lze dále ověřit
výpočtem Strouhalova čísla. Jeho hodnota se při vývinu vírové cesty pohybuje v rozmezí
Sh » 0,17 - 0,22 . Po dosazení je:
Sh =
f × d 105 × 20 × 10 -3
=
= 0,2
10
v
( 11.2.2)
Z hodnoty vypočteného Strouhalova čísla pro tento případ lze předpokládat, že byla
skutečně naměřena frekvence odtrhávání vírů ve vírové cestě.
117
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
Výsledky fyzikálního experimentu:
u [m.s-1]
Re [1]
FFT max [Hz]
Sh [1]
10
1,1.104
105
0,21
tab. 11.1 Přehled výsledků fyzikálního experimentu
Pro potřeby matematického modelu byly z fyzikálního experimentu zjištěny tyto
parametry:
·
velikost měřicí sekce:
170x140x270 mm
·
průměr vloženého válce:
20 mm
·
rychlost vzduchu:
10 ms-1
·
okolní teplota vzduchu:
22 oC
·
hustota vzduchu:
1,225 kgm-3
·
viskozita vzduchu:
1,7894.10-5 Pa.s
Pro následnou verifikaci výsledků s matematickým modelem je ještě zapotřebí
vyhodnotit možnou chybu fyzikálního experimentu. Výsledky měření jednodrátkovou sondou
mohou být jak mírně nadhodnoceny, tak mírně podhodnoceny. K nadhodnocení dochází
z důvodů ochlazovacího účinku složek rychlosti ve směru y a z. Tento ochlazovací účinek je
sčítán s ochlazovacím účinkem složky x a tak jej nadhodnocuje. Instalovaný suport
umožňuje natočení sondy pouze v jednom směru, což umožňuje měření pouze dvou složek
rychlosti.
Ochlazovací
účinek složek
rychlosti
y
a
z proto
nebyl
vyhodnocován.
Nadhodnocení, které touto cestou vzniká, bylo s největší pravděpodobností kompenzováno
umístěním sondy do měřicího prostoru. Důvodem bylo obtížné stanovení kolmosti
žhaveného drátku sondy k proudu tekutiny (vzduchu). Sonda byla vždy do suportu
umísťována ručně a ručně byl také regulován úhel svírající žhavený drátek s tekutinou.
V ideálním případě by tento úhel měl být 90o. Tato hodnota však nešla ověřit jinak než
vizuální kontrolou. Mírné odchylky od hodnoty 90o však způsobují přesně opačný efekt, než
způsobují složky rychlosti y a z, a to zmenšení průmětu plochy žhaveného drátku v rovině
kolmé na proud vzduchu tj. snížení ochlazovacího efektu dominantní složky rychlosti ve
směru x. Ověření těchto hypotéz však bylo problematické, proto byla odhadnuta chyba
měření na 5 %.
11.3.
Matematický model
Vzhledem k tomu, že výsledky matematického modelu velmi závisí na vytvořené síti
oblasti, bylo vytvořeno v programu Gambit 2.2 a 2.3 několik testovacích verzí sítě. I přes
118
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
možnost volby trojúhelníkové sítě byla nakonec zvolena síť čtyřúhelníková a to z důvodu
lepší konvergence při výpočtu. Hlavní problém při tvorbě sítě nastal u síťování přechodu
okolo válce, kde docházelo k tvorbě jak příliš zdeformovaných, tak příliš velkých buněk. Toto
se nakonec podařilo odstranit rozdělením oblasti okolo válce na dílčí podoblasti společně
s ručním zadáním uzlů sítě na hranicích těchto podoblastí. Nejlepší výsledky výpočtu byly
dosaženy na sítí viz obr. 11.7. Další poznatky z tvorby a testování různých sítí pro tuto oblast
jsou následující:
·
síť musí být symetrická vzhledem k ose válce ve směru proudění, není-li dodržena
tato podmínka může dojít k nesouměrnosti vírové cesty.
·
okolo stěny válce je vhodné vytvořit mezní vrstvu zhuštěnou směrem k této stěně
vzhledem k zachycení zárodku vírových struktur.
·
síť je vhodné vytvořit v úplavu jemnější a naopak okolo stěn oblasti hrubou. Tímto
způsobem lze ušetřit čas při výpočtu.
obr. 11.7 Obrázek sítě + detail kolem válce
Síť byla nejprve vytvořena jako dvourozměrná a poté z ní byla vytvořena její
trojrozměrná varianta. Její rozměry odpovídaly fyzikálnímu experimentu – s výjimkou výšky
z u 3D úlohy – ta byla zmenšena na 50 mm z důvodu výpočtové náročnosti rozsáhlejší
(vyšší) sítě. Pro kompenzaci tohoto zmenšení byla na spodní a horní části oblasti užita
podmínka symetrie. Další okrajové podmínky byly nastaveny takto: vsup do oblasti – velocity
inlet; výstup z oblasti – pressure outlet; stěna levá – wall; stěna pravá – wall; válec – wall.
Pro prvotní verifikaci s fyzikálním experimentem byla úloha modelována jako
dvojrozměrná. Okrajové a fyzikální podmínky byly nastaveny dle fyzikálního experimentu, tj.
rychlost vzduchu 10 ms-1, okolní teplota vzduchu 22 oC, hustota vzduchu 1,225 kg.m-3,
viskozita vzduchu 1,7894.10-5 Pa.s.
U matematického modelu jsou fyzikální veličiny vyhodnocovány v celé oblasti
proudového pole. U případu časově závislé úlohy (jako je obtékání válce) jsou fyzikální
veličiny vyhodnocovány průběžně v čase. Pro verifikaci s výsledky fyzikálního experimentu
119
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
provedeného pomocí měřicí metody CTA jsou však důležité zejména hodnoty analogické
s měřením, tj. získané jako průměr ze všech zaznamenaných hodnot po celou dobu výpočtu.
Na těchto časově „nezávislých“ záznamech nejsou patrné vírové struktury. To je způsobeno
symetrickým odtrháváním vírů, které v časově nezávislém vyhodnocení zaniká, viz obr. 11.8.
Tento jev je logický a mimo jiné potvrzuje korektnost vytvořené sítě a vizuální kontrolu
výpočtu.
obr. 11.8 Střední hodnota velocity magnitude 2D RNG-kε modelu [m.s-1]
Zdali je výpočet již „stabilizovaný“ lze také určit pomocí monitorovaných fyzikálních
veličin. Ve všech případech matematických modelů byla monitorována rychlost v úplavu za
válcem. Na obr. 11.9 je uveden příklad takovéhoto „stabilizovaného“ záznamu pořízeného
v monitorovacím bodě umístěném v úplavu za válcem. Je na něm patrný záznam pravidelné
amplitudy a frekvence rychlosti ve směru proudění.
obr. 11.9 Příklad záznamu rychlosti z monitorovacího bodu
Při matematickém modelování úlohy obtékání válce byly zaznamenávány a
vyhodnocovány profily rychlosti a intenzity turbulence a další a průběh rychlosti v čase
zaznamenávaný v monitorovacím bodě X = [25mm,10mm] (odpovídá záznamu z fyzikálního
experimentu). Z grafu z této časové řady bylo metodou FFT vyhodnocováno výkonové
spektrum, ze kterého lze určit při jaké frekvenci dochází k nejintenzivnějšímu odtrhávání
vírových struktur a z této frekvence lze posléze vypočítat Strouhalovo číslo.
Všechny výpočty matematického modelu proběhly s časovým krokem 1.10-3 s 30ti
120
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
maximálně možnými iteracemi v jednom časovém kroku. Toto nastavení bylo zvoleno jako
minimální s ohledem na délku a přesnost výpočtu – tj. množství časových kroků v jedné
periodě. Frekvence odtrhávajících se vírů byla určena jak z fyzikálního experimentu, tak
z předběžných matematických modelů. Lze tak určit, kolik časových kroků bude počítáno
během jedné periody odtržení víru T =
1
1
=
= 0.095 a počet časových kroků je pak
f 105
pcs =
T
0.095
=
= 9.5 , kde pcs je počet časových kroků a vck je velikost časového
vck 0.001
kroku.
Vzhledem k hustotě sítě nebylo dosaženo v jednotlivých časových krocích
konvergence výpočtu, přičemž konvergenční kritérium bylo nastaveno u všech veličin na
hodnotu 1.10-4.
Pro srovnání jsou v následující kapitole uvedeny také výsledky trojrozměrného
matematického modelu LES pro hodnotu rychlosti 10 m.s-1 a hodnotě časového kroku 1.10-4
s 40ti iteracemi v jednom časovém kroku – čemuž odpovídá 95 časových kroků při jedné
periodě odtržení vírové struktury. Tyto výsledky jsou také vyhodnoceny tabulkami.
11.4.
Vyhodnocení
výsledků
fyzikálního
matematických modelů
Vizuální srovnání výsledků matematických modelů
a) 2D model
RNG-kε
b) 3D model
RNG-kε
c) 3D model
LES
časový krok
1.10-3 [s]
121
experimentu
a
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
u=
d) 3D model
LES
ux
uvstupni
časový krok
1.10-4 [s]
obr. 11.10 Vyobrazení okamžitých hodnot podílu rychlosti tekutiny ve směru x a rychlosti
tekutiny na vstupu do oblasti
Z vizuálního porovnání výsledků jednotlivých modelů na obr. 11.10 je patrné, že při
použití jak 2D modelu, tak 3D modelů vzniká za obtékaným válcem úplav s odtrhávajícími se
vírovými strukturami. Nejintenzivnější odtrhávání je patrné z modelu 3D LES při použití
časového kroku 1.10-4. Ostatní modely vykazují zavíření menší. Je však také patrný rozdíl
mezi modely LES a RNG k-ε. U modelu LES je oblast zpětného proudění za válcem
mnohem zřetelnější.
Další zajímavé, co se týče volby 3D modelování, je srovnání matematických modelů
turbulence při řezu válce osou, rovinou z. Z pohledu na obr. 11.11 je vidět, že při použití
matematického modelu RNG k-ε nejsou rozvinuty turbulentní struktury v ose z, což by podle
teoretických předpokladů mělo nastat. U modelu LES je vývin těchto struktur patrný, přičemž
nejintenzivnější zavíření je opět identifikováno u modelu LES při použití časového kroku
1.10-4.
a) 3D model
RNG-kε
b) 3D model
LES
časový krok
1.10-3 [s]
c) 3D model
LES
časový krok
u=
1.10-4 [s]
ux
uvstupni
obr. 11.11 Okamžité hodnoty podílu rychlosti tekutiny ve směru x a rychlosti tekutiny na
vstupu do oblasti v řezu osou válce
Z tohoto srovnání vyplývá, že vhodnějším modelem pro 3D úlohy bude model LES, neboť
122
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
dvourozměrné modely tohoto případu mají v řešení určitou odchylku od skutečnosti díky
chybějící složce rychlosti v rovině z. To, že u modelu 3D RNG k-ε dochází k menšímu
zpětnému proudění a tudíž i zavíření za válcem je také patrné z vyobrazení vířivosti v ose
z na obr. 11.12.
a) 3D model
RNG-kε
b) 3D model
LES
časový krok
1.10-3 [s]
c) 3D model
u=
LES
ux
uvstupni
časový krok
1.10-4 [s]
+1500 [1/s]
-1500 [1/s]
obr. 11.12 Okamžité hodnoty podílu rychlosti tekutiny ve směru x a rychlosti tekutiny na
vstupu do oblasti a vířivosti v ose z
Porovnání
profilů
rychlosti
a
intenzit
turbulence
u
experimentu
a
matematických modelů
Srovnání profilu rychlosti a intenzity turbulence mezi fyzikálním experimentem a
matematickým modelem je možné pomocí středních hodnot, tj. střední hodnoty rychlosti
v čase a střední hodnoty směrodatné odchylky rychlosti v čase. U fyzikálního experimentu
jsou tyto hodnoty získány z naměřených dat funkcí programu MiniCTA – Data Reduction. Při
této operaci je automaticky proveden přepočet naměřených hodnot napětí přes kalibrační
křivku na hodnoty rychlosti.
U matematického modelu je pro získání těchto veličin nutné zapnout funkci „Data
sampling for time statistics“, pomocí které jsou zaznamenávány v celé oblasti proudového
123
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
pole modelu hodnoty vybraných fyzikálních veličin, tj. hodnoty střední rychlosti ve všech
směrech a jejich směrodatných odchylek a také hodnoty velikosti rychlosti (velocity
magnitude) a jejich směrodatných odchylek.
Pro srovnání byla vybrána vzdálenost x/D = 2.5 v úplavu za válcem. Srovnání je
provedeno na obr. 11.14 a obr. 11.15.
D
obr. 11.13 Schématické zobrazení měření
Profil rychlosti ve vzdálenosti x/D =2,5
u mean/u vstup [1]
CTA
2D RNG-kE
3D LES 1.10-3
3D RNG-kE
3D LES 1.10-4
2
1,8
1,6
1,4
1,2
1
0,8
0,6
0,4
0,2
0
-2
-1,5
-1
-0,5
0
0,5
1
1,5
2
y/D [1]
obr. 11.14 Profil střední rychlosti
Profil intenzity turbulence x/D = 2,5
CTA
2D RNG-kE
3D LES 1.10-3
3D RNG-kE
3D LES 1.10-4
0,9
u rms/ u mean [1]
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
-2
-1,5
-1
-0,5
0
0,5
1
1,5
2
y/D [1]
obr. 11.15 Profil intenzity turbulence
Ze srovnání profilů rychlosti je patrné mírné nadhodnocení výsledků 3D modelů
124
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
s časovým krokem 1.10-3. Toto nadhodnocení může být způsobeno krátkou dobou výpočtu.
3D model s časovým krokem 1.10-4 a 2D model vykazuje velmi dobrou shodu s fyzikálním
experimentem. Mírné podhodnocení může být v případě 2D modelu způsobeno chybějící
složkou rychlosti v rovině z.
Ze srovnání profilů rychlosti je zřetelné výrazné nadhodnocení intenzity turbulence u
matematického modelu LES. To může být způsobeno opět krátkou dobou výpočtu, také příliš
brzkým zapisováním hodnot rychlosti pro časové středování, ale také samotným modelem
LES. U 3D modelu RNG-kε totiž dochází naopak k podhodnocení výsledků intenzity
turbulence. Z vizuálního srovnání modelů viz obr. 11.10 a obr. 11.11 vyplývá, že model
RNG-kε nezachytává všechny vírové struktury. Proto lze také podhodnocení profilu intenzity
turbulence tohoto modelu očekávat. Je zajímavé, že k takovémuto poklesu nedochází také u
2D modelu. Je také otázkou, jak by vypadal tento profil, kdyby byl použit časový krok 1.10-4.
To z časových důvodů nebylo realizováno.
Vzájemné srovnání výkonového spektra experimentu a modelů
Ze
vzájemného
srovnání
výkonového
spektra
fyzikálního
experimentu
a
matematického modelu lze vyčíst, zda se jednotlivé metody liší ve frekvenci, při které
dochází k nejintenzivnějšímu odtrhávání vírových struktur.
F = 105
a) měření
CTA
F = 105
b) 2D model
RNG-kε
125
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
F = 105
c) 3D model
LES
časový
krok
1.10-3 [s]
d) 3D model
F = 105
LES
časový
krok
1.10-4 [s]
interval dat
0,12-0,24
[s]
obr. 11.16 Srovnání průběhu výkonového spektra fyzikálního experimentu a
matematického modelu
Při srovnání výkonových spekter obr. 11.16 je na první pohled patrné, že ač se
jednotlivé průběhy přenosu energie vírovými strukturami liší, frekvence, při které dochází
k nejintenzivnějšímu přenosu energie a tudíž v důsledku k odtrhávání největších vírových
struktur se shoduje, a to na hodnotě 105 Hz.
U 2D-RNG kε modelu je výkonové spektrum výrazně ovlivněno chybějící třetí složkou
rychlosti v rovině z. Průběh odtrhávání vírů je tudíž „zidealizovaný“ a jsou tak patrné vyšší
harmonické frekvence nejintenzivnějšího odtrhávání.
U 3D LES modelu je patrná poměrně dobrá shoda s fyzikálním experimentem,
přičemž hlavní rozdíl je v příkřejším sklonu křivky spektra. Tento jev je způsoben rychlejším
poklesem přenosu energie směrem k vyšším frekvencím, což je způsobeno zanedbáním
vírových struktur, a to jak zjednodušující definicí matematického modelu, tak také díky
velikosti nejmenších buněk sítě, které jsou větší než nejmenší vírové struktury a nejsou tudíž
v modelu zachyceny. U fyzikálního experimentu se tedy při vyšších frekvencích přenáší více
energie než u matematických modelů. Tyto hodnoty však mohou být také mírně
nahodnocené vlivem zašumění vstupného signálu senzoru žárového anemometru CTA.
Ze vzájemného srovnání spekter 3D LES modelu s časovým krokem 1.10-3 a 1.10-4
-3
patrné, že model s časovým krokem 1.10
je
má zřetelnější maximum u modelu s časovým
126
Měření a výpočet turbulence při obtékání válce v aerodynamickém tunelu
krokem 1.10-4
toto maximum ještě není tak patrné. Je jasné, že by tento model potřeboval
ještě čas na úplné ustálení proudění. Což je zásadní nevýhodou modelu s tímto časovým
krokem – doba výpočtu.
Závěrem z tohoto srovnání lze říci, že je-li potřeba zjistit frekvence, při kterých
dochází k největším přenosům energie a není-li zapotřebí příliš klást důraz na co
nejpřesnější výpočet samotného proudění, je z časových důvodů u úlohy obtékání jednoho
válce výhodnější použití 2D modelu. Ze srovnání výkonových spekter je ještě parné, že
největší shodu vykazuje model 3D LES s časovým krokem 1.10-4. Je však třeba konstatovat,
že výkonové spektrum tohoto modelu bylo provedeno z malého množství dat (časový interval
0,12-0,14 s) tato úloha nebyla dále počítána z důvodů její časové náročnosti. Proto byl také
upřednostněn pro další modely časový krok 1.10-3.
127
Literatura
12. Literatura
[1]
REYNOLDS, O.: An Experimental Investigation of the Circumstances which
determine whether the Motion of Water shall be direct or sinuous, and the Law of Resistance
in Parallel Channels. Phil. Trans. Roy. Soc., 1883, Vol. 174, p. 935.
[2]
KOZUBKOVÁ, M., DRÁBKOVÁ, S., JANČÍK, P. , SKÝBA, V. Závěrečná zpráva o
řešení projektu COST OC 715.60 COST Action 715: „Meteorology applied to urban air
pollution problems“, 715.60 „Numerical Modelling of the Small Scale Urban Air Flow and
Pollutant Dispersion“. Ostrava: VŠB-TUO, 2004, 23 s.
[3]
LESIEUR, M. Turbulence in Fluids. Third Revised and Enlargened Version. Kluwer
Academic Publishers, Dordrecht, 1990, ISBN 0-7923-4416-2, 395 p.
[4]
CIOFALO, M. Large Eddy Simulation of Turbulent Flows with Heat. A State of the Art
Review. Palermo: Dipartemento di Ingegneria Nucleare-Universita di Palermo, 1993, 95 p.
[5]
VAN DYKE, M.: Albom těčenij židkosti i gaza. Moskva: MIR, 1986, 180 p.
[6]
NIEWSTADT, F.T.M. et al. Direct and Large Eddy simulations of turbulence in fluids.
Future Generation Computer System, 1994, p.189-205.
[7]
DVOŘÁK, R.: Vnitřní aerodynamika, ČVUT Praha 1997
[8]
RODI, W., MURAKAMI, S. Turbulence Models for Practical Applications, výtah
z „Seisan-Kenkyu“, sv. 41, č.8, 9, 12, sv.42, č. 1,2, Institute of Industrial Science, University
of Tokyo 1989.
[9]
ZUBER, I. Model turbulence k-e kontra Reynolds stress model (jedna pověra
v hydrodynamice. Osobní sdělení
[10]
ŠŤÁVA, P., JANALÍK, J., KOZUBKOVÁ, M.DRÁBKOVÁ, S.: Aspekty matematického
modelování turbulentníbo proudění, Seminář 3 , Ostrava 1994
[11]
RODI,W.,
LESCHZINER,M.A.,
SCHEUERER,G.,
SCHONUNG,B.:
Numerische
Berechnung turbulenter Stroemungen in Forschung und Praxis, Sonderforschungsbereich
210, TU Karlsruhe,1992
[12]
STULL, R.B. An Introduction to Boundary Layer Meteorology, Kluwer Academic
Publishers, 1994, 666 str.
[13]
TESAŘ,V.: Mezní vrstvy a turbulence, skripta ČVUT Praha 1991
[14]
COANTIC, M. Turbulent Fluid Flow, An Introduction to its Physics and Modelling. In A
series of lecture notes for the seminar „Industrial Mathematics and Mathematical Modelling“.
Plzeň: University of Western Bohemia, 1994, 26 p.
[15]
FLUENT: Fluent 6.3.26 - User’s guide Fluent Inc. 2006. VŠB-TU Ostrava
<URL:http://sp1.vsb.cz/DOC/Fluent_6.3/html/ug/main_pre.htm>
[16]
ROACHE,P.J.: Computationa Fluid Dynamics, Mir, Moskva, 1980
[17]
SUNG-EUN KIM, DIPANKAR CHOUDHURY: Computations of Complex Turbulent
128
Literatura
Flows and Heat Transfer Using Two-Layer Based Wall Functions, Fluent Inc., 1995
[18]
BIRD, R.,B., STEWART, W., E., LIGHTFOOT, E., N. Přenosové jevy. Academia,
Praha 1968
[19]
ŠŤÁVA, P.: Turbulence a ekologické úlohy, VŠB-TU Ostrava, fakulta strojní, 1998
[20]
KOZUBKOVÁ, M., DRÁBKOVÁ, S. Numerické modelování proudění – FLUENT I.
[Online]. c2003. Ostrava: VŠB – TUO, 116 s, poslední revize 6.1.2005 [cit. 2005-01-06].
Dostupné z: <URL: http://www.338.vsb.cz/seznam.htm>.
[21]
YAKHOT, V., ORSZAG, S. A. Renormalization Group Analysis of Turbulence: I. Basic
Theory. Journal of Scientific Computing, 1(1): p.1-51, 1986.
[22]
FABIÁN, P. Metody matematického a fyzikálního experimentu v proudění tekutin.
Dizertační práce. 2008. Ostrava: VŠB – TUO, 115 s.
[23]
ZAVILA, O. Matematické modelování turbulentního proudění, šíření tepla a polutantu
v mezní vrstvě atmosféry se zaměřením na tunelové stavby. Dizertační práce. 2007.
Ostrava: VŠB – TUO, 114 s.
[24]
MARÁK, R. Analýza pojistného ventilu matematickými metodami. Diplomová práce.
2007. Ostrava: VŠB – TUO, 50 s.
[25]
DRÁBKOVÁ, S. a kol. Mechanika tekutin. Ostrava: VŠB-TU Ostrava, 2007. 248 s.
(Elearningová učebnice). ISBN 978-80-248-1508-4.
[26]
MICHALCOVÁ, V. Numerické modelování zatížení budov při kvazistatickém působení
větru. Dizertační práce. 2007. Ostrava: VŠB – TUO, 104 s.
[27]
Incropera, P. a kol.: Fundamentals of Heat and Mass Transfer. John Wiley&Sons.
2007. 997 p. ISBN-13: 978-0-471-45728-2, ISBN-10: 0-471-45728-0.
[28]
FABIÁN, P., KUNDYS, J., KOZUBKOVÁ, M. Měření žárovým anemometrem [online].
Návody do cvičení, Kat. hydromechaniky a hydraul. zařízení, FS VŠB-TU Ostrava, 2004,
Scripta electronica. 36 s. Dostupné z : <URL http://www.338.vsb.cz/seznam.htm> [cit. 200601-06].
[29]
JANALÍK, J. Měření tekutinových mechanizmů [online], Ostrava: VŠB – TUO, 2003.
116 s. Scripta electronica, Dostupné z: <URL: http://www.338.vsb.cz/seznam.htm> [cit.
2006-01-06].
[30]
BOJKO, Marian; KOVÁŘ Ladislav. Mathematical Modelling of the Pure Oxygen Blowing
Process on Surface of Liquid Bath in Metallurgical Aggregates. Sborník vědeckých prací VŠB-TU
Ostrava, řada strojní, 2007, roč. LIII, č.1, s. 7 – 14. ISBN 978-80-248-1633-3, ISSN 1210-0471.
[31]
MALÝ, Rostislav;
BLEJCHAŘ, Tomáš. Využití CFD při vývoji plazmové technologie.
TRANSFER 2007, 2007, str.109-112. ISBN 978-80-8075-236-1, ISSN 1336-9695.
[32]
LENHARD, R., JANDAČKA, J.
Simulácia parametrov pasívného stropného chladiaceho
konvertora. In Sborník XXVII. setkání kateder mechaniky tekutin a termomechaniky. Plzeň: ZČU
v Plzní. 24.-27. června 2008. 6 str. ISBN 978-80-7043-666-0
129
Literatura
[33]
GAŠPAROVIČ, P., ČARNOGURSKÁ, M. Aerodynamic Optimisation of Centrifugal Fan
Casing using CFD. J. of applied science in the thermodynamics and fluid mechanics. VOl. 2. No.
1/2008, 6 p. ISSN 1802-9388.
130
Příloha – sčítací pravidla
13. Příloha - sčítací pravidla
Při vyjádření proměnných o třech případně devíti složkách (složky rychlostí, napětí
apod.) je vhodné využít speciální zkrácené označení, známé jako Einsteinova sumace, kdy
pouze jedním členem lze vyjádřit všech devět napětí. V reálném případě se očekává, že
například některé síly budou působit ve všech souřadnicových směrech, jiné pouze v
některých. Záměrem je, abychom byli schopni oddělit takové závislosti. V této kapitole budou
nejdříve definovány speciální výrazy, pak pravidla pro počítání s nimi a na závěr příklady a
vazba na vektorové vyjádření [12].
13.1.
Definice a pravidla
Nechť m,n,q jsou celočíselné proměnné nabývající hodnot 1,2 nebo 3. Nechť am
představuje obecně vektor rychlosti, xm představuje souřadnici, dm je jednotkový vektor
(vektor délky 1 a směru shodného s jedním ze tří kartézských souřadnic).
m = 1, 2, 3
a1 = u
x1 = x
n = 1, 2, 3
a2 = v
x2 = y
q = 1, 2, 3
a3 = w
x3 = z
Proměnná a
-bez indexu
= skalár
-s jedním indexem
= vektor
-se dvěma indexy
= tensor
( 13.1.1)
Jednotkový vektor:
e1=i
e2=j
e3=k
( 13.1.2)
Kroneckerovo delta
ì+ 1 pro m = n
ì+ 1 pro m = n
d mn = í
d mn = í
î0 pro m ¹ n
î0 pro m ¹ n
( 13.1.3)
Jednotkový tenzor
e mnq
ì+ 1 pro mnq = 123,231,312
ï
= í- 1 pro mnq = 321,213,132
ï0 pro ostatn’ kombinace
î
( 13.1.4)
Základní pravidla:
1. Jestliže se objevují v zápise dva indexy v témže členu, znamená to, že se jedná o součet
členů přes všechny hodnoty (1,2,3) opakovaného indexu.
131
Příloha – sčítací pravidla
2. Jestliže jeden index se předpokládá jako nesčítací index (volný index), pak se tento index
musí objevit jako nesčítací ve všech členech rovnice. Tudíž taková rovnice stručně
představuje 3 rovnice pro každou hodnotu nesčítacího indexu.
3. Tentýž index se nemůže objevit více než dvakrát v jednom členu.
13.2.
Příklady
Příklad 1: (užití pravidla 1)
an
¶bm
¶b
¶b
¶b
= a1 m + a2 m + a3 m =
¶x n
¶x1
¶x2
¶ x3
¶b
¶b
¶b
= u m +v m +w m
¶x
¶y
¶z
( 13.2.1)
Příklad 2: (užití pravidla 1)
d 2 n an = d 21a1 + d 22a2 + d 23a3 =
=0+a2 +0 =
=v
( 13.2.2)
Příklad 3: Rozepište následující výraz (užití pravidla 2)
am = bm + d mnc n
řešení:
a1 = b1 + d1n cn = b1 + c1
a2 = b2 + d 2 n cn = b2 + c2
( 13.2.3)
a3 = b3 + d 3 n cn = b3 + c3
Příklad 4: Užijte zkrácený zápis pro pohybovou rovnici a rozepište:
¶am
¶a
1 ¶p 1 é ¶t mn ù
+ bn m = -d m 3g + fce mn 3 bn +
¶t
¶x n
r ¶x m r êë ¶xn úû
Nejdříve se rozepíše suma pro sčítací index
¶am
¶a
¶a
¶a
+ b1 m + b2 m + b3 m = -d m 3g + fce m13 b1 + fce m 23 b2 + fce m 33 b3 ¶t
¶x1
¶x2
¶x3
-
1 ¶p 1 é ¶t m1 ¶t m 2 ¶t m 3 ù
+
+
+
r ¶x m r êë ¶x1
¶x2
¶x3 úû
Člen em33 je roven nule, protože obsahuje dva shodné indexy.
Pro m=1:
132
( 13.2.4)
Příloha – sčítací pravidla
¶a1
¶a
¶a
¶a
+ b1 1 + b2 1 + b3 1 = -d13 g + fce 113 b1 + fce 123 b2 ¶t
¶x1
¶x2
¶x3
-
1 ¶p 1 é ¶t 11 ¶t 12 ¶t 13 ù
+
+
+
r ¶x1 r êë ¶x1 ¶x2 ¶x3 úû
V této rovnici je d13 a e113 rovno také nule podle definice.
Pro m=2:
¶a2
¶a
¶a
¶a
+ b1 2 + b2 2 + b3 2 = + fce 213 b1 ¶t
¶x1
¶x2
¶x3
-
1 ¶p 1 é ¶t 21 ¶t 22 ¶t 23 ù
+
+
+
r ¶x2 r êë ¶x1 ¶x 2 ¶x3 úû
Člen e213 je roven -1.
Pro m=3:
¶a3
¶a
¶a
¶a
+ b1 3 + b2 3 + b3 3 = -d 33 g ¶t
¶x1
¶x2
¶x3
-
1 ¶p 1 é ¶t 31 ¶t 32 ¶t 33 ù
+
+
+
r ¶x3 r êë ¶x1 ¶x2 ¶x3 úû
Při užití substituce ( 13.1.1) lze tyto rovnice vyjádřit obvyklejším způsobem:
¶u
¶u
¶u
¶u
1 ¶p 1 é ¶t xx ¶t xy ¶t xz ù
+u
+v
+w
= +fcv +
+
+
¶t
¶x
¶y
¶z
r ¶x r êë ¶x
¶y
¶z úû
1 ¶p 1 é ¶t yx ¶t yy ¶t yz ù
¶v
¶v
¶v
¶v
+u
+v
+w
= -fcu +
+
+
¶t
¶x
¶y
¶z
r ¶y r êë ¶x
¶y
¶z úû
¶w
¶w
¶w
¶w
1 ¶p 1 é ¶t zx ¶t zy ¶t zz ù
+u
+v
+w
= -g +
+
+
¶t
¶x
¶y
¶z
r ¶z r êë ¶x
¶y
¶z úû
( 13.2.5)
Srovnáním zápisu ( 13.2.5) a ( 13.2.4) lze pozorovat obrovskou výhodu při užití
Einsteinových sumačních pravidel a navíc existuje určitá analogie s vektorovým značením,
takže soustavu rovnic ( 13.2.5) lze vyjádřit takto:
¶ui
¶u
1 ¶p 1 é ¶t ij ù
+ u j i = -d i 3 g + fce ij 3u j + ê
ú
¶t
¶x j
r ¶xi r êë ¶x j úû
13.3.
( 13.2.6)
Srovnání s vektorovým označením
Vektory jsou reprezentovány třemi složkami, tenzory devíti složkami. Tudíž existuje
spojitost mezi Einsteinovu sumací a definicí vektoru. Vektory mohou být přepsány v
sumačním označení, stejně jako vektorové operace.
Vektor
u=umem
133
Příloha – sčítací pravidla
Skalární součin
em.en=dmn
Vektorový součin
emxen=emnqdq
Operátor Ñ
Ñ ( ) = em
( 13.3.1)
¶( )
¶xm
Na příkladech budou demonstrovány předchozí zápisy.
· skalární součin
rr
u.v = (emum )(
. env n ) = (em .en )umv n = d mnumv n = umv m
( 13.3.2)
První úprava spočívá v substituci vektorů jejich sumačním zápisem, dále jsou vytknuty
souřadnice vektorů a provedeno skalární pronásobení jednotkových vektorů.
· vektorový součin
r r
u xv = (emu m )x (env n ) = (em xen )u mv n = e mnq u mv n eq
( 13.3.3)
Výsledkem je vektor. Poznamenejme, že emnq=eqmn=enqm
· divergence vektoru
¶ ö
¶u
¶u
¶um
r æ
÷÷.(en un ) = (em .en ) n = d mn n =
Ñu = çç em
¶xm
¶x m ¶x m
è ¶xm ø
13.4.
( 13.3.4)
Substanciální derivace
Pohyb tekutin se definuje v určitém vymezeném dostatečně malém prostoru, jímž
tekutina prochází a nazývá se konečný objem ohraničený plochami. V tomto objemu jsou
odvozeny základní rovnice zachování někdy také nazývané jako rovnice změny, protože
popisují změnu rychlosti, teploty a koncentrace vzhledem k místu a času systému, tj.
všechny hledané veličiny jsou obecně funkcí času a souřadnic z = z (t , x, y , z )
Nejdříve ale je třeba se seznámit se třemi druhy derivace obecné proměnné z podle
času, kterých se bude používat.
·
Parciální derivace podle času
¶z
je změna veličiny z podle času při konstantních
¶t
prostorových souřadnicích
·
Totální derivace podle času
dz ¶z ¶z ¶x ¶z ¶y ¶z ¶z ¶z ¶z ¶x j
=
+
+
+
=
+
dt
¶t ¶x ¶t ¶y ¶t ¶z ¶t
¶t ¶x j ¶t
je
změna veličiny z podle času a souřadnic (Lagrangeovo pojetí)
· substanciální derivace podle času
¶z
¶z
¶z ¶z
¶z
Dz ¶z
=
+u
+v
+w
=
+ uj
¶t
¶x
¶y
¶z
¶t
Dt
¶x j
je
zvláštní druh totální derivace (derivace sledující pohyb), kde u, v , w jsou složky místní
134
Příloha – sčítací pravidla
rychlosti (Eulerovo pojetí).
135
Příloha - Příloha - základní statistické metody
14. Příloha - základní statistické metody
Základním postupem k řešení turbulentního problému je stochastický přístup, proto je
užitečné mít základní znalosti ze statistiky. V následující kapitole se budou definovat
základní statistické pojmy a metody, jako je střední hodnota, rozptyl (variance), standardní
odchylka, kovariance a korelace.
14.1.
Střední hodnota
Pro časové, prostorové a smíšené (časové i prostorové) průměry existují tři způsoby
vyjádření. Časový průměr je definován v určitém bodě prostoru a skládá se ze součtu resp.
integrálu přes časovou periodu T. Pro libovolnou proměnnou u(t,s), která je funkcí času a
prostoru platí:
t
u (s ) =
1 N -1
å u(i , s )
N i =0
T
nebo
t
u (s ) =
1
u(t , s )dt
T t =ò0
( 14.1.1)
kde t=iDt. Prostorový průměr je definován v určitém čase a je dán součtem nebo integrálem
přes prostorovou oblast S:
S
u (s ) =
1 N -1
å u (t , j )
N j =0
S
nebo
S
u (s ) =
1
u(t , s )ds
S s ò= 0
( 14.1.2)
kde s=jDs. Smíšený průměr sestává ze sumy všech experimentů
e
u (t , s ) =
1 N -1
å ui (t, s )
N i =0
( 14.1.3)
obr. 14.1 Homogenní turbulence generovaná laboratorně za mřížkou [5]
Je-li turbulence homogenní (statisticky je shodná v každém bodě prostoru, viz obr. 14.1,
pak každý senzor bude měřit tutéž veličinu při prostorovém zprůměrování. Časové průměry
136
Příloha - Příloha - základní statistické metody
jsou často užívány a počítány ze senzorů umístěných v jediném určitém bodě. U turbulence,
která je homogenní a stacionární (statisticky se nemění v čase), jsou časové, prostorové a
smíšené průměry stejné. Tento jev se nazývá ergodičnost.
e
( ) = t ( ) =S ( ) = ( )
( 14.1.4)
Pravidla pro počítání se střední hodnotou
Nechť a a b jsou dvě proměnné závislé na čase a c je konstanta. Pak platí
(a + b ) = a + b
(c.a ) = c (a )
(a) = a
( 14.1.5)
æ da ö d a
ç ÷=
è dt ø dt
obr. 14.2Schematické srovnání lokální derivace a derivace středované složky.
14.2.
Reynoldsova pravidla
Pravidla o časovém průměrování budou aplikována na proměnné, které lze rozložit
na část časově středovanou a fluktuační složku a = a + a¢ a b = b + b¢ , přičemž platí
T
(a ) = 1 ò adt = (a + a¢) = (a ) + a¢ = a + a¢Þa ¢ = 0
T
0
( 14.2.1)
pro takto rozepsané veličiny platí Reynoldsova pravidla
a=a
a a¢ = (aa¢) = aa¢ = 0
a + b = (a + a¢) + (b + b¢) = (a + b ) + (a¢ + b¢) = (a + b ) + (a¢ + b¢) = a + b
ab = (a + a¢)(
. b + b¢) = (ab + a¢b + ab¢ + a¢b¢) = (ab ) + (a¢b ) + (ab¢) + (a¢b¢) = (ab ) + (a¢b¢)
( 14.2.2)
kde a ¢b ¢ je korelační moment fluktuačních složek, který je obecně nenulový podobně jako
137
Příloha - Příloha - základní statistické metody
a¢ 2 ,a¢b¢ 2 ,a¢ 2 b¢ 2 . Tato vlastnost je typická pro turbulenci, což je odlišné od mnoha lineárních
teorií vlnění, kdy nelineární členy se často zanedbávají.
14.3.
Rozptyl, směrodatná odchylka, turbulentní intenzita,
kovariance, korelace
Významnou statistickou veličinou daného souboru dat při známé střední hodnotě je
rozptyl definovaný jako
s a2 =
1 N -1
(ai - a )2
å
N i =1
( 14.3.1)
Lepším odhadem rozptylu populace daného souboru dat je
1 N -1
s =
(ai - a)2
å
N - 1 i =1
2
a
( 14.3.2)
Je-li počet měření N velký, jak je u turbulence zvykem, pak se oba vzorce shodují.
Poznamenejme, že turbulentní část (fluktuaci) turbulentní proměnné lze vyjádřit jako
a¢ = a - a . Pak po substituci do vzorce ( 14.3.2) je
s a2 =
1 N -1 2
ai¢ = a¢2
å
N - 1 i =1
( 14.3.3)
Tedy určíme-li průměr z druhé mocniny fluktuace jako u¢2, v ¢2 ,w ¢2 ,T ¢2, q¢2 , můžeme jej
interpretovat jako rozptyl.
Směrodatná odchylka (obr. 14.3) je definována jako druhá odmocnina rozptylu:
( )
s a = a¢2
1/ 2
( 14.3.4)
obr. 14.3 Směrodatná odchylka a střední hodnota
Směrodatná odchylka má stejný rozměr jako původní proměnná. Na obr. je ukázán vztah
mezi turbulentní veličinou a směrodatnou odchylkou, která může být vysvětlena jako
vzdálenost původních dat od průměru. Proto je používána k odhadu intenzity turbulence
jako bezrozměrné veličiny definované
138
Příloha - Příloha - základní statistické metody
I=
sa
a
( 14.3.5)
Ve statistice je kovariance mezi dvěma proměnnými definovaná jako
cov ar (a, b ) =
1 N -1
1 N -1
(
a
a
)(
b
b
)
=
å i
å (ai¢ )(bi¢ ) =a¢b¢
i
N i =0
N i =0
( 14.3.6)
Kovariance určuje stupeň závislosti mezi dvěma veličinami. Někdy je užitečné používat
bezrozměrnou kovarianci jako lineární korelační koeficient
rab =
a¢b¢
s as b
( 14.3.7)
Tato veličina je v absolutní hodnotě menší než 1. Jeli rovna 1, pak veličiny jsou korelovány,
jeli rovna -1, pak jsou opačně korelovány. Veličiny, které spolu vůbec nesouvisejí, mají r=0.
139
Příloha - Příloha - rozklad v Taylorovu řadu
15. Příloha - rozklad v Taylorovu řadu
15.1.
Definování problému
Modelovou rovnicí přenosu je linearizovaná jednorozměrná rovnice s konvektivním a
difúzním členem, zapsaná v konzervativním tvaru [16]
¶z
¶ (uz )
¶ 2z
=+a
¶t
¶x
¶ x2
( 15.1.1)
æ ¶ uz
¶u
¶z ö
÷
çç
= z.
+ u.
¶x
¶ x ÷ø
è ¶x
nebo v nekonzervativním tvaru
¶z
¶z
¶ 2z
= -u
+a
¶t
¶x
¶ x2
( 15.1.2)
kde
z - je vírová funkce nebo jakákoliv jiná skalární veličina (T,z )
a - je zobecněný koeficient difúze
u - je linearizovaná rychlost konvekce.
Jestli se neřekne jinak, tak u je konstantní vzhledem k souřadnici x, ačkoliv rovnice může být
užita k zkoumání stability i v případě, že rychlost u je funkcí x, tj. u=u(x).
15.2.
Taylorův rozvoj
Základní diferenční schemata se získají použitím rozkladu v Taylorovu řadu. Použije
se pravoúhlá síť dle obr. 15.1.
t
Dx
n+1
Dt
n
n-1
i-1
i
i+1
x
L
obr. 15.1 Definice sítě
Index i odpovídá Dx , n odpovídá časovému kroku Dt . Předpokládá se, že funkce f a její
derivace jsou spojité do potřebného řádu. Pak Taylorova řada pro směr x je
140
Příloha - Příloha - rozklad v Taylorovu řadu
1 æ ¶ 2f ö
æ ¶f ö
2
fi +1,n = fi ,n + ç ÷ (xi +1,n - xi ,n ) + çç 2 ÷÷ (xi +1,n - xi ,n ) + ...
2 è ¶x øi ,n
è ¶x øi ,n
1 æ ¶ 2f ö
æ ¶f ö
= fi ,n + ç ÷ Dx + çç 2 ÷÷ Dx 2 + 0 Dx 3
2 è ¶x øi ,n
è ¶x øi ,n
( )
( 15.2.1)
eventuelně
1 æ ¶ 2f ö
æ ¶f ö
2
fi -1,n = fi ,n - ç ÷ (x i -1,n - xi ,n ) + çç 2 ÷÷ (x i -1,n - xi ,n ) + ...
2 è ¶x øi ,n
è ¶x øi ,n
1 æ ¶ 2f ö
æ ¶f ö
= fi ,n - ç ÷ Dx + çç 2 ÷÷ Dx 2 + 0 Dx 3
2 è ¶x øi ,n
è ¶x ø i ,n
( )
( 15.2.2)
a je zřejmé, že má o řád vyšší přesnost. Nyní se odvodí diferenční aproximace pro druhou
derivaci sečtením ( 15.2.1) a ( 15.2.2)
æ ¶ 2f ö
fi +1,n + fi -1,n = 2fi ,n + çç 2 ÷÷ Dx 2 + 0( Dx 4 )
è ¶x ø i ,n
Z toho
fi -1,n - fi ,n 1 æ ¶ 2f ö
æ ¶f ö
=
- çç 2 ÷÷ Dx + 0 Dx 2
ç ÷
Dx
2 è ¶x øi ,n
è ¶x øi ,n
( )
( 15.2.3)
nebo méně přesně
fi +1,n - fi ,n
æ ¶f ö
+ 0(Dx )
ç ÷ =
Dx
è ¶x øi ,n
( 15.2.4)
Tedy pro první derivaci při diferenční aproximaci vpřed dle ( 15.2.4) se dostane
fi +1,n - fi ,n
æ ¶f ö
ç ÷ =
Dx
è ¶x ø i ,n
( 15.2.5)
Jestliže se použije Taylorova rozvoje (C.2.2), dostane se pro první derivaci
diferenční aproximace vzad
fi ,n - fi -1,n
æ ¶f ö
ç ÷ =
Dx
è ¶x øi ,n
( 15.2.6)
Tzv. centrální diferenční aproximace první derivace se získá jako rozdíl rozkladů
(C.2.1) a (C.2.2) a
fi +1,n - fi -1,n
æ ¶f ö
+ 0 Dx 2
ç ÷ =
2D x
è ¶x øi ,n
( )
Tedy centrální diference první derivace je
fi +1,n - fi -1,n
æ ¶f ö
ç ÷ =
2Dx
è ¶x øi ,n
( 15.2.7)
141
Příloha - Příloha - rozklad v Taylorovu řadu
( 15.2.8)
æ ¶ 2f ö
f
+f
- 2fi ,n
çç 2 ÷÷ = i +1,n i -1,2n
Dx
è ¶x ø i , n
( )
s řádem přesnosti 0 Dx 2 . Podobný postup platí pro nezávislou proměnnou y a poloviční
krok
Dx Dt
,
atd. Pokud se použije centrální diferenční aproximace na modelovou rovnici
2 2
(C.1.2), získá se
z i ,n +1 - z i ,n -1
z
u.z i +1,n - uz i -1,n
+ z i -1,n - 2z i ,n
=+ a i +1,n
2.Dt
2D x
Dx 2
( 15.2.9)
kde lze explicitně vyjádřit z ni +1 pomocí hodnot na předchozím časovém kroku n. Ale lze
dokázat, že toto schema je pro Dt ñ 0 a a ñ 0 nestabilní, tzn. že přivádí k vzniku chaotických
řešení nemajících vztah k diferenciální rovnici. Jestliže místo centrálních diferencí v
nestacionárním členu se použije diferenční schema vpřed, je:
uz
- uz i -1,n
+ z i -1,n - 2z i ,n
z i , n +1 - z i , n
z
= - i +1,n
+ a i +1,n
Dt
2Dx
Dx 2
(15.2.10)
Takové schema se nazývá schema FTCS (Forward Time Central Space). Lze dokázat, že
toto schema je stabilní alespoň při omezujících podmínkách na Dt, u, a a Dx. Při užití
diferencí vpřed i vzad (upwind) na prostorovou první derivaci lze získat schéma stabilní vždy.
142
Číslo skladové:
2308
Určeno pro posluchače:
1. r. nav. mgr. FS
Autor:
Katedra, institut:
Název:
300
Milada Kozubková, doc. RNDr., CSc.
hydromechaniky a
hydraulických zařízení
338
Modelování proudění tekutin,
FLUENT, CFX
Místo, rok, vydání:
Ostrava, 2008, 1. vydání
Počet stran
153
VŠB – TECHNICKÁ UNIVERZITA
Vydala:
OSTRAVA
17. listopadu 15/2172
708 33 Ostrava-Poruba
Výroba:
katedra hydromechaniky a
hydraulických zařízení
Náklad:
20 CD
Tématická skupina:
17
Download

Kozubková, M. - Katedra hydromechaniky a hydraulických zařízení