Tvarová optimalizace rozváděcí skříně topení osobního automobilu Ing. Tomáš Mužík 1. Úvod Úkolem této práce je navrhnout opatření pro zrovnoměrnění hmotnostního toku proudu vzduchu na výstupech z traktu topení pomocí numerických metod. To znamená numericky vyřešit proudové pole v traktu topení automobilu a rozložení hmotnostních průtoků na jeho výstupech a dále se zabývat zrovnoměrněním těchto průtoků. V první části práce je proveden výpočet proudového pole pomocí komerčního software Fluent 6.3. Ve druhé části práce je řešena rovnoměrnost hmotnostních průtoků na výstupech pomocí tvarové optimalizace traktu topení propojením komerčních software Fluent 6.3, Sculptor a ISight. Řeší se jen proudění v rozvodné skříni, ale v praxi se řeší rozvody topení v celém voze. Tam je pak požadavek na to, aby z rozvodné skříně vycházel hmotnostní průtok proudu vzduchu ze všech kanálů rovnoměrně. Zadání práce bylo formulováno ve spolupráci s firmou Škoda Auto a.s., která dodala rozměry geometrie (obr. 1.1).
210
Výstup 3
Výstup 4
Výstup 2
60
210 Výstup 1
97
85
73,5
145 400
145
Vstup
210 210 Obr. 1.1 Použitá geometrie 2. Výpočet proudového pole pomocí komerčního software Fluent 6.3 2.1 Základní přístup k řešení V této části je obsažen stručný postup, vedoucí od geometrie až k numerickým výstupům, získaným výpočtem v software Fluent 6.3. A dále pak postprocessing těchto dat v software EnSight. Postup prací se dá shrnout do následujících několika bodů: 1
1. Tvorba geometrie a povrchové sítě. Geometrie i povrchová síť byly vytvořeny v software Gambit 2.2.3. 2. Načtení dodané povrchové sítě do programu Tgrid 3.6.8, který je preprocesorem software Fluent 6.3. Vytvoření objemové sítě v programu Tgrid 3.6.8. 3. Provedení výpočtu v programu Fluent 6.3. 4. Získání numerických výstupů ze software Fluent 6.3 a export dat pro vyhodnocení do software Ensight 8.0. 5. Vyhodnocení proudění v software EnSight 8.0 2.2 Metoda výpočtu Výpočty byly prováděny v software Fluent 6.3. V této části bude popsán systém rovnic, pomocí nějž popisuje tento software proudění. Proudění je popsáno soustavou Navierových - Stokesových rovnic, což představuje rovnice kontinuity a rovnice bilance hybnosti, která je doplněna rovnicemi transportními pro kinetickou energii turbulence a pro rychlost disipace. Prouděni je uvažováno za nestlačitelné a vazké a celý model je stacionární. Toto zjednodušení je možné pro hodnoty Machova čísla M do 0,17. Rovnice kontinuity pak bude mít tvar: ∂u i = 0, (1) ∂xi kde ui je rychlost v i-tém směru. Rovnici bilance hybnosti lze zapsat ve tvaru: ∂ (u i u j ) ∂p ∂τ ij ρ =− + + Fi , (2) ∂xi ∂x j ∂x j kde p je statický tlak, Fi jsou vnější objemové síly a τij je tenzor vazkých napětí, který se dá zapsat ve tvaru: ∂u ∂u j 2 ∂u i − µ (3) τ ij = µ i + δ , 3 ∂x ij ∂ ∂ x x i i j
kde µ je dynamická viskozita a δij je Kroneckerův tenzor. Transportní rovnice pro turbulentní proudění kinetickou energii k a rychlost disipace ε, model turbulence rke: u ∂ (ρk ) + ∂ (ρu j k ) = ∂ µ + t ∂k + G k + Gb − ρε − YM + S k , (4) σ k ∂x j ∂t ∂x i ∂x j ∂ (ρε ) + ∂ (ρu j ε ) = ∂ ∂t ∂x j ∂x j
µt µ + σε
∂ε ε2 ε + C S − C + C1ε C 3ε Gb + S ε ρ ε ρ 1 2 k k + νε ∂x j (5)
kde η C1 = max 0.43, η + 5 k η=S
ε
S = 2 S ij S ij
2
µ t reprezentuje turbulentní viskozitu, Gk reprezentuje generování kinetické energie pomocí rychlostních gradientů, Gb reprezentuje generování kinetické energie pomocí vztlaku, YM reprezentuje příspěvek fluktuační dilatace ve stlačitelném proudění celkovou disipační rychlostí, C1ε , C 2ε , C 3ε , S k , S ε jsou empirické konstanty, σ k , σ ε jsou turbulentní Prandtlova čísla pro k a ε. Hodnoty konstant pro turbulentní model rke: C1ε = 1,44 C 2ε = 1,9 C µ = 0,09
σ k = 1,0 σ ε = 1,2 2.3 Výpočet v blízkosti stěn
Modelování proudění u stěny ovlivňuje přesnost výpočtu v celé oblasti. V blízkosti stěny se řešené veličiny velice rychle mění, především 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. Proudění v blízkosti stěny je zde modelováno stěnovou funkcí (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í. V této práci byla použita nerovnovážná stěnová funkce. Kdyby nebyla použita, bylo by třeba podrobně modelovat proudění u stěny včetně vazké podvrstvy v souvislosti s jemností sítě. V případě velkých Reynoldsových čísel uplatnění stěnových funkcí podstatně snižuje nároky na výpočet a poskytuje tak časově nenáročné a přesto dostatečně přesné řešení. Stěnové funkce představují soubor poloempirických vztahů a funkcí, pomocí nichž lze přemostit vzdálenost mezi stěnou a buňkou v blízkosti stěny. Zahrnují zákon stěny pro střední rychlost a teplotu a vztahy pro turbulentní veličiny v blízkosti stěny. Běžně se užívají tři typy stěnových funkcí a to: standardní stěnové funkce, nerovnovážné stěnové funkce a dvouvrstvý model (ten je vhodný pro proudění s malými Reynoldsovými čísly). Pro výpočet proudění v blízkosti stěn byla použita nerovnovážná stěnová funkce (Non - Equilibrium Wall Function), která vychází z teorie Laundera a Spaldinga [1]. Je vhodná pro proudění u stěny vystavené účinkům velkého tlakového gradientu, kde se nedá očekávat splnění podmínky lokální rovnováhy. Rovnice vyjadřující nerovnovážnou stěnovou funkci [10]: 1 1 1 1 ~ UC µ4 k 2 1 ρ C µ4 k 2 y = ln E (6) , τω k µ ρ kde
3
y y − y v y v2 1 dp y v ~ ln + + U =U − , 2 dx ρκ k y v ρκ k µ kde y v je dimenzionální tloušťka vazké podvrstvy
yv ≡
µy v* 1 4
ρC µ k
1 2 P
,
(7)
(8)
kde y v* = 11,225 je konstanta. Nerovnovážná stěnová funkce rozděluje úlohu na dvě části. Na vazkou podvrstvu a plně turbulentní vrstvu. To je popsáno následujícím schématem: 2νk y 2 y 2 , y < yv k , y < y v 0, y < y v , k = y v P ,ε = 3 , (9) τt = τ ω , y > y v k2 , y > yv k P , y > y v Cl y
kde −
3 4
C l = κC µ .
(10) 2.4 Nastavení řešiče
Protože je dosahováno ve všech částech výpočtové oblasti nízkých rychlostí, je proudění považováno za nestlačitelné. Diskretizační schéma pro rovnice popisující hybnost, kinetickou energii a rychlost disipace se řeší numerickým schématem upwindu 1. řádu přesnosti, pro výpočet několika prvních hodnot. Toto schéma je přepnuto na schéma 2. řádu přesnosti pro získání přesnějšího řešení. Není možné z důvodů horší konvergence nasadit na úlohu již od začátku schéma 2. řádu. Schématem prvního řádu se nechá spočítat jen několik desítek prvních iterací na začátku výpočtu, aby se nám vyvinulo proudové pole a ve všech buňkách sítě jsme měli hodnoty. Pak po nastavení schématu 2. řádu je toto schéma schopno úlohu dopočítat. Rovnice soustavy jsou řešeny odděleně (segregated solution method). To znamená že hodnoty veličin jsou nastaveny z předchozího kroku (iterace), v případě počátečního stavu výpočtu jsou nastaveny inicializací. Nejprve jsou řešeny rovnice hybnosti ve všech směrech, čímž se získají nové hodnoty pro složky rychlosti. Poté se řeší rovnice pro turbulenci. Vazba mezi rychlostí a tlakem je zaručena algoritmem SIMPLE. 2.5 Popis řešení
V programu Fluent 6.3 bylo řešeno rozložení hmotnostních průtoků na výstupech z traktu topení při vstupní rychlosti 5,3 m/s. Hodnoty materiálových konstant (hustota, dynamická viskozita, tlak) byly uvažovány takové, jaké jsou standardně nastaveny ve Fluentu a odpovídají reálnému prostředí. Pro vzduch jsou tyto hodnoty: kg ρ = 1,225 3 m µ = 1,789 ⋅ 10 −5 Pa ⋅ s p = 101324 Pa . 4
2.6 Výsledky řešení
Nejprve se zaměřme na rozdělení hmotnostních toků na jednotlivých výstupech. Při porovnání s experimentem provedeným firmou Škoda Auto, a.s., uvedeným v následující tabulce (Tab. 2.1). je vidět, že došlo k velké shodě výsledků. To bezesporu dokládá velký význam numerických metod v řešení takovýchto problémů. Tab. 2.1 Hmotnostní průtok [kg/hod] vstup Měření Výpočet
288,1
výstup 1
výstup 2
výstup 3
výstup 4
63,5
93,5
71,9
59,2
288,0723 -63,9369 -93,9112 -72,1107 -58,1137
Je vidět že hmotnostní toky na výstupech jsou mezi sebou velice odlišné a proto je nutné tento problém řešit. Tím se zabývá druhá část této práce. Ještě zde pro doplnění ukážeme několik hlavních řezů proudovým polem pro lepší představu o problému.
Odtržení
Kontury rychlostí, řez výstupem 2
Kontury rychlostí, řez výstupem 3
Proudnice v traktu topení
Pole vektorů rychlostí v traktu topení
Obr. 2.1 Proudové pole traktu topení 3. Tvarová optimalizace traktu topení
Pro dosažení dobrých výsledků optimalizace bylo třeba vhodně zvolit optimalizační metodu. Při špatné volbě optimalizační metody by došlo k získání velice nepřesných výsledků, nebo k neúměrně velké výpočtové době.
5
3.1 Volba optimalizační metody
V našem případě je nutné zohlednit především výpočetní dobu optimalizace. Čas potřebný ke korektnímu výpočtu cílové funkce na daném tvaru je velký a z tohoto důvodu je téměř vyloučené použít postupné vyhledávání globálních extrémů, tj. především stochastických metod a generických algoritmů, jejichž výpočetní náročnost je zatím neúnosná. Vzhledem k tomu, že známe analytickou podobu cílové funkce (jak byla vytvořena, bude popsáno později), je pro nás nejvhodnější gradientová metoda. Ta je vhodná pro úlohy, kde lze předpokládat výskyt jednoho extrému, jelikož je schopna nalézt pouze nejbližší lokální extrém. Byla zvolena a použita gradientová metoda LSGRG2 (Generalized Reduced Gradient Methods). Tato metoda je vhodná k řešení nelineárních, nespojitých problémů což přesně vystihuje náš případ. Metoda využívá k nalezení extrému směr největšího spádu. 3.2 Metoda největšího spádu
Metoda největšího spádu pro nalezení lokálního minima funkce f (x) postupuje následovně. Začíná v nějakém inicializačním bodě x 0 a počítá ∇f ( x o ) (gradient funkce). Udělá krok ve směru největšího spádu, − ∇f ( x 0 ) , užitím kroku o délce α 0 , k získání nového bodu x1 . Tento krok se opakuje, dokud se procedura nezastaví na nějakém uspokojivém kroku. Tento proces se dá popsat následující relací: dáno x 0 xi +1 = xi − α i ∇f ( xi ) i = 0,1,2,... (11) kde α i > 0 . Proces bude konvergovat přinejmenším k lokálnímu minimu funkce f (x) jestliže α i je voleno tak, že f ( x i +1 ) < f ( x i ) (12) pro ∀ i, jestliže se funkce snižuje v každém kroku. Protože je funkce od začátku sestupná, má daný směr − ∇f ( x i ) , kde vždy existuje α i > 0 tak aby vyhovovala (12). 3.3 Délka kroku a optimální gradient
Jedna z cest, jak nalézt α i vyhovující (12) je zvolit α i minimalizováním funkce g (α ) = f [x i − α∇f ( x i )] . (13) Všimněme si že x i a ∇f ( x i ) jsou známé vektory tak, že jedinou proměnou je zde α . Adaptaci na metodu největšího sestupu, která užívá (13), nazýváme metodou optimálního gradientu, popsanou takto dáno x 0 s i = −∇f ( x i ) Zvolme α = α i minimalizováním v (13), x i +1 = x i + α i s i (14) Nastavme i = i + 1 a opakujme. Geometricky je α i voleno minimalizováním funkce f (x) podle směru s i , začínající v x i . Na minimu, dg ′ = s i ∇f ( x i + α i s i ) = 0 , (15) dα α =α i
6
kde vektor xi + α i si musí být tangentou k vrstevnici v α = α i , dg je pak pro malé změny dα rovno nule. Jestliže ∇f ( x i +1 ) je normálou k té samé vrstevnici, následné kroky jsou na sebe vzájemně kolmé. 3.4 Stop kritéria
Některá možná stop kritéria jsou ukázána v následujícím: ∂f 1. Až minimum = 0 , zastaví se když ∂x i
(a)
∂f <ε ∂x i
i=1,2,3,… 2
∂f < ε (b) ∑ i −1 ∂x i 2. Zastaví se, když změny ve funkci jsou menší než nějaká limitní hodnota η, například, f ( x i +1 ) − f ( x i ) < η n
3.5 Cílová funkce
Jako cílová funkce byla pro optimalizaci zvolena rovnoměrnost hmotnostních průtoků na výstupech. Vzhledem k tomu, že potřebujeme funkci minimalizovat, jak bylo popsáno dříve, byla vytvořena jako rozdíl hmotnostních průtoků na vstupu a výstupech. Tato funkce byla zvolena takto: 4 m & vstup rovnoměrnost = ∑ (16) − m& i , 4 i =1 kde m& vstup vyjadřuje hmotnostní průtok na vstupu a m& i vyjadřuje hmotnostní průtoky na jednotlivých výstupech. 3.6 Postup vlastní optimalizace
Tvarová optimalizace traktu topeni automobilu byla provedena v komerčním software ISight, s využitím gradientové metody LSGRG2, o níž bylo pojednáno dříve. Proudové pole v jednotlivých krocích optimalizačního procesu bylo řešeno programem Fluent. Software byl nastaven stejně jako v předchozí úloze. Pro deformaci sítě byl použit software Sculptor. Zde bylo třeba nadefinovat také body, kterými se pohybovalo. Nebylo totiž možné, deformovat celou geometrii traktu topení vzhledem k jeho zabudování do automobilu. Tvar krabice musel zůstat nezměněn a poloha a tvar výstupů také. Po této omezující podmínce nám zbylo jen měnit tvar kanálů. Ve Sculptoru bylo třeba nadefinovat tak zvaný ASD objem, což je objem, ve kterém je možno pohybovat geometrií. Jak byl tento objem zvolen, je zřejmé z Obr. 3.1. Jelikož nebylo možné měnit tvar celé geometrie, bylo se třeba v tomto ASD omezit na několik bodů (9 kolem každého z kanálů), kterými se dále pohybovalo. Tyto body kolem jednoho z kanálů jsou znázorněny také na Obr. 3.1 (směr pohybu je znázorněn šipkami). Dále bylo třeba omezit rozsah pohybu těchto bodů, vzhledem k deformaci sítě. Kdyby se síť zdeformovala příliš mnoho, došlo by k vzniku příliš šikmých buněk a na takovéto síti bychom nebyli schopni získat korektní výpočet.
7
Obr. 3.1 3.6.1 Průběh optimalizačního procesu
Optimalizační proces je realizován propojením tří programů a to Fluent, Skulptor a ISight. Fluent je zde užíván pro výpočet proudového pole a hmotnostních průtoků pomocí CFD metod. Skulptor obstarává morfing výpočtové sítě. Isight provádí vlastní optimalizaci. Pro lepší představu je toto propojení znázorněno na Obr. 3.2. Tyto tři software si mezi sebou vyměňují informace pomocí textových souborů ve formátu *.jou. Tento proces je prováděn dokud není nalezeno minimum, nebo nenarazí na hraniční kriteria zadaná uživatelem. V nasem případě je to omezení z důvodu přílišné deformace sítě, jak již bylo řečeno.
CFD Processing •Výpočet proudového pole •Výstupem jsou parametry proudění (hmotnostní průtoky)
Isight
Fluent
Optimalizační nástroj Gradientová metoda LSGRG2
Sculptor
Deformace sítě pomocí ASD objemů Obr. 3.2 Optimalizační proces
3.7 Výsledky optimalizačního procesu
Výsledkem optimalizačního procesu je zdeformovaná geometrie, která se značně liší od geometrie původní. Toto srovnání je patrné z Obr. 3.3. Je zřejmé, že kanály proti rychlostnímu vstupu se zúžily a kanály po jeho bocích se rozšířily. 8
Obr. 3.3 Geometrie před a po optimalizaci Velice pozitivně se tyto tvarové změny promítly do rozložení hmotnostních toků na výstupech z traktu. Jednotlivé toky před a po optimalizaci jsou zaznamenány v Tab. 3.1 a pro srovnání s průměrnou (optimální) hodnotou hmotnostního průtoku je dále přiložen Graf 3.1. Tab 3.1 Vstup [kg/h]
Výstup1 [kg/h]
Výstup2 [kg/h]
Výstup3 [kg/h]
Výstup4 [kg/h]
Před optimalizací
288,072
64,009
93,932
72,078
58,051
Po optimalizaci
288,072
68,538
80,711
71,974
66,847
Graf 3.1 100
Před optimalizací
90
Optimalizovaná varianta
80
Průměrná hodnota
70 60
4
m& vstup
i =1
4
průměr = ∑
50 40 30
− m& i
20 10 0 Výstup 1
Výstup 2
Výstup 3
Výstup 4
Pro větší názornost je zde třeba ukázat, jak se změnilo proudové pole. To je zobrazeno v Obr. 3.4. z tohoto obrázku je patrné, ve kterých částech geometrie došlo k rozšíření a ve kterých k zúžení. Je patrné, že výstup 2 se velice zúžil. To by mohlo způsobit značnou tlakovou ztrátu. Tím se ale tato úloha nezabývá. V dalších řešeních bude potřebné brát na tento fakt ohled. Dále je vidět, že v tomto místě dochází k značnému urychlení proudu, což by mohlo způsobit poměrně vysokou aerodynamickou hlučnost. Tím by bylo třeba se také zabývat.
9
zúžení Kontury rychlostí, řez výstupem 2
Kontury rychlostí, řez výstupem 3
Kontury tlaků, řez výstupem 3
Kontury tlaků, řez výstupem 2
Pole vektorů rychlostí v traktu topení
Obr. 3.2 Proudové pole zoptimalizované varianty 4. Závěr
Byla navržena úloha tvarové optimalizace na geometrii traktu topení automobilu s cílem dosáhnout rovnoměrnosti hmotnostních průtoků na výstupech z traktu topení. K tomu bylo použito jako nástroje gradientové metody LSGRG2, která byla vyhodnocena jako nejvhodnější optimalizační schéma pro tuto úlohu. Tuto optimalizační úlohu se podařilo úspěšně vyřešit a získat uspokojivé výsledky. Vzhledem k omezením, kterými byla úloha vytyčena (nutnost zachovat tvar pro usazení do automobilu, omezení kvalitou sítě), nebylo možné dosáhnout úplné rovnoměrnosti. Přesto však jsou výsledky velice uspokojivé a pro praxi použitelné. U reálných úloh však není možné řešit jen část traktu topení jako v této práci, ale je nutné řešit rozvody topení jako komplet v celém automobilu, čímž nám přibývají další problémy, které je nutno řešit, jako například vznik protitlaku, s kterým se tato práce nezabývala. 5 Použitá literatura
[1] [2] [3] [4] [5] [6] [7] [8]
Manuály pro užití komerčního software Fluent, 2005 Manuály pro užití komerčního software Tgrid, 2005 Manuály pro užití komerčního software Sculptor, 2005 Manuály pro užití komerčního software ISight, 2005 LEDERER P.: Teorie a optimalizace mechanických systémů. Ediční středisko ČVUT, Praha, 1987 KOZUBKOVÁ, M., DRÁBKOVÁ, S., ŠŤÁVA, P. : Matematické modely stlačitelného a nestlačitelného proudění, Metoda konečných objemů, Ostrava 1999 MUŽÍK T.: Numerická simulace proudění v traktu topení automobilu společnosti ŠkodaAuto, Závěrečný projekt, ČVUT v Praze, Fakulta strojní, 2006 LASDON, L. S. : Optimization Theory for Large Systems, New York, 1972 10
