NUMERICKÁ ANALÝZA PROCESů

1. 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

2. Řešení transportních rovnic NAP10 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 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)

3. w e A W P E B Řešení transportních rovnic NAP10 U 1D problémů je kontrolní objem vlastně jen úsečka Difuzní tok hranicí kontrolního objemu e Konvektivní tok hranicí kontrolního objemu e

4. Řešení transportních rovnic NAP10 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í. 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.

5. Řešení transportních rovnic NAP10 Dosá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 Toto schema splňuje požadavek minimax pro libovolné Reynoldsovo číslo buňky

6. Řešení transportních rovnic NAP10 Naprosté 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, Upwind prvního řádu Upwind druhého řádu 1 f 0 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).

7. Řešení NS rovnic NAP10 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. Beckman

8. Proudová funkce-vířivost NAP10 Metoda proudové funkce a vířivosti (2D proudění nestlačitelné kapaliny) aplikovaná na NS rovnice a rovnici kontinuity 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 Výraz v závorkách je zetová složka  vektoru vířivosti Složky rychlostí u,v můžeme vyjádřit pomocí jediné skalární funkce tak, aby byla automaticky splněna rovnice kontinuity Vířivost lze vyjádřit proudovou funkcí přímo z definice

9. Proudová funkce-vířivost NAP10 Vířivost popisuje rotaci kapaliny a může být charakterizována alternativně vektorem nebo tenzorem 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

10. Proudová funkce-vířivost NAP10 Uvažujme tuhé těleso, rotující kolem osy z úhlovou rychlostí z y r  x vířivost je tedy v každém bodě tuhého rotujícího tělesa konstantní

11. Proudová funkce-vířivost 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) Čá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. 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

12. stěna je proudnice takže i w je konstanta y H x Osa kanálu je proudnice takže třeba =0 Proudová funkce vířivost NAP10 Proudová funkce  Vstup kanálu,  roste od nuly do w třeba =uy Na výstupu předpokládáme např. téměř vyvinutý profil a N-1=N M N-1 N 2 1 2 1 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 Pokud by byl na vstupu rychlostní profil vyvinutý, byl by první člen nula. Když ne, aproximujeme ho rozvojem proudové funkce Stejný postup se aplikuje u stěny kde je předepsaná rychlost UM (v tomto případě nulová)

13. UM M y H N-1 N 2 1 2 x 1 Příklad – tok v dutině 1/5 NAP10 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ý). 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

14. UM y N W P E S x Příklad – tok v dutině 2/5 NAP10 Metoda kontrolních objemů aplikovaná na transport vířivosti (upwind prvního řádu). Uvažujeme čtvercovou síť =x=y Eliptická rovnice pro proudovou funkci 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ů)

15. Metoda střídavých směrůADI 3/5 NAP10 ADI-Alternating Directions Implicit y x-implicitní půlkrok ve kterém se přesně splní okrajové podmínky vlevo i vpravo. 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 y-implicitní půlkrok, ve kterém se přesně splní okrajové podmínky dole i nahoře. výpočet první řady y-implicitní. Postupuje se tentokrát zleva doprava. x

16. Příklad – tok v dutině 4/5 NAP10 Matlab % 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-implicit for 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); end end %stream function y-implicit for 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); end end % vorticity boundary conditions for 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-implicit for 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); end end %vorticity y-implicit for 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); end end end

17. Příklad – tok v dutině 5/5 NAP10 Re=1000 Re=1    u v

18. Primitivní proměnné 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: • Oscilace tlaků a rychlostí, tzv. šachovnicový vzor (checkerboard pattern) způsobený vynecháním centrálních hodnot tlaků při aproximaci gradientů • Jak počítat tlaky z rovnice kontinuity v níž se u nestlačitelné kapaliny tlak vůbec nevyskytuje. Bellin

19. Primitivní proměnné NAP10 Příklad: Stacionární 2D proudění (neznámé jsou rychlostiu,v a tlakp Rovnice prou Rovnice prov To je rovnice prop?

20. 0 0 0 10 10 10 10 10 10 10 10 10 0 0 0 0 0 0 0 0 10 Primitivní proměnné NAP10 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

21. 0 0 0 0 10 10 10 10 10 10 10 10 10 0 0 0 0 0 0 0 10 Primitivní proměnné NAP10 Ú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. Kontrolní objem pro hybnost ve směrux Kontrolní objem pro hybnost ve směru y Kontrolní objem pro rovnici kontinuity

22. Metoda SIMPLE NAP10 SIMPLE neznamená jednoduchá, ale Semi Implicit Method for Pressure Linked Equations. První krok každé iterace spočívá v odhadu tlaku p* 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)

23. Metoda SIMPLE NAP10 “správné” řešení Zanedbáme!!! odečti Dosazeníu*+u’ a v*+v’ do rovnice kontinuity Řešení rovnic pro korekce tlaků

24. Metoda SIMPLE NAP10 Třetí krok: aktualizace tlaků a rychlostí z vypočtených korekcí Tyto zpřesněné tlaky se použijí jako nová výchozí iterace Další iterace (návrat k prvnímu kroku SIMPLE). Iterace se opakují tak dlouho dokud nejsou korekce tlaků dostatečně malé.

25. Metoda SIMPLE NAP10 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).

26. Příklad FLUENT tok v dutině u v NAP10 Re=1 pro porovnání jsou zde uvedeny dříve získané výsledky metodou proudové funkce a vířivosti (Re=1 ale hrubší síť) v u