72 | Pomurska obzorja 1/ 2014/ 2 Tehnika Zoran Žunič * Mešana metoda robnih in končnih elementov za reševanje hitrostno vrtinčne formulacije Navier-Stokesovih enačb 1. Uvod Pri razvoju na področju numeričnih metod za obravnavanje tokov viskoznih tekočin se običajno osredotočamo na aproksimacijske metode za reševanje Navier-Stokesovih enačb. Večina teh pristopov uporablja le eno od aproksimacijskih metod, na primer metodo končnik razlik, metodo končnih prostornin, metodo končnih elementov ali metodo robnih elementov. Vsaka od teh metod ima svoje prednosti in slabosti, zato se je začel razvoj tudi na področju istočasne uporabe večih metod za reševanje istega sistema enačb. Z uporabo hitrostno-vrtinčne formulacija so se v preteklosti ukvarjali že Liu [1], ki je uporabil metodo končnih razlik, Guevremont in dr. [2], [3] ter Wong in Baker [4], ki so uporabili metodo končnih elementov. Na področju robnih elementov so delovali Žagar in Škerget [5], Hriberšek in Škerget [6], Ramšak in dr. [7] ter Ravnik in dr. [8]. Nekaj raziskovalcev se lotilo tudi uporabe večih metod, med katere spadajo Young in dr. [9] ki so uporabljali metodo robnih elementov in metodo končnih elementov vendar s formulacijo z osnovnimi spremenljivkami, ter Brown in Ingber [10], ki sta notranje hitrosti računala z posplošeno Helmhotzovo metodo dekompozicije, transport vrtinčnosti pa z Galerkinovo metodo končnih elementov. Prikazano delo je kombinacija reševanja vektorske Poissonove hitrostne enačbe z metodo robnih elementov za izračun robnih vrtinčnosti in uporabe metode končnih elementov za reševanje vektorske Poissonove hitrostne enačbe za notranje hitrosti in transportne enačbe vrtinčnosti. Prvi pristop omogoča implicitni izračun robnih vrtinčnosti, medtem ko drugi pristop reši težavo s prevelikimi zahtevami po računalniškem pomnilniku. 2. Hitrostno vrtinčna formulacija Tok gibanja zvezne tekočine opisujeta enačbi ohranitve gibalne količine in ohranitve mase, ki sta zapisani z osnovnimi spremenljivkami (tlaki in hitrostmi). Pri reševanju problemov nestisljivih tokov, pri katerih predpostavimo, da je gostota konstantna, se pojavi težava zaradi odsotnost tlaka in gostote v enačbi ohranitve mase. Težavo rešimo z vpeljavo pojma vrtinčnosti, 0 ,              v (1) s pomočjo katerega preoblikujemo osnovne enačbe v hitrostno – vrtinčno obliko (Škerget in dr. [7], Ravnik in dr. [8]). Dobimo Poissonovo enačbo za hitrost (enačbo kinematike) 0 2          v , (2) in ob predpostavki Newtonske tekočine še transportno enačbo za vrtinčnost (enačbo kinetike)              2 ) ( ) (           v v t . (3) Iščemo rešitev enačb (2) in (3), ki zadovoljuje začetne in robne pogoje v računskem območju. 3. Metoda robnih elementov Enačbo (2) zapišemo kot singularno robno integralsko enačbo za vektor hitrosti s pomočjo Greenovih teoremov (Wrobel [11])                   d u d v n u d u n v v c * * * ) ( ) ( ) ( ) ( ) (             (4) POVZETEK V prispevku je predstavljena numerična metoda za reševanje nestisljivega toka viskozne newtonske tekočine z uporabo hitrostno- vrtinčne formulacije ohranitvenih enačb. Izpeljali smo robno integralsko enačbo kinematike v tangentni obliki, s čemer je bil dosežen natančnejši izračun robnih vrednosti vrtinčnosti. Po metodi končnih elementov smo izpeljali še integralsko enačbo kinematike, s katero smo izračunali neznane območne hitrosti in pa enačbo kinetike, s katero smo izračunali neznane območne vrtinčnosti. Vse dobljene sisteme diskretnih enačb smo povezali v računski algoritem, ki z metodo Picardove iteracije nesklopljeno reši vezan sistem nelinearnih enačb. Da bi pokazali robustnost in uporabnost predstavljene metode in istočasno izločili možnost napak v programu, smo uporabili primer gnanega toka v prostorski kotanji. Ključne besede: prostorski tok viskozne tekočine, Navier-Stokesove enačbe, hitrostno-vrtinčna formulacija, metoda robnih elementov, metoda končnih elementov. * AVL – AST d.o.o. / Ulica kneza Koclja 22, SI-2000 Maribor E-naslov: zoran.zunic@amis.net Dr. Zoran Žunič, Lorbekova 15, 2341 Limbuš, Slovenia; Tel.: +386-2-250-85-16 Zoran ŽUNIČ: MEŠANA METODA ROBNIH IN KONČNIH ELEMENTOV ZA REŠEVANJE… Pomurska obzorja 1/ 2014/ 2 | 73 kjer * u predstavlja eliptično Laplaceovo osnovno rešitev,  je izvorna točka na robu računskega območja, c je geometrijski koeficient in n  normala na rob računskega območja. Enačbo (4) z Gausovim divergenčnim teoremom pretvorimo v tangentno obliko, s čimer se znebimo odvodov hitrostnega in vrtinčnega polja (Škerget in dr. [7], Ravnik in dr. [8]), ki jih prenesemo na osnovno rešitev                    d u d v n u d v n u v c * * * ) ( ) ( ) ( ) (             (5) 4. Metoda končnih elementov Z uporabo Galerkinove metode utežnih ostankov in uporabo Greenovega prvega teorema v enačbi (2) dobimo šibko obliko Poissonove enačbe za hitrostni vector 0 ) ( ) ( ) (                     d N d v n N d v N          (6) kjer nam N predstavlja utežno funkcijo. Enak postopek ponovimo z enačbo (3), da dobimo šibko obliko prenosne enačbe za vrtinčnost                                     d N t d n N d N d v N d v N d N t 1 1 ) ( ) ( ) ( ) ( 1                        (7) kjer je t  velikost časovnega koraka in podpis 1   vrednost v prejšnjem časovnem koraku. Velikosti integralov potrebnih za reševanje enačb (6) in (7) so veliko manjše od velikosti integralov v enačbi (5), kar nam omogoča hitrejše reševanje problemov pri enaki gostoti računske mreže, oziroma obravnavanje večjih (kompleksnejših) problemov. 5. Računski postopek Računski postopek je potekal v naslednjem vrstnem redu: 1. Izberemo začetno hitrostno polje in izračunamo začetno vrtinčno polje z enačbo (1). Začetni čas =0. 2. Časovna zanka. 3. Zanka zaradi nelinearnosti. 4. Kinematika toka: a. izračunamo robne vrtinčnosti (5) z metodo robnih elementov, b. izračunamo notranje hitrosti (6) z metodo končnih elementov. 5. Kinetika toka, transport vrtinčnosti: a. izračunamo notranje vrtinčnosti (7) z metodo končnih elementov, b. podrelaksiramo rešitev. 6. Preverimo napako rešitve, če je večja od tolerance gremo na korak 3. 7. Konec časovnega koraka: a. shranimo vrednosti hitrosti in vrtinčnosti, b. če je čas manjši od predpisanega gremo na korak 2. 8. Konec izračuna. 6. Primerjava rezultatov Gnan tok v prostorski kotanji je eden najpogosteje uporabljanih testnih primerov za testiranje laminarnih tokov viskozne Newtonske tekočine. Kljub enostavni geometriji, se v njem pojavljajo vsi fenomeni kompeksnih tokov, kot so vrtinci, sekundarni vrtinci, prostorski vzorci, kaotično gibanje, nestabilnosti in turbulence (Shankar in Despande [12]). V našem primeru smo uporabili kotanjo v obliki kocke s stranicami dolžine 1 enote. Reynoldsovo število, ki je v tem primeru edino karakteristično število, je bilo zasnovano na dolžini stranice kocke in hitrosti vlečenja zgornje ploskve (Re = vx L/ ) in je znašalo Re=100 in 1000. Velikost uporabljenih mrež je znašala 8 3 , 16 3 in 32 3 elementov. Vse mreže so bile zgoščene proti zunanjim stenam s faktorjem 8 (slika 1a). Predpisali smo naslednje robne pogoje:  zgornja ploskev (𝑧 = 1): 𝑣 𝑥 = 1, 𝑣 𝑦 = 𝑣 𝑧 = 0,  ostale ploskve: 𝑣 𝑥 = 𝑣 𝑦 = 𝑣 𝑧 = 0. Vse testne primere smo izračunali z začetnimi pogoji 𝑣 𝑥 = 𝑣 𝑦 = 𝑣 𝑧 = 0. Slika 1. a) Računska mreža, b) Hitrostni profili za Re=1000. Zoran ŽUNIČ: MEŠANA METODA ROBNIH IN KONČNIH ELEMENTOV ZA REŠEVANJE… 74 | Pomurska obzorja 1/ 2014/ 2 Slika 2. Izopovršine hitrosti |v  |=0.13. Levo Re=100, desno Re=1000. Na sliki 1b je prikazana primerjava profilov hitrosti, dobljenih na različno gostih mrežah, z rezultati Yanga in dr. [13]. Primerjali smo hitrostne profile vx in vy po navpični in vodoravni srednjici kocke. Analiza z zgoščevanjem mreže je pokazala, da je najredkejša mreža pregroba, medtem ko se gostejši mreži dobro ujemata z referenčnimi. Slika 2 prikazuje izopovršine absolutne hitrosti |v  |=0.13 za različne vrednosti Reynoldsovega števila. Vidimo, da se hitro vrteče se jedro toka pomika proti zunanjim stenam, obenem pa se debelina hitrega pasu zmanjšuje. Slika 3 prikazuje prikazuje vrtinčno polje na različnih presekih kotanje (x = 0.5, y = 0.5 in z = 0.5) za vrednost Reynoldsovega števila Re = 1000. Iz izolinij vrtinčnosti je razvidno, da je tokovno polje pri izbrani vrednosti Reynoldsovega števila simetrično preko xz ravnine. Slika 3. Izolinije vrtinčnosti za Re=1000. Zgoraj levo  x, zgoraj desno  y, spodaj  z . Zoran ŽUNIČ: MEŠANA METODA ROBNIH IN KONČNIH ELEMENTOV ZA REŠEVANJE… Pomurska obzorja 1/ 2014/ 2 | 75 7. Zaključki V prispevku smo prikazali postopek reševanja Navier- Stokesovih enačb zapisanih s hitrostno-vrtinčno formulacijo. Robne vrtinčnosti smo izračunali implicitno z uporabo metode robnih elementov, medtem ko smo notranje hitrosti in notranje vrtinčnosti izračunali z uporabo metode končnih elementov. Na tak način smo ohranili visoko natačnost izračuna robnih vrtinčnosti, prihranili pa smo pri količini potrebnega računalniškega pomnilnika. Za preverjanje pravilnosti in prikaz natančnosti uporabljene metode smo uporabili primer gnanega toka v prostorski kotanji. Primerjava dobljenih rezultatov z rezultati iz literature je pokazala odlično ujemanje. Literatura 1. C. H. Liu, Numerical solution of three-dimensional Navier- Stokes equations by a velocity-vorticity method, Int. J. Numer. Meth. Fluids 35, 533–557, 2001. 2. G. Guevremont, W. G. Habashi, Finite element solution of the Navier-Stokes equations by a velocity-vorticity method, Int. J. Numer. Meth. Fluids 10, 461–475, 1990. 3. G. Guevremont, W. G. Habashi, P. L. Kotiuga, M. M. Hafez, Finite element solution of the 3D compressible Navier-Stokes equations by a velocity-vorticity method, Journal of Computational Physics 107 (1), 176–187, 1993. 4. K. L. Wong, A. J. Baker, A 3D incompressible Navier- Stokes velocity-vorticity weak form finite element algorithm, Int. J. Numer. Meth. Fluids 38, 99–123, 2002. 5. I. Žagar, L. Škerget, Boundary elemens for time dependent 3-D laminar viscous fluid flow, Journal of mechanical engineering 35 (10/12), 160–163, 1989. 6. M. Hriberšek, L. Škerget, Iterative methods in solving Navier-Stokes equations by the boundary element method, Int. J. Numer. Methods Engnrg 39, 115–139, 1996. 7. M. Ramšak, L. Škerget, M. Hriberšek, Z. Žunič, A multidomain boundary element method for unsteady laminar flow using stream function – vorticity equations, Eng. Anal. Bound. Elem. 29, 1–14, 2005. 8. J. Ravnik, L. Škerget, M. Hriberšek, Two-dimensional velocity-vorticity based LES for the solution of natural convection in a differentially heated enclosure by wavelet transform based BEM and FEM, Eng. Anal. Bound. Elem. 30, 671–686, 2006. 9. D. L. Young, J. L. Huang, T. I. Eldho, Simulation of laminar vortex shedding flow past cylinders using a coupled BEM and FEM model, Comput. Methods Appl. Mech. Engrg. 190, 5975–5998, 2001. 10. M. J. Brown, M. S. Ingber, Parallelization of a vorticity formulation for the analysis of incompressible viscous fluid flows, Int. J. Numer. Meth. Fluids 39, 979–999, 2002. 11. L. C. Wrobel, The Boundary Element Method, John Willey & Sons, LTD, 2002. 12. P. N. Shankar, M. D. Deshpande, Fluid mechanics in the driven cavity, Annu. Rev. Fluid Mech. 32, 93–136, 2000. 13. J. Y. Yang, S. C. Yang, Y. N. Chen, C. A. Hsu, Implicit weighted ENO schemes for the three-dimensional incompressible Navier-Stokes equations, Journal of Computational Physics 146, 464–487, 1998.