MODELOWANIE DYNAMIKI PIERŚCIENIA WIROWEGO METODĄ CZĄSTEK WIROWYCH Z WYKORZYSTNIEM OBLICZEŃ RÓWNOLEGŁYCH NA KARTACH GRAFICZNYCH

MODELOWANIE INŻYNIERSKIE ISSN 1896-771X 44, s. 145-150, Gliwice 2012

MODELOWANIE DYNAMIKI PIERŚCIENIA WIROWEGO METODĄ CZĄSTEK WIROWYCH

Z WYKORZYSTNIEM OBLICZEŃ RÓWNOLEGŁYCH NA KARTACH GRAFICZNYCH

ANDRZEJ KOSIOR,HENRYK KUDELA

Instytut Inżynierii Lotniczej, Procesowej i Maszyn Energetycznych, Politechnika Wrocławska e-mail: henryk.kudela@pwr.wroc.pl, andrzej.kosior@pwr.wroc.pl

Streszczenie. W artykule przedstawiono metodę wiru w komórce (VIC – Vortex in Cell) służącą do rozwiązywania równań ruchu płynu. Ze względu na długi czas obliczeń przy użyciu jednego procesora w pracy została przedstawiona implementacja metody VIC z wykorzystaniem obliczeń równoległych w środowisku kart graficznych. Uzyskano dzięki temu 46-krotne skrócenie czasu obliczeń. Przedstawione zostały dwa przypadki testowe. Pierwszy prezentuje poprawność obliczeń prowadzonych z wykorzystaniem kart graficznych. Drugi pokazuje wpływ lepkości na prędkość przemieszczania się pierścienia wirowego.

1. WSTĘP

Znaczenie wirowości dla badań zjawisk hydrodynamicznych jest trudne do przecenienia.

Wszystkie przepływy rzeczywiste mają wirowość różną od zera. Jedną z prostszych, a jednocześnie fascynujących, trójwymiarowych struktur wirowych obserwowanych w eksperymentach, jest pierścień wirowy. Znajduje ona liczne próby wykorzystania w zagadnieniach technicznych, takich jak odwierty podwodne, gaszenie pożarów lub modelowanie zjawisk pogodowych typu szkwał [6]. Wzajemna interakcja pierścieni wirowych prowadzi do ciekawych zjawisk hydrodynamicznych, takich jak gra wirów (wzajemne przeciąganie się pierścienie wirowych) czy też rekonekcja (zmiana topologii dwóch zderzających się pierścieni). Do modelowania numerycznego ruchu pierścienia wirowego wybrano metodę cząstek wirowych. Jest to metoda wyjątkowo dobrze nadająca się do modelowania zjawisk wirowych, a także bardzo dobrze nadaje się do prowadzenia obliczeń w środowisku wieloprocesorowym. W obliczeniach wykorzystuje się cząstki, które są nośnikami wirowości. Cząstki wirowe umożliwiają w łatwy i efektywny sposób analizę ewolucji pól wirowych. Rozwiązywanie równań ruchu cieczy dla zagadnień trójwymiarowych, niezależnie od użytej metody, wiąże się z długimi czasami obliczeń.

Przyrost mocy obliczeniowej komputerów, związany ze zwiększeniem częstotliwości zegara taktującego pracę procesora, w ostatnich latach zdecydowanie się zmniejszył. Dlatego dla przyspieszenia obliczeń zmuszani jesteśmy do stosowania wieloprocesorowych środowisk obliczeniowych i opracowywania algorytmów właściwych dla obliczeń równoległych.

Niespodziewanie okazało się, że budowane dla gier komputerowych karty graficzne, posiadające niekiedy kilkaset procesorów strumieniowych, można wykorzystać do obliczeń

naukowych. W pracy przedstawiono wyniki numeryczne dla ruchu pierścienia wirowego dla cieczy doskonałej (nielepkiej) oraz badania dotyczące wpływu lepkości na prędkość przemieszczania się pierścienia.

2. RÓWNANIA RUCHU

Równania nieściśliwego, lepkiego płynu można zapisać w następujący sposób:



u



u u

u     



 

 p t

1 (1)

0





 u (2)

gdzie u



u ,,v w



jest wektorem prędkości,  - gęstością płynu, p - ciśnieniem,  - kinematycznym współczynnikiem lepkości, w równaniu (1) działając obustronnie operatorem



 można przekształcić je w równanie Helmholtza dla ewolucji wirowości :



u



ω



ω



u ω

ω     



 

t (3)

gdzie ωu. Z warunku nieściśliwości (2) wynika istnienie wektora potencjału A A

u (4)

Zakładając dodatkowo, że wektor A jest nieściśliwy



 A 0



, jego składowe można wyznaczyć, rozwiązując równanie Poissona

3 , 2 , 1

, 





A_i _i i

(5) Rozwiązując równanie (5), można obliczyć prędkość ze wzoru (4).

W metodzie cząstek wirowych używany jest algorytm dekompozycji lepkościowej.

Rozwiązanie uzyskuje się w dwóch krokach: najpierw rozwiązywane jest równanie Eulera (dla cieczy nielepkiej).



u



ω



ω



u ω   





t (6)

W drugim kroku symulowany jest efekt lepkości poprzez rozwiązanie równania dyfuzji.

ω ω



 

 

t (7)

Do rozwiązania równania dyfuzji można użyć różnych metod, takich jak Particle Strength Exchange (PSE), metodę Monte Carlo lub metody różnic skończonych. W obecnej pracy zastosowano metodę różnic skończonych ze schematem niejawnym do przybliżonego rozwiązania równania (7).

3. OPIS METODY WIR W KOMÓRCE DLA PRZYPADKU CIECZY NIELEPKIEJ

W pracy została wykorzystana metoda dekompozycji lepkościowej. W metodzie tej rozwiązanie znajdowane jest w dwóch podkrokach. Najpierw rozwiązywane są równania ruchu dla płynu nielepkiego. W tym celu obszar obliczeniowy został pokryty trójwymiarową

siatką strukturalną



^j₁^x,j₂^{y, j}₃^z

 

^j₁^{, j}₂^{, j}₃ ^1,2,...,N



^{, gdzie xyz h. Ta sama}

siatka będzie wykorzystywana do rozwiązania równania Poissona. Ciągłe pole wirowości zostało zastąpione dyskretnym rozkładem delt Diraca. [1, 4]

    

p p



N

p

p x x x

α x

ω 









1 (8)

gdzie α_p oznacza masę (cyrkulację) cząstki wirowej α_p 



_p1,_p2,_p3



w punkcie



_p1, _p2, _p3



p  x x x

x . i-ta składowa  jest zdefiniowana przez wyrażenie: _i



x1,x2x3



d h³ _i

 

_p , _p V_p, V_p h³

V i i

p











^{ x  x x}



(9) gdzie V_pjest objętością komórki z indeksem p.

Z twierdzenia Helmholtza [7] wiadomo, że wirowość jest unoszona przez płyn, a więc równanie ruchu cząstek wirowych ma postać:



^t



dt d

p

p u x ,

x 

(10) Należy także wziąć pod uwagę fakt, że w trójwymiarowym polu wirowym intensywność cząsteczek jest zmieniana przez efekt zwany rozciąganiem linii wirowych (prawa strona równania (6). A więc:

 



p



p

p t

dt

dα u x α





 ,

(11) Prędkość w węzłach została obliczona w wyniku rozwiązania równania Poissona (5) metodą różnic skończonych przy użyciu (4). Układ równań (8), (9) został rozwiązany metodą Rungego–Kutty 4. rzędu.

W drugim podkroku algorytmu symulowany był wpływ lepkości na cząsteczki poprzez rozwiązanie równania:

ω ω



 

 

t (12)

4. REDYSTRYBUCJA MASY CZĄSTEK WIROWYCH NA WĘZŁY SIATKI

W metodzie wirów w komórce cząsteczki mają tendencję do zbierania się w obszarach dużych gradientów prędkości. Może to prowadzić do niedokładności spowodowanych zbytnim zbliżaniem się cząsteczek do siebie. Aby temu zapobiec, cząsteczki umieszcza się z powrotem na węzłach siatki numerycznej (tzw. remeshing). Wykonuje się to przy pomocy interpolacji:

3

~ ~ 







 





^j _h ^{p h}

p pn j

x

 x





(13) gdzie j jest indeksem węzła siatki, a p jest indeksem cząsteczki. Jakość interpolacji zależy od własności jądra φ. W tej pracy zostało wykorzystane następujące jądro interpolacyjne:

 

   































x x

x x

x x x

2 2 1

1 0

0

2 / 1 2

2 / 3 5 2

2 2 3



(14) Jądro (14) jest rzędu m = 4. W przypadku trójwymiarowym 



x,y,z





     

x y z . Proces redystrybucji musi być dokonywany w każdym kroku czasowym, aby możliwe było rozwiązanie równania Poissona (5).

4. REALIZACJA NA GPU

Z powodu długiego czasu obliczeń powstała potrzeba zastosowania obliczeń równoległych. Wybrana została architektura CUDA wykorzystująca karty graficzne (GPU - Graphics Processing Unit). Szczegóły implementacji znajdują się w [3]. W celu sprawdzenia poprawności implementacji porównano otrzymane wyniki obliczeń na pojedynczym procesorze (CPU - Central Processing Unit) i na karcie graficznej. Wyniki wykazały poprawność implementacji. Na GPU obliczenia były 46 razy szybsze niż na CPU.

Rys.1. Porównanie wyników dla ewolucji pierścienia wirowego po n = 100 krokach czasowych dla obliczeń prowadzonych na CPU (Intel Core i7 960 - góra) i GPU (NVIDIA

GTX480 - dół). Obliczenia na GPU były 46 razy szybsze.

5. PRZEPŁYW NIELEPKI I LEPKI

Przeprowadzono badania symulujące wpływ lepkości na zachowanie się ewolucji pierścienia wirowego. Liczbę Reynoldsa przepływu definiowano jako Re/ . Ewolucja pierścienia wirowego w przepływie nielepkim jest widoczna na rys. 1. a w przepływie lepkim na rys. 2. Na rys. 4. i 5. widoczne są odpowiednio przemieszczenie i prędkość pierścienia wirowego w czasie z zależności od lepkości.

Z zaprezentowanych wyników widać, że prędkość translacji pierścienia wirowego w płynie lepkim maleje z czasem. Jest to zgodne z wynikami teoretycznymi [2]. Dla przepływu nielepkiego prędkość translacji pierścienia wirowego jest stała, co widać na rys. 4. Pokrywa się ona ze wzorem teoretycznym Kelvina [5]. Wahania w prędkości przemieszczania wynikają z ustalania się rozkładu wirowości wewnątrz rdzenia pierścienia wirowego. Do obliczeń przyjęto równomierny rozkład wirowości. Wraz z upływem czasu rozkład wirowości wewnątrz rdzenia przyjmuje rozkład podobny do rozkładu Gaussa.

Rys.2. Ewolucja pierścienia wirowego w przepływie nielepkim dla Γ = 1.0. Dół - t =

0, góra - t = 10

Rys.3. Ewolucja pierścienia wirowego w przepływie lepkim (Re = 1000) dla Γ = 1.0.

Dół - t = 0, góra - t = 10

Rys.4. Zależność położenia środka pierścienia wirowego od czasu dla przepływu

nielepkiego i lepkiego

Rys.5. Zależność prędkości środka pierścienia wirowego od czasu dla przepływu

nielepkiego i lepkiego

6. WNIOSKI

Z zaprezentowanych wyników można wysnuć wnioski, że metoda cząstek wirowych jest odpowiednia do badania zagadnień związanych z ewolucją trójwymiarowych struktur wirowych. Mają one fundamentalne znaczenie dla zrozumienia procesów powstawania i rozwoju turbulencji. Zbudowany program obliczeniowy pozwala na symulację trójwymiarowych przepływów lepkich. W pracy wykorzystano karty graficzne do obliczeń równoległych. Pozwoliło to na uzyskanie 46-krotnego przyspieszenia obliczeń w stosunku do pojedynczego rdzenia. Rozkład wirowości w przekroju pierścienia wirowego jest ważnym czynnikiem mającym wpływ na jego zachowanie.

LITERATURA

1. Cottet G. H., Koumoutsakos P. D.: Vortex Methods: Theory and practice. Cambridge University Press, 2000, p. 16 – 25.

2. Kaplansky F.B., Rudi Yu. A.: Evolution of a viscous vortex ring. “Fluid Dynamics”

2001, Vol. 36, No. 1.

3. Kosior A., Kudela H.: Parallel computations on GPU in 3D using the vortex particle method. http://dx.doi.org/10.1016/j.compfluid.2012.01.014,

4. Kudela H., Regucki P.: The vortex-in-cell method for the study of three-dimensional flows by vortex methods. Tubes, Sheets and Singularities in Fluid Dynamics. „Fluid Mechanics and Its Applications”. Kluwer Academic Publisher, 2009, Vol. 7, p. 49 – 54.

5. Lim T. T., Nickels T. B.: Vortex rings. Fluid vortices. Editor: Sheldon I. Green.

Dordrecht: Kluwer Academic Publishers, 1995.

6. Shariff K., Leonard A.: Vortex rings. “Annu. Rev. Fluid Mech” 1992, Vol. 24, p. 235- 279.

7. Wu J. Z., Ma H. Y., Zhou M. D.: Vorticity and vortex dynamics. Berlin, Heidelberg:

Springer, 2006.

Zadanie współfinansowane ze środków Unii Europejskiej w ramach Europejskiego Funduszu Społecznego

MODELING OF THE VORTEX RING DYNAMICS BY VORTEX PARTICLE METHOD

USING PARALLEL COMPUTATION ON GRAPHICS CARDS

Summary. The paper presented the Vortex in Cell (VIC) method for solving the fluid motion equations in 3D. Due to the long time of computation on single processor the parallel implementation of the VIC method was presented. The speed-up for the entire VIC method implementation on the GPU was 46 times.

Two test cases were presented. First one, shows correctness of the parallel implementation on GPU. Second example shows influence of viscosity on the translation velocity of the vortex ring.