NUMERICKÁ ANALÝZA PROCESů
NAP10
Řešení transportních rovnic metodou kontrolních objemů.
Řešení Navierových Stokesových rovnic v primitivních proměnných (metoda SIMPLE) i použitím proudové funkce a vířivosti.
Laminární tok v dutině (FLUENT a vlastní program MATLAB)
Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010
NAP10 Řešení transportních rovnic
( ) [ ( ) ]
( ) ( )
CV CV
A A CV
u dV S dV
n u dA n dA S dV
Asi nejčastěji používanou metodou řešení transportních rovnic je metoda kontrolních objemů (např.Fluent).
Např. pro stacionární případ a transportní rovnici zapsanou v konzervativním tvaru
se integrují konvektivní, difúzní i zdrojové členy přes kontrolní objem CV
( ) ( )u S
Gaussova věta umožňuje nahradit integrály přes objem integrály přes povrch kontrolního objemu
CV – kontrolní objem (buňka, cell)
A – povrch CV (faces)
NAP10 Řešení transportních rovnic
( ) ( )d d du Sdx dx dx
A W P E B
w e
( ) ee e e
PE
F u Dx
( ) ww w w
WP
F u Dx
( ) ( )A A CV
n u dA n dA S dV
U 1D problémů je kontrolní objem vlastně jen úsečka
( ) ( )e w e Ew P w Pe WF F D D S
Difuzní tok hranicí kontrolního objemu eKonvektivní tok hranicí
kontrolního objemu e
NAP10 Řešení transportních rovnic
2 2w e
P P W W E E W w E eF F
a a a a D a D
1 1( ) ( )2 2e P E w P W
( ) ( )e w e Ew P w Pe WF F D D S
Největší problém spočívá ve stanovení hodnot e w na hranicích kontrolního objemu. To nejsou přímo počítané uzlové hodnoty, protože se řeší soustavy rovnic jen pro neznámé W P E v uzlech (těžištích kontrolních objemů)
Hodnoty na hranici je třeba interpolovat, nejpřirozenější se zdá lineární interpolace
V tomto případě bude rovnice kontrolního objemu (pro nulový zdrojový člen S)
Výsledná aproximace (tzv. centrální schema) má druhý řád přesnosti, ale všechny její koeficienty aP,aW,aE budou kladné jen když je konvektivní tok relativně malý: Jen tehdy splňuje schema podmínku minimaxu zajišťující absenci oscilací.
| | 2e e PEe
e e
F u xPeD
Pe je Pecletovo číslo kontrolního objemu, dá se zmenšit zjemňováním sítě. Při transportu hybnosti se nazývá Reynoldsovo číslo buňky.
NAP10 Řešení transportních rovnicDosáhnout toho, aby Pecletovo (Reynoldsovo) číslo buňky bylo ve všech kontrolních objemech menší než 2, lze při vyšších rychlostech proudění a nízké viskozitě dosáhnout jen za cenu extrémně husté sítě (uvažujte třeba proudění vody s rychlostí 2 m/s – velikost kontrolního objemu nesmí být větší než 10-6 m, mikrometr). Nesplnění této podmínky vede k nekorektním oscilacím řešení, které většinou ani nekonverguje. Východiskem je používání schemat, které podmínku minimaxu (boundedness) splňují i při vyšších nebo dokonce libovolných hodnotách Reynoldsova čísla buňky, někdy ovšem za cenu snížení řádu přesnosti.
Schemata tohoto typu se nazývají UPWIND (protiproudá), interpolují hodnoty na hranici dle směru proudění (interpoluje se z hodnot odkud tekutina přitéká). Nejjednodušší je protiproudé schéma prvního řádu přesnosti
max ,0 max ,0P P W W E E W w w E e ea a a a D F a D F
pro směr proudění vpravo ( 0) jinak pro směr proudění vlevo ( 0) jinak
e P e e E
w P w w W
FF
Toto schema splňuje požadavek minimax pro libovolné Reynoldsovo číslo buňky
NAP10 Řešení transportních rovnicNaprosté stejné principy a postupy platí i u 2D a 3D transportních rovnic. Například ve Fluentu můžete vybírat centrální i protiproudá schemata,
= V A
dV dA
0
1f
fA
f
ff AV
1)(
21
10 f
Upwind prvního řádu
Upwind druhého řádu 0f s 0 f
Metody UPWIND přepisu konvektivních členů se dají realizovat i v metodě konečných prvků tím, že se použijí asymetrické váhové funkce w(x,y,z), jejichž těžiště je posunuté proti směru proudění (varianta Galerkin-Petroff).
NAP10
Beckman
NS rovnice se řeší nejčastěji metodou konečných objemů (FVM), konečných diferencí (FD) nebo konečných prvků (FEM). Používají se odlišné varianty pro
Stlačitelné proudění (konečná rychlost zvuku, diskontinuity rychlostí, rázové vlny)
Proudění téměř nestlačitelné kapaliny (nekonečná rychlost zvuku)
Převažují explicitní metody (hyperbolické rovnice) uvedené v předchozích přednáškách (metoda charakteristik, Lax Wendroff) a další (McCormack).
Převažují implicitní varianty (rovnice jsou parabolické a eliptické)
Transformace NS rovnic do nových proměnných (vířivost a proudová funkce). Je tím eliminován tlak a automaticky je přesně splněna rovnice kontinuity. Výhodné zejména pro 2D proudění.
Formulace v primitivních proměnných (rychlosti a tlak). Tlak představuje omezující podmínku, která má zajistit splnění rovnice kontinuity, jenomže tlak se v rovnici kontinuity neobjevuje. Pro řešení se využívají metody tlakových korekcí SIMPLE, SIMPLER, SIMPLEC a PISO.
Řešení NS rovnic
NAP10 Proudová funkce-vířivostMetoda proudové funkce a vířivosti (2D proudění nestlačitelné kapaliny) aplikovaná na NS rovnice a rovnici kontinuity
vyp
yvv
xvu
uxp
yuv
xuu
2
2
1
1
Z těchto bilancí hybnosti lze eliminovat tlak tím, že první rovnici derivujeme dle y, druhou dle x, odečteme a upravíme užitím rovnice kontinuity
2( ) ( ) ( )u v u v u vu vx y x y y x y x
Výraz v závorkách je zetová složka vektoru vířivosti
,u vy x
wvuzyx
kji
u
2u vx y
Složky rychlostí u,v můžeme vyjádřit pomocí jediné skalární funkce tak, aby byla automaticky splněna rovnice kontinuityVířivost lze vyjádřit proudovou funkcí přímo z definice
2
0
yv
xu
yu
xv
NAP10 Proudová funkce-vířivostVířivost popisuje rotaci kapaliny a může být charakterizována alternativně vektorem nebo tenzorem
)(21
0
0
0
21))((
21
j
i
i
jij
T
xu
xu
zv
yw
zu
xw
zv
yw
yu
xv
zu
xw
yu
xv
uu
( , , ) jk kij
i
uw v u w v uuy z z x x y x
Všimněte si, že antisymetrický tenzor má jen 3 nezávislé složky, které de facto tvoří 3 prvky vektoru vířivosti . Ta ½ v definici tenzoru souvisí s rozkladem tenzoru gradientu rychlosti na symetrický tenzor rychlosti deformace a antisymetrický tenzor vířivosti
))((21))((
21 TT uuuuu
NAP10 Proudová funkce-vířivostUvažujme tuhé těleso, rotující kolem osy z úhlovou rychlostí z
r
x
y
ff
ff
xrv
yru
cos
sin
fyu
xv 2
vířivost je tedy v každém bodě tuhého rotujícího tělesa konstantní
NAP10
Místo trojice parciálních rovnic pro u,v,p máme jen dvojici rovnic pro neznámé , , které má navíc výhodu, že rychlosti odvozené z proudové funkce splňují rovnici kontinuity naprosto přesně. Transport vířivosti je parabolická rovnice a definice vířivosti eliptická (tzv.Poissonova rovnice)
| , |w w w wu vy x
2u vx y
Z předepsaných rychlostí na hranici (např. na stěně) je tedy třeba odvodit i okrajovou podmínku pro vířivost.
Jak se to dělá uvedeme na příkladu toku kanálem
2
Čára konstantní hodnoty je proudnice, např. na stěně proto platí okrajová podmínka =const. Na stěně ovšem musí být splněna i podmínka lpění kapaliny, jsou tedy předepsané obě složky rychlostí. Pro diferenciální rovnici druhého řádu však lze předepsat jen jedinou okrajovou podmínku.
Proudová funkce-vířivost
NAP10
( | , | )w w w wu vy x
x
y
Osa kanálu je proudnice takže třeba =0
stěna je proudnice takže i w je konstanta
H
w dyyu0
)(
H
Vstup kanálu, roste od nuly do w třeba =uy
1 2 N-1 N
1
2
M
Na výstupu předpokládáme např. téměř vyvinutý profil a
N-1=N
yu
xv
Proudová funkce
Vířivost
Na ose je vířivost nulová (to plyne přímo z definice).
Na vstupu je zadaný rychlostní profil u1(y) a nulová příčná složka v1(y)=0
yu
x
12
2
1
Pokud by byl na vstupu rychlostní profil vyvinutý, byl by první člen nula. Když ne, aproximujeme ho rozvojem proudové funkce 2
21
2
12
21
21
12 21
21 x
xx
xx
x
2 1 11 2
2( ) ux y
Stejný postup se aplikuje u stěny kde je předepsaná rychlost UM (v tomto případě nulová)
)(212 yU
y MMMM
Proudová funkce vířivost
NAP10 Příklad – tok v dutině 1/5
x
y
H
1 2 N-1 N
1
2
M
UM
Tok v uzavřené čtvercové dutině, jejíž víko se pohybuje konstantní rychlostí (to je často používaný benchmark). Systém nemá žádný přítok a odtok, takže proudová funkce je na celé hranici nulová (rozdíl hodnot proudové funkce je obecně objemový průtok, v tomto případě nulový).
)(212 yU
y MMMM
Netriviální jsou pouze okrajové podmínky pro vířivost na stěně
Vířivost na víku
Vířivost na stěnách dutiny 1 1 2 1 1 1 22 2 2
2 2 2( ) ( ) ( )N N Nx x y
NAP10 Příklad – tok v dutině 2/5Metoda kontrolních objemů aplikovaná na transport vířivosti (upwind prvního řádu). Uvažujeme čtvercovou síť =x=y
x
y
W P E
UM
N
S
))0,max(())0,max(())0,max(())0,max(()4||||( 22222
P
NP
SP
EP
WPP
Pvvuuvu
Eliptická rovnice pro proudovou funkci
PNSEWP
)(1422
Tyto soustavy algebraických rovnic spolu s okrajovými podmínkami je třeba řešit iteračně. Výpočet lze zefektivnit tzv. metodou střídavých směrů, kdy se postupně řeší soustavy rovnic pro jednotlivé řádky eliminační metodou (s tridiagonální maticí soustavy, stejně jako u předchozích 1D problémů)
2
SNPu
2W E
Pv
NAP10 Metoda střídavých směrů ADI 3/5
x
y
výpočet první řady x-implicitní krok (hodnoty na hranici jsou známé, z předchozí iterace se použijí pouze modré uzly). A tak se postupuje k
druhé, třetí,… řadě až nahoru.
x
y
výpočet první řady y-implicitní. Postupuje se tentokrát zleva doprava.
x-implicitní půlkrok ve kterém se přesně splní okrajové podmínky vlevo i vpravo.
y-implicitní půlkrok, ve kterém se přesně splní okrajové podmínky dole i nahoře.
ADI-Alternating Directions Implicit
NAP10 Příklad – tok v dutině 4/5Matlab% tok v dutině. metoda vířivosti a pokutové funkice. Upwind. Metoda% střídavých směrů. h=0.01; um=0.0001; visc=1e-6;% parametry site, a casovy krok relax=0.5; n=21; d=h/(n-1); d2=d^2; niter=50; psi=zeros(n,n); omega=zeros(n,n); u=zeros(n,n); v=zeros(n,n); a=zeros(n,1); b=ones(n,1); c=zeros(n,1); r=zeros(n,1); for iter=1:niter for i=1:n u(i,n)=um; end for i=2:n-1 for j=2:n-1 u(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*d); v(i,j)=(psi(i-1,j)-psi(i+1,j))/(2*d); end end
%stream function x-implicitfor j=2:n-1 for i=2:n-1 a(i)=-1/d2;b(i)=4/d2;c(i)=-1/d2; r(i)=(psi(i,j-1)+psi(i,j+1))/d2+ omega(i,j); end r(1)=0;r(n)=0; ps=tridag(a,b,c,r,n); for i=2:n-1 psi(i,j)=(1-relax)*psi(i,j)+relax*ps(i); endend%stream function y-implicitfor i=2:n-1 for j=2:n-1 a(j)=-1/d2;b(j)=4/d2;c(j)=-1/d2;r(j)=(psi(i-1,j)+psi(i+1,j))/d2+ omega(i,j); end r(1)=0;r(n)=0; ps=tridag(a,b,c,r,n); for j=2:n-1 psi(i,j)=(1-relax)*psi(i,j)+relax*ps(j); endend% vorticity boundary conditionsfor i=1:n omega(i,n)=relax*omega(i,n)+(1-relax)*2/d2* (psi(i,n)-psi(i,n-1)-um*d); omega(i,1)=relax*omega(i,1)+(1-relax)*2* (psi(i,1)-psi(i,2))/d2; omega(1,i)=relax*omega(1,i)+(1-relax)*2*(psi(1,i)-psi(2,i))/d2; omega(n,i)=relax*omega(n,i)+(1-relax)*2*(psi(n,i)-psi(n-1,i))/d2;end
%vorticity x-implicitfor j=2:n-1 for i=2:n-1 up=u(i,j);vp=v(i,j); a(i)=-visc/d2-max(up,0)/d; b(i)=4*visc/d2+(abs(up)+abs(vp))/d; c(i)=-visc/d2-max(-up,0)/d; r(i)=omega(i,j-1)*(visc/d2+max(vp,0)/d)+ omega(i,j+1)*(visc/d2+max(-vp,0)/d); end r(1)=omega(1,j);r(n)=omega(n,j); ps=tridag(a,b,c,r,n); for i=2:n-1 omega(i,j)=(1-relax)*omega(i,j)+relax*ps(i); endend%vorticity y-implicitfor i=2:n-1 for j=2:n-1 up=u(i,j);vp=v(i,j); a(j)=-visc/d2-max(vp,0)/d; b(i)=4*visc/d2+(abs(up)+abs(vp))/d; c(j)=-visc/d2-max(-vp,0)/d; r(j)=omega(i-,j)*(visc/d2+max(up,0)/d)+ omega(i+1,j)* (visc/d2+max(-up,0)/d); end r(1)=omega(i,1);r(n)=omega(i,n); ps=tridag(a,b,c,r,n); for j=2:n-1 omega(i,j)=(1-relax)*omega(i,j)+relax*ps(j); endend end
NAP10 Příklad – tok v dutině 5/5
5 10 15 20 25 30
5
10
15
20
25
30
-9
-8
-7
-6
-5
-4
-3
-2
-1
0x 10
-8
5 10 15 20 25 30
5
10
15
20
25
30
-1
0
1
2
3
4
5
6
7
8
9x 10
-5
u
5 10 15 20 25 30
5
10
15
20
25
30
-0.4
-0.35
-0.3
-0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
5 10 15 20 25 30
5
10
15
20
25
30
-3
-2
-1
0
1
2
3x 10
-5
v
Re=1
10 20 30 40 50 60 70
10
20
30
40
50
60
70
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0x 10
-6
Re=1000
NAP10
Řešení NS rovnic v primitivních proměnných (u,v,p) je nejčastější, koneckonců ho používá i Fluent.
Musí se ale vypořádat se dvěma problémy:
1. Oscilace tlaků a rychlostí, tzv. šachovnicový vzor (checkerboard pattern) způsobený vynecháním centrálních hodnot tlaků při aproximaci gradientů
2. Jak počítat tlaky z rovnice kontinuity v níž se u nestlačitelné kapaliny tlak vůbec nevyskytuje.
Primitivní proměnné
Bellin
NAP10 Primitivní proměnné
2
2
( ) ( )
( ) ( )
0
u
v
u uv p u u Sx y x x x y y
uv v p v v Sx y y x x y yu v
x y
Příklad: Stacionární 2D proudění (neznámé jsou rychlosti u,v a tlak p
Rovnice pro u
Rovnice pro v
To je rovnice pro p?
NAP10 Primitivní proměnné
10
10
10
10
10
10
10
10 10
1000
0
0
0
0
0
00 0
0
Problém šachovnicového vzoru. Když napíšeme třeba bilanci hybnosti kontrolního objemu ve směru x zjistíme, že alternativní střídání hodnot tlaků v kontrolních objemech nemá vůbec žádný vliv na rovnováhu ve směru x
Problém se objeví vždy, když jsou všechny primitivní
proměnné lokalizovány v centrech kontrolních objemů.
W P E
2
( ) ( ) uv S
x y x x xu u up
yuy
10 10| 02 2
E WP
p ppx x x
NAP10 Primitivní proměnné Úplně stejný problém se objeví i u rovnice kontinuity. I tato rovnice vůbec necítí divoké (a nesprávné) oscilace rychlostí v buňkách. Numerické řešení je potom superpozicí skutečného rozložení rychlostí a tlaků a víceméně libovolného šachovnicového vzoru, který podle zákona schválnosti nakonec vždy převládne. Řešením je použití posunutých kontrolních objemů pro bilance hybností tak, aby se uzlové tlaky ocitly na hranicích kontrolních objemů (staggered grid). Fluent řeší tento problém speciálním typem interpolace rychlostí Rhie-Chow.
10
10
10
10
10
10
10
10 10
1000
0
0
0
0
0
00 0
0
Kontrolní objem pro hybnost ve směru y
Kontrolní objem pro rovnici kontinuity
Kontrolní objem pro hybnost ve směru x
NAP10 Metoda SIMPLE
, , 1 ,
* * * *, , ,( )
I j nb I J I JI j nb I j I jneighbours
a v va A p p b
SIMPLE neznamená jednoduchá, ale Semi Implicit Method for Pressure Linked Equations.
První krok každé iterace spočívá v odhadu tlaku p*
, 1, ,
* * * *, , ( )
i J nb I J I Ji J nb i J iJneighbours
a u Au a p p b
Stanoví se uzlové rychlosti u* ve směru x a rychlosti v* ve směru y z transportních rovnic hybnosti (použije se třeba schema typu UPWIND nebo CENTRAL dle hodnot Re).
Tyto rychlosti splňují bilance hybnosti, ale ne rovnici kontinuity, proto následuje druhý krok, korekce tlaků p’ a korekce rychlostí u’ v’ (další folie)
NAP10 Metoda SIMPLE * * *' ' 'p p p u u u v v v
' ' 'IJ IJ nb nb IJneighbours
a p a p b
, 1, ,
* * * *, , ( )
i J nb I J I Ji J nb i J iJneighbours
a u Au a p p b
, , , 1, ,( )i J i J nb nb i J I J I J iJ
neighbours
a u a u A p p b “správné” řešení
, 1, ,
' ' ' ', , ( )
i J nb I J I Ji J nb i Jneighbours
a au u A p p
odečti Zanedbáme!!!
, 1, , , 1 ,
' ' ' ' ' ', , ,( ) ( )
i J I J I J I J I Ji J I j I jd p p d p pu v
Dosazení u*+u’ a v*+v’ do rovnice kontinuityŘešení rovnic pro
korekce tlaků
NAP10 Metoda SIMPLE
*, , ,
*, , , 1, ,
*, , , , 1 ,
'
( ' ' )
( ' ' )
I J I J I J
i J i J i J I J I J
I j I j I j I J I J
p p p
u u d p p
v v d p p
Tyto zpřesněné tlaky se použijí jako nová výchozí
iterace
Třetí krok: aktualizace tlaků a rychlostí z vypočtených korekcí
Další iterace (návrat k prvnímu kroku SIMPLE). Iterace se opakují tak dlouho dokud nejsou korekce tlaků dostatečně malé.
NAP10 Metoda SIMPLE
Kromě základní varianty SIMPLE existují i další: SIMPLEC, SIMPLER, PISO. Nemyslím si, že je důležité abyste tyto algoritmy znali, důležité je vědět to, že ovlivňují rychlost výpočtu (počet iterací potřebných na korekci tlaku), ale ne výslednou přesnost numerického řešení. Ta totiž závisí jen na zvolené interpolační metodě a tudíž na přesnosti diferenční náhrady (třeba na tom, zda použijete CENTRAL nebo UPWIND).
NAP10 Příklad FLUENT tok v dutině Re=1
u v
5 10 15 20 25 30
5
10
15
20
25
30
-1
0
1
2
3
4
5
6
7
8
9x 10
-5
u
5 10 15 20 25 30
5
10
15
20
25
30
-3
-2
-1
0
1
2
3x 10
-5
v
pro porovnání jsou zde uvedeny dříve získané výsledky metodou proudové funkce a vířivosti (Re=1 ale hrubší síť)
