41
Numerické řešení stlačitelného proudění v kanále a mříži Autor: Milan ŽALOUDEK Vedoucí práce: Doc. Ing. Jaroslav FOŘT, CSc. ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE FAKULTA STROJNÍ
Numerické řešení stlačitelného proudění v kanále a mříži
Autor: Milan ŽALOUDEK
Vedoucí práce: Doc. Ing. Jaroslav FOŘT, CSc.
ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZEFAKULTA STROJNÍ
Výchozí rovnice• zákon zachování hmoty
• zákon zachování hybnosti
• zákon zachování energieEulerovy rovnice
0W F G
t x y
uW
v
e
2
u
u pF
uv
e p u
2
v
uvG
v p
e p v
2 211
2p e u v
uzavírací vztah
normování 0 0, ,p L
Matematická formulace úlohysystém nelineárních hyperbolických rovnic slabé řešení
hledáme funkci W(x,y,t) na oblasti 2+ :
• W K() – třída funkcí, ve které připouštíme existenci spočetně mnoha křivek, podél nichž funkce W nabývá různých konečných
limitních hodnot zleva a zprava
tyto křivky nazýváme rázové vlny popř. nespojitosti I. druhu
mimo tyto křivky je funkce W spojitá
• W splňuje rovnici
pro libovolné t2>t1 a libovolnou oblast D s dostatečně hladkou hranicí
• W splňuje počáteční podmínky W(t=0)=W0
• W splňuje okrajové podmínky
2
1
0t
t x y
t D
W F G dxdydt
Okrajové podmínky
• objevují se celkem 4 základní druhy
(vstup, výstup, stěna, periodicita)
• při formulaci vycházíme z jednodimenzionální analýzy
• podzvuková rychlost v normálovém směru ke vstupní hranici
k výstupní hranici
Obecná výpočtová oblast rovinného kanáluObecná výpočtová oblast lopatkové mříže
Numerické řešení úlohy• cell-centered (hodnoty proměnných v těžišti objemu)
diskretizace základních rovnic:
R: obdélníkové pravidlo numerického integrování
• výpočtová oblast je pokryta strukturovanou čtyřúhelníkovou sítí
2
1
Greenova v.
v.o stř.hodnotě 10
i
t
t x y tit D D
W F G dxdydt W Fdy Gdx
1n ni i
t
W WW
t
4
1, ,
1,
n ni j i j k k k k
ki j
tW W F y G x
L: Eulerova dopředná aproximace
Aproximace toku v 1D pomocí numerické metody AUSM
schéma je založeno na struktuře řešení Riemannova problému
Eulerovy rovnice v 1D:
vlastní čísla , kde
Machovo číslo rozhoduje o druhu režimu a o počtu kladných a záporných vlastních čísel
tok F rozdělíme na advektivní a tlakovou část
dále upravíme
tento tok aproximujeme pomocí hodnot z L a R jako jejich vhodnou kombinaci
proto použité hodnoty formálně přeznačíme
AUSM = zkratka ang. Advection Upstream Splitting Method
0W F
t x
W u
e
2
u
F u p
e p u
1 2 3, ,u u a u a pa
uM
a
2
0
0
u
F u p
e p u
0 0
0 0
a
F M ua p M p
e p a
, ,LR LR LRM M p p
0
0LR LR LRF M p
Aproximace toku F na hranici mezi i-tou a (i+1)-ní buňkou k určení MLR , pLR používáme tzv. rozkládající (splitting) polynomy
Označení rozkládajících polynomů buňka i buňka i+1
Machova čísla
tlaku
Požadavky na rozkládající polynomy (platí stejně pro M+/- i p+/-):
• M+ , M- spolu se svými prvními derivacemi spojitě závislé na Machově čísle
• vyjádřeny polynomem nejnižšího možného stupně
v každé buňce vyčíslíme jednotlivé polynomy a jejich následnou kombinací získáme MLR a pLR
,LR L R LR L RM M M p p p
,M M
W W
M
W
M
W
výraz LR je dán proměnnými pouze jedné hraniční buňky a to podle znaménka výsledného MLR
0
0
LR LR L L
LR LR R R
M W
M W
• vlastní čísla jsou kladná, vlastní čísla jsou záporná
LM RM
RpLp
Rozkládající polynomy:
buňka i buňka i+1
podzvukový režim
|M|<1
nadzvukový režim
|M|>1
211
4L LM M 211
4R RM M
1
2L L LM M M 1
2R R RM M M
21 2
4L
L L L
pp M M
2L LL
LL
M Mpp
M
21 2
4R
R R R
pp M M
2R RR
RR
M Mpp
M
Rozšíření numerického schématu na 2D• zavedeme kladnou orientaci hran (a, b)
• vektor jednotkové vnější normály
• aproximovaný numerický tok přepíšeme
• zavedeme matici rotace
• Eulerovy rovnice jsou invariantní vůči rotaci dvourozměrný případ je založen na jednodimenzionálním rozkladu v normálovém směru
• o režimu proudění rozhoduje Machovo číslo v normálovém směru
cos ,sin ,y x
ns s
,F y G x F G n s
1 0 0 0
0 cos sin 0
0 sin cos 0
0 0 0 1
T
1, ( ) ... , , ,T
n n nF G n T F TW TW u v e
,u v nM
a
2
2
0 0 0 0
0 0,
0 0
( ) 0 ( ) 0 0 0
L
u v
u p uv u py y x x y xF G n M
uv v p v ps s s s s s
e p u e p v e p
0
0
xR LR LR
y
nM p
n
Zvyšování řádu přesnosti v prostorových proměnných
• základní numerické schéma AUSM je prvního řádu přesnosti v prostoru
• původní schéma používalo pouze hodnoty v těžištích
11/ 2 1/ 2 1/ 2 1/ 2; ;n n L R L R
i i i i i i
tW W F W W F W W
x
1/ 2 1/ 2R L
i iW W • zvyšování řádu přesnosti je založeno na náhradě těchto hodnot „přesnějšími“
lineární rekonstrukce
+limiter
Lineární rekonstrukce• strukturovaná čtyřúhelníková síť rekonstrukce ve 2 nezávislých směrech
• na každé buňce 2 lokálně jednodimenzionální rekonstrukce
• několik způsobů, jak spočítat rekonstrukci v buňce i
1) upwind
2) downwind
3) centrálně
1
1
i iU
W W
x
1
2
i iD
W W
x
1 1
1 2
i iC
W W
x x
• nové hodnoty1/ 2
1/ 2
2
2
R ii i X
L ii i X
xW W
xW W
Výpočet rekonstruovaných hodnot v buňce i
Limitery (omezovače)od prováděných úprav požadujeme neoscilativní chování
rekonstruované hodnoty na hranici i-1/2 musí ležet mezi Wi-1 a Wi
rekonstruované hodnoty na hranici i+1/2 musí ležet mezi Wi a Wi+1
samotná rekonstrukce 1), 2) nebo 3) nezaručí neoscilativní chování
doplnění o vhodný limiter
minmod limiter
MC limiter
superbee limiter
Barthův limiter
minmod minmod ;D U
minmod ;minmod 2 ;2MC C D U
definujme funkce minmod a maxmod
0
minmod ; 0
0 0
a a b a b
a b b b a a b
a b
maxmod ;a a b
a bb b a
superbee maxmod minmod ;2 ;minmod 2 ;D U D U
Barth ...
speciální úprava rekonstrukce v okrajových buňkách
Numerická aproximace okrajových podmínekPředpokládáme: pracovní médium uloženo ve velkém zásobníku ze zásobníku
dopraveno izoentropicky na vstupní hranici proudění výpočtovou oblastí výstup do prostředí se známým tlakem
VSTUP + VÝSTUP
- podle druhu zadáme vždy vhodný
počet parametrů
- zbývající veličiny extrapolujeme
ze 2 sousedních buněk
STĚNA - podmínka neprostupnosti stěny
PERIODICITA - periodická hranice má regulární buňky na obou stranách
- numerický tok počítáme podle vzorce pro regulární hranici
, 0u v n
Výsledkyvlastnosti vyvinutého programu byly testovány na několika případech
• koleno
(kanál konstantního průřezu
s otočením proudu o 90)
• GAMM kanál
• lopatková mříž DCA 8%
Koleno• výpočtová síť 16035 buněk
• izočáry Machova čísla (zvýrazněná izočára M=1, přírůstek M=0.02)
1. řád přesnosti 2
0
0,620p
p
2
0
0,591p
p
vyšší řád přesnosti
(minmod limiter)1. řád přesnosti
numerické schéma Ron-Ho-Ni
výsledek převzatý z [1]
[1] Halama J.: 2D stacionární nevazké proudění v kanále, sem. práce z Vnitřní aerodynamiky, ČVUT, 1996
GAMM kanál• hrubá síť 90 30 buněk, jemná síť 150 45 buněk
• tlakový poměr odpovídá výstupnímu Machovu číslu M2=0,675
1. řád přesnosti – hrubá síť vyšší řád přesnosti
(MC limiter – hrubá síť)
vyšší řád přesnosti
(MC limiter – jemná síť)
GAMM kanál – průběh Machova čísla podél stěn
porovnání s výsledky převzatými z [2]
jiná síť, jiné numerické schéma
(TVD MacCormack, Implicit WENO)
vlastní výsledky:
• jemná síť
• minmod limiter
• Barthův limiter
[2] Kozel K., Fűrst J.: Numerické metody řešení problémů proudění I
ČVUT Praha, 2001
lopatková mříž DCA 8%
Parametry výpočtu:
• vyšší řád přesnosti s minmod limiterem
• výpočtová síť 120 40 buněk
• tlakový poměr odpovídá výstupnímu Machovu číslu M2=0,833
• úhel nabíhajícího proudu =0,9
ZávěrCílem práce bylo vyvinout a odladit vlastní numerický program,
pro řešení nevazkého stlačitelného proudění, založený na
numerickém schématu AUSM.
Tento cíl byl splněn.
Dosahované výsledky jsou ve shodě s jinými numerickými
výsledky i s fyzikálními předpoklady proudění.
Další vývoj programu:
• implementace a testování dalších variant AUSM schématu
• přechod na stlačitelné vazké proudění
• rozšíření na 3D úlohy
Děkuji za pozornost
Děkuji za pozornost
Výchozí rovnice
• zákon zachování hmoty
• zákon zachování hybnosti
• zákon zachování energie
0wt
::::::::::::::
Dwp f
Dt
:::::::::::::: ::::::::::::: :
:de
q p w Qdt
::::::::::::: :::::::::::::::::::::::::::::::::::::::: :::
hustota
rychlost
energie
tlak
tečné napětí
objemová síla
hustota objemového toku
rychlost deformace
čas
::::::::::::::::::::::::::::
w::::::::::::::
p
qf::::::::::::::
e
t
zjednodušující předpoklady
• rovinné proudění
• nevazká tekutina
• nulové hmotové síly
• žádné zdroje tepla
Eulerovy rovnice:
,w u v
0
0f ::::::::::::::
0Q
0W F G
t x y
uW
v
e
2
u
u pF
uv
e p u
2
v
uvG
v p
e p v
2 211
2p e u v
uzavírací vztah
normování 0 0, ,p L
Matematická formulace úlohysystém nelineárních hyperbolických rovnic
slabé řešení
hledáme funkci W(x,y,t) na oblasti 2+ :
• W K()
• W splňuje rovnici
pro libovolné t2>t1 a libovolnou oblast D s dostatečně hladkou hranicí
• W splňuje počáteční podmínky W(t=0)=W0
• W splňuje okrajové podmínky
2
1
0t
t x y
t D
W F G dxdydt
Okrajové podmínky
• celkem se objevují 4 základní druhy
(vstup, výstup, stěna, periodicita)
• vycházíme z jednodimenzionální analýzy
• podzvukový vstup i výstup
Obecná výpočtová oblast rovinného kanáluObecná výpočtová oblast lopatkové mříže
Numerické řešení úlohy• metoda konečných objemů (FVM)
• cell-centered (hodnoty proměnných v těžišti objemu)
R: obdélníkové pravidlo numerického integrování
L: Eulerova dopředná aproximace
• výpočtová oblast je pokryta strukturovanou čtyřúhelníkovou sítí
2
1
10
i
t
t x y tit D D
W F G dxdydt W Fdy Gdx
1n ni i
t
W WW
t
4
1, ,
1,
n ni j i j k k k k
ki j
tW W F y G x
suma na R straně je aproximována pomocí numerického schématu
Numerické schéma AUSM v 1D
schéma je založeno na struktuře řešení Riemannova problému
Eulerovy rovnice v 1D:
vlastní čísla , kde
Machovo číslo rozhoduje o druhu režimu a o počtu kladných a záporných vlastních čísel
tok F rozdělíme na advektivní a tlakovou část
dále upravíme
tento tok aproximujeme pomocí hodnot z L a R jako jejich vhodnou kombinaci
proto použité hodnoty formálně přeznačíme
AUSM = zkratka ang. Advection Upstream Splitting Method
0W F
t x
W u
e
2
u
F u p
e p u
1 2 3, ,u u a u a pa
uM
a
2
0
0
u
F u p
e p u
0 0
0 0
a
F M ua p M p
e p a
, ,LR LR LRM M p p
0
0LR LR LRF M p
Aproximace toku F na hranici mezi i-tou a (i+1)-ní buňkou k určení MLR , pLR používáme tzv. rozkládající (splitting) polynomy
Označení rozkládajících polynomů buňka i buňka i+1
Machova čísla M+ M-
tlaku p+ p-
Požadavky na rozkládající polynomy (platí stejně pro M+/- i p+/-):
• M+ , M- spolu se svými prvními derivacemi spojitě závislé na Machově čísle
• vyjádřeny polynomem nejnižšího možného stupně
• vlastní čísla jsou kladná, vlastní čísla jsou záporná
v každé buňce vyčíslíme jednotlivé polynomy a jejich následnou kombinací získáme MLR a pLR
,LR LRM M M p p p
,M M
W W
M
W
M
W
výraz LR, je dán proměnnými pouze jedné hraniční buňky a to podle znaménka výsledného MLR
0
0
LR LR L L
LR LR R R
M W
M W
Rozkládající polynomy:
buňka i buňka i+1
podzvukový režim
|M|1
nadzvukový režim
|M|>1
211
4 LM M 211
4 RM M
1
2 L LM M M 1
2 R RM M M
21 2
4L
L L
pp M M
2L LL
L
M Mpp
M
21 2
4R
R R
pp M M
2R RR
R
M Mpp
M
Rozšíření numerického schématu na 2D• zavedeme kladnou orientaci hran (a, b)
• vektor jednotkové vnější normály
• aproximovaný numerický tok přepíšeme
• zavedeme matici rotace
• Eulerovy rovnice jsou invariantní vůči rotaci dvourozměrný případ je založen na jednodimenzionálním rozkladu v normálovém směru
• o režimu proudění rozhoduje Machovo číslo v normálovém směru
cos ,sin ,y x
ns s
,F y G x F G n s
1 0 0 0
0 cos sin 0
0 sin cos 0
0 0 0 1
T
1, ( ) ... , , ,T
n n nF G n T F TW TW u v e
,u v nM
a
2
2
0 0 0 0
0 0,
0 0
( ) 0 ( ) 0 0 0
L
u v
u p uv u py y x x y xF G n M
uv v p v ps s s s s s
e p u e p v e p
0
0
xR LR LR
y
nM p
n
Numerická aproximace okrajových podmínek VSTUP - médium uloženo ve velkém zásobníku, kde má klidové parametry p0, 0
- ze zásobníku dopraveno izoentropicky na vstupní hranici
- zadáváme klidové parametry p0, 0, úhel náběhu - z proudového pole extrapolujeme Machovo číslo Min
1 1
1 22
0 0
1cos 1
2in in
uM M
p
1
12
0
11
2 inM
2 212
0
1 11
1 2 2in
e u vM
p
1 1
1 22
0 0
1sin 1
2in in
vM M
p
, , ,in in in inp a M
0
0
in in inin
yF y G x M s p
x
Numerická aproximace okrajových podmínek VÝSTUP - médium vystupuje z výpočtové oblasti do prostředí se známým tlakem p2
- tento tlak je dán poměrem
- z proudového pole extrapolujeme první 3 složky vektoru W
- zadáváme tlakový poměr
- čtvrtou složku W dopočítáme podle
STĚNA - idealizovaný model nevazké stěny (žádná mezní vrstva, rychlostní profil...ap.)
- podmínka neprostupnosti stěny
- tlak na stěně pwall nahrazujeme tlakem v nejbližší buňce přilehlé ke stěně
PERIODICITA - periodická hranice má regulární buňky na obou stranách
- numerický tok počítáme podle vzorce pro regulární hranici
2 2
2
0 0
1
1 2out out
out
u vpe
p p
2
0
p
p
, 0u v n
0
0
xwallwall
y
nF y G x p
n
po úpravách
Zvyšování řádu přesnosti v prostorových proměnných
• základní numerické schéma AUSM je prvního řádu přesnosti v čase a prostoru
• původní schéma používalo pouze hodnoty v těžištích
11/ 2 1/ 2 1/ 2 1/ 2; ;n n L R L R
i i i i i i
tW W F W W F W W
x
1/ 2 1/ 2R L
i iW W • zvyšování řádu přesnosti je založeno na náhradě těchto hodnot „přesnějšími“
lineární rekonstrukce
+limiter
Lineární rekonstrukce
• strukturovaná čtyřúhelníková síť rekonstrukce ve 2 nezávislých směrech
• 2 lokálně jednodimenzionální úlohy
• několik způsobů, jak spočítat rekonstrukci v buňce i
1) upwind
2) downwind
3) centrálně
1
1
i iU
W W
x
1
2
i iD
W W
x
1 1
1 2
i iC
W W
x x
• nové hodnoty - levá hranice buňky i :
- pravá hranice buňky i :
1/ 2
1/ 2
2
2
R ii i X
L ii i X
xW W
xW W
Limiterypožadujeme neoscilativní chování
rekonstruované hodnoty na hranici i-1/2 musí ležet mezi Wi-1 a Wi
rekonstruované hodnoty na hranici i+1/2 musí ležet mezi Wi a Wi+1
samotná rekonstrukce 1), 2) nebo 3) nezaručí neoscilativní chování
doplnění o vhodný limiter
minmod limiter
MC limiter
superbee limiter
Barthův limiter
minmod minmod ;D U
minmod ;minmod 2 ;2MC C D U
definujme funkce minmod a maxmod
0
minmod ; 0
0 0
a a b a b
a b b b a a b
a b
maxmod ;a a b
a bb b a
superbee maxmod minmod ;2 ;minmod 2 ;D U D U
Barth ...
Rekonstrukce v okrajových buňkách
• předpokládáme, že hodnota daná okrajovou podmínkou je přesná není třeba ji upravovat hodnota na okrajových hranách je bez rekonstrukce
• na následující hraně používáme jednostrannou rekonstrukci bez limiteru
Rekonstrukce s použitím minmod limiteru Porovnání rekonstrukcí s různými limitery
Výsledkyvlastnosti vyvinutého numerického programu byly testovány na několika
případech
• koleno
• GAMM kanál
• lopatková mříž DCA 8%
Koleno• výpočtová síť 16035 buněk
• izočáry Machova čísla (zvýrazněna izočára M=1, přírůstek M=0.02)
1. řád přesnosti2
0
0,620p
p
2
0
0,591p
p 2
0
0,586p
pvyšší řád přesnosti
(minmod limiter)
2
0
0,591p
p
GAMM kanál• hrubá síť 90 30 buněk, jemná síť 150 45 buněk
• tlakový poměr odpovídá výstupnímu Machovu číslu M2=0,675
1. řád přesnosti vyšší řád přesnosti
hrubá síť (MC limiter – hrubá síť)
vyšší řád přesnosti
(MC limiter – jemná síť)
GAMM kanál – průběh Machova čísla podél stěn
porovnání s předchozími výsledky (TVD, WENO)
jiná síť, jiné numerické schéma
vlastní výsledky:
• jemná síť
• minmod limiter
• Barthův limiter
lopatková mříž DCA 8%
Parametry výpočtu:
• vyšší řád přesnosti s minmod limiterem
• výpočtová síť 120 40 buněk
• tlakový poměr odpovídá výstupnímu Machovu číslu M2=0,833
• úhel nabíhajícího proudu =0,9
