MODELOWANIE INŻYNIERSKIE ISSN 1896-771X 44, s. 245-251, Gliwice 2012
BADANIE DYNAMIKI PIERŚCIENIA WIROWEGO ZABURZONEGO FALAMI KELVINA
O SKOŃCZONEJ AMPLITUDZIE
Z WYKORZYSTANIEM TRÓJWYMIAROWEJ METODY „WIR W KOMÓRCE”
PAWEŁ REGUCKI
Instytut Inżynierii Lotniczej, Procesowej i Maszyn Energetycznych, Politechnika Wrocławska e-mail: pawel.regucki@pwr.wroc.pl
Streszczenie. Jedno z fascynujących i wciąż badanych zjawisk związane jest z ewolucją fal Kelvina (sinusoidalnych zaburzeń) występujących na obwodzie pierścienia wirowego. Numeryczne wyniki prezentowane przez Kiknadze i Mamaladze [4] oraz Barenghi et al. [1] pokazują, że prędkość translacji pierścienia wirowego w warunkach panujących w nadciekłym helu He4 (brak lepkości oraz infinitezymalnie mały promień rdzenia pierścienia - nić wirowa) zależy od częstotliwości i amplitudy fal Kelvina. Interesująca jest więc weryfikacja, czy podobne zachowanie można zaobserwować w granicy klasycznej cieczy nielepkiej. Wyniki numerycznego modelowania prezentowane przez autora odnoszą się do przypadku klasycznego, nielepkiego pierścienia wirowego ze skończonym promieniem wewnętrznym rdzenia. Ewolucja pierścienia była modelowana z zastosowaniem trójwymiarowej metody „wir w komórce” [2].
W prezentowanych wynikach fale Kelvina rozwijają się wzdłuż obwodu pierścienia, prowadząc do skomplikowanej ewolucji pola wirowości. Pomimo złożonej dynamiki wyniki numeryczne wskazują, że prędkości translacji zarówno pierścienia zaburzonego jak i niezaburzonego falą Kelvina pozostają takie same.
1. WSTĘP
Pierścienie są najprostszymi strukturami wirowymi, które można w relatywnie prosty sposób wygenerować w trójwymiarowych przepływach [2, 6, 11]. Pomimo ich geometrycznej prostoty dynamika tych struktur jest interesującym przykładem nieliniowej interakcji obszarów o skoncentrowanej wirowości. Jedno z fascynujących i wciąż badanych zjawisk odnosi się do ewolucji pierścienia wirowego z występującymi na jego obwodzie falami Kelvina (sinusoidalnymi zaburzeniami) [8]. Wyniki numeryczne prezentowane w pracach Kiknadze i Mamaladze [4] oraz Barenghi et al. [1] pokazały, że prędkość przemieszczania pierścienia wirowego zaburzonego falami Kelvina w supercieczy silnie zależy od częstotliwości i amplitudy tych fal. Wzrost wartości amplitudy fali powodował, że pierścień wirowy zwalniał, a nawet, dla dostatecznie dużej jej wartości, mógł zmienić kierunek ruchu na przeciwny. Obliczenia zostały przeprowadzone przy założeniu przepływu nielepkiego z infinitezymalnie małym promieniem wewnętrznym pierścienia, co odzwierciedla warunki
panujące w nadciekłym He4 [3] (brak lepkości i średnica nici wirowej rzędu angstremów).
W celu wyznaczenia pola prędkości generowanego przez włókno wirowe Kiknadze i Mamaladze [4] zastosowali metodę przybliżenia lokalnej indukcji (LIA - local induction approximation), zaś Barenghi et al. [1] zastosował prawo Biota-Savarta (BSL).
Interesująca była zatem weryfikacja, czy podobne zachowanie można zaobserwować w granicy klasycznego, nielepkiego pierścienia wirowego ze skończoną wartością promienia rdzenia.
2. TRÓJWYMIAROWA METODA „WIR W KOMÓRCE”
Równania opisujące ewolucję pola wirowości w trójwymiarowym, nielepkim, nieściśliwym przepływie mają postać [2, 7, 9]:
+ u = u
t
, (1)
u = 0
(2)
gdzie ( 1, 2, 3) oznacza pole wirowości, zaś u = (u , u , u ) 1 2 3 - pole prędkości. Warunek nieściśliwości płynu (2) gwarantuje istnieje potencjału wektorowego A
[7]:
u = A
. (3)
Składowe pola wektorowego A
wyznaczane są z rozwiązania równania Poissona przy dodatkowym założeniu, że A = 0
:
Ai i, i 1, 2,3.
(4)
W metodzie “wir w komórce” ciągłe pole wirowości jest zastępowane dyskretnym rozkładem cząstek wirowych niosących informacje o trzech składowych pola wirowości [2, 5]:
1
(x) (x ) (x x )
N
p p p
p
, (5)
gdzie ( )x jest trójwymiarową deltą Diraca, a p oznacza intensywność cząstki wirowej
1, 2, 3
p p p p p
p=1,...,N w położeniu xp xp
x , x , xp1 p2 p3
. Obszar obliczeniowy pokryty jest siatką numeryczną
NxNyNz
o równomiernym kroku przestrzennym h, zaś i-ta składowa wektora p jest definiowana wyrażeniem
xpVp, Vp h3
:
1 2 3
3(x ) x , x , x x (x )
p
i p i i p
V
d h
(6)Metoda “wir w komórce” pozwala na rozwiązanie równań Eulera dla trójwymiarowego, nielepkiego przepływu. Z teorii Helmholtza wynika, że wirowość unoszona jest przez ciecz, a zatem ewolucja cząstek wirowych odbywa się zgodnie z równaniem:
xp u(x , ),p
d t
d t
(7)
u(x , ) .
p
p p
d t
d t
(8)
Prawa strona równania (8) może być wyrażona na podstawie tożsamości wektorowej
u
u
u
T0.5 u
u
T
. W obliczeniach zastosowano
wyrażenie u
T
, ponieważ lepiej zachowuje niezmienniki ruchu dla przepływu
nielepkiego [2].
T
, ponieważ lepiej zachowuje niezmienniki ruchu dla przepływu nielepkiego [2].W celu rozwiązania równania Poissona (4) intensywności cząstek p są redystrybuowane na węzły siatki numerycznej (k, l, m) z wykorzystaniem metody „objętości węzła” [2]:
,( ) x
dla 0
x
0 dla 0
n n
p i klm p
p n
n n klm
i klm klm
n klm
J J
J
, (9)
gdzie Jklmn
ph3klmn
xp jest nazywany objętością węzła, a górny indeks „n” określa krok czasowy tn n t . Zastosowana metoda redystrybucji jest mniej czuła na położenie cząstki w oczku siatki numerycznej i lepiej przeciwdziała samorzutnemu grupowaniu się cząstek wirowych w obszarach o dużych gradientach prędkości. Jako funkcję wybrano trójwymiarową funkcję sklejaną trzeciego stopnia klm
x k
x l y m
z. Jednowymiarowa funkcja sklejana trzeciego stopnia ma postać:
3 2
3 2
1 2
x x dla x 1,
2 3
1 4
x x x 2 x dla 1 x 2 ,
6 3
0 dla x 2.
(10)
Podsumowując, można stwierdzić, że obliczenia w jednym kroku czasowym tn n t przebiegały następująco [10]:
- redystrybucja intensywności cząstek wirowych na węzły siatki numerycznej zgodnie z równaniem (9),
- rozwiązanie równania Poissona (4) z okresowymi warunkami brzegowymi oraz wyznaczenie w węzłach siatki pola prędkości z wyrażenia (3),
- interpolacja prędkości z węzłów siatki na położenia cząstek wirowych z wykorzystaniem interpolacji Lagrange’a drugiego rzędu oraz wyznaczenie nowych pozycji cząstek zgodnie z równaniem (7) z wykorzystaniem w tym celu metody Rungego-Kutty drugiego rzędu,
- uaktualnienie intensywności cząstek w nowych położeniach zgodnie z równaniem (8).
3. WYNIKI OBLICZEŃ NUMERYCZNYCH
Jako obszar obliczeniowy wybrano sześcian o wymiarach 101010 pokryty prostokątną siatką o równomiernym kroku h = 0.1 w każdym kierunku (rys. 1A). Wybór takiej wartości h związany był z optymalizacją czasu obliczeń numerycznych. W pracy [10] przeprowadzono analizę wpływu kroku siatki h na modelowanie złożonej dynamiki struktur wirowych (np. gry wirów zaprezentowanej w [5]), uzyskując bardzo dobrą zgodność otrzymanych obliczeń z wynikami literaturowymi [2]. Ponadto uśredniający charakter metody „objętości węzła”
wymagał, aby w elementarnym sześcianie h3 znajdowało się od kilku do kilkunastu cząstek wirowych, co było bardzo dobrze spełnione dla h = 0.1. Krok czasowy przyjęto t = 0.02.
Pierścień wirowy został podzielony na 100 przekroi (rys. 1B) a w każdym przekroju rozmieszczono równomiernie 121 cząstek wirowych (rys. 1C). Finalnie, pojedynczy pierścień był aproksymowany 12100 cząstkami [10].
Rys.1. A) Obszar obliczeniowy z pierścieniem wirowym; B) podział pierścienia na 100 przekroi; C) początkowy rozkład 121 cząstek wirowych w pojedynczym przekroju Najpierw zamodelowano ruch nielepkiego, niezaburzonego pierścienia wirowego w celu porównania wartości prędkości translacji U otrzymanych w wyniku obliczeń numerycznych z formułą teoretyczną UT [7] (gdzie a R, - wewnętrzny i zewnętrzny promień pierścienia, - cyrkulacja):
8 1
4 ln 4
T
U R
R a
. (11)
Otrzymano bardzo dobrą zgodność wyników w szerokim zakresie parametrów a R , , , które zestawiono w tabeli 1.
Tabela 1. Porównanie prędkości translacji z obliczeń numerycznych U z formułą (11)
a R
U UT
0.30 2.0 0.80 0.110 0.120
0.30 1.5 1.0 0.181 0.185
0.30 1.0 1.6 0.394 0.390
0.20 2.0 0.80 0.120 0.133
0.20 1.5 1.0 0.215 0.217
0.20 1.0 1.6 0.430 0.441
Następnie obwód pierścienia wirowego został zaburzony N falami Kelvina o amplitudzie A (długość fali zaburzenia wynosi 2R/N) [1]. Początkowe położenie zaburzenia opisano we współrzędnych walcowych
r, , z
:
x cos cos cos( )
y sin cos sin( )
z sin
R A N
R A N
A N
(12)
Obliczenia przeprowadzono dla szerokiego zakresu parametrów N = 4, 6, ..., 14 i stosunku A/R = 0.04, 0.08, 0.12, ..., 0.20. Przykładowe wyniki obliczeń zostały przedstawione na rys. 2 dla parametrów pierścienia wirowego R2.5,a0.2, 1.5.
Rys.2. Widok w rzucie z boku i przodu położenia początkowego i końcowego tf = 10 (odpowiednio czerwony (ciemny) i niebieski (jasny) kolor): A, D) niezaburzonego pierścienia
wirowego; B, E) oraz C, F) zaburzonego pierścienia wirowego dla parametrów odpowiednio:
A/R = 0.04, N = 4 oraz A/R = 0.04, N = 14
Przedstawiony na rys. 2 pierścień wirowy został zaburzony falami Kelvina o parametrach odpowiednio: N = 4, A/R = 0.04 oraz N = 14, A/R = 0.04 (odpowiednio na rys. 2 B oraz C).
Wyniki obliczeń numerycznych zaburzonych pierścieni zestawiono z jego niezaburzonym odpowiednikiem (rys. 2 D). W trakcie ewolucji początkowo małe zaburzenia rozwijają się, deformując kształt pierścienia wirowego, co przedstawiono na rys. 2 E,F. Pomimo skomplikowanej dynamiki końcowe położenia pierścieni: zaburzonego i niezaburzonego po czasie t = 10 są identyczne, niezależnie od doboru początkowych parametrów zaburzenia.
Ewolucję izopowierzchni stałej wirowości dla pierścienia wirowego zaburzonego falą Kelvina o parametrach A/R = 0.04, N = 4 przedstawiono na rys. 3. Można zauważyć, że wraz z upływem czasu początkowo sinusoidalne zaburzenie zaczyna w znaczący sposób deformować kształt pierścienia wirowego, prowadząc do bardzo złożonego dynamicznie pola wirowości.
Rys.3. Widok z boku izopowierzchni stałej wirowości pierścienia wirowego zaburzonego falą Kelvina o parametrach A/R = 0.04, N = 4 dla czasu: A) t = 0, B) t = 5, C) t = 10
4. PODSUMOWANIE
Przedstawione przykłady obliczeń numerycznych potwierdzają, że w granicy klasycznego, nielepkiego pierścienia wirowego ze skończonym promieniem wewnętrznym rdzenia ewolucja fal Kelvina nie wpływa na jego prędkość przemieszczenia. Rozwój początkowo sinusoidalnego zaburzenia prowadzi z czasem do powstania na obwodzie pierścienia skomplikowanej struktury. Warto odnotować, że fale Kelvina wywołują również ruch cząstek wirowych wzdłuż obwodu pierścienia oraz samorzutne ich grupowanie w obszarach o wysokim gradiencie prędkości. Zjawisko to może powodować niestabilność rozwiązania numerycznego i przerwanie procesu obliczeń. W algorytmie obliczeniowym przeciwdziałano temu efektowi, stosując na etapie redystrybucji metodę „objętości węzła”.
Zaprezentowane wyniki wskazują, że trójwymiarowa metoda cząstek wirowych jest efektywnym i skutecznym narzędziem do modelowania dynamiki struktur wirowych nawet wtedy, jeżeli do obliczeń wprowadzonych zostaje kilkanaście tysięcy cząstek wirowych.
Aproksymacja ciągłego pola wirowości dyskretnym rozkładem cząstek, niosących informację o trzech składowych pola wirowości, wymaga w trakcie procesu obliczeniowego stałej kontroli bezźródłowości pól: prędkości i wirowości. Wartości te w trakcie symulacji numerycznej zmieniały się w zakresie od 10-4 do 10-2 na końcu obliczeń. W trakcie obliczeń numerycznych monitorowane były również niezmienniki ruchu dla cieczy nielepkiej: energia kinetyczna oraz enstropia.
Autor dziękuje prof. Carlo F. Barenghi za konstruktywne i intrygujące pytania, które dały impuls do realizacji niniejszej pracy.
LITERATURA
1. Barenghi C.F., Hanninen R., Tsubota M.: Anomalous translational velocity of vortex ring with finite-amplitude Kelvin waves. “Phys. Rev. E“ 2006, 74, 4, 046303(5).
2. Cottet G.-H., P. Koumoutsakos P.: Vortex methods: theory and practice. New York: Cambridge University Press, 2000.
3. Donnelly R.J.: Quantized vortices in Helium II. Cambridge: Cambridge University Press, 1991.
4. Kiknadze L., Mamaladze Yu.: The waves on the vortex ring in HeII.” J. Low Temp. Phys.” 2002, 126, 1-2, p. 321-326.
5. Kudela K., Regucki P.: The vortex-in-cell method for the study of three-dimensional vortex structures. Tubes, sheets and singularities in fluid dynamics. Series: Fluid Mechanics and Its Applications. Dordrecht: Kluwer Academic Publisher 2002, 71, p. 49-54.
6. Lim T.T., Nickels T.B.: Vortex rings. In: Fluid vortices. Series: Fluid Mechanics and Its Applications. Dordrecht: Kluwer Academic Publ., 1995.
7. Marshall J.S.: Inviscid incompressible flow. New York: John Wiley & Sons, Inc., 2001.
8. Maxworthy T.: Some experimental studies of vortex rings. “J. Fluid Mech.” 1977, 81, 3, p. 465- 495.
9. Quartapelle L.: Numerical solution of the incompressible Navier-Stokes equations. Birkhauser Verlag, 1993.
10. Regucki P.: Modelling of three dimensional flows by vortex methods. Rozprawa doktorska.
Wrocław: Pol. Wrocł., 2003.
11. Saffman P.G.: Vortex dynamics. Cambridge: Cambridge University Press, 1992.
INVESTIGATION OF DYNAMICS OF VORTEX RING WITH FINITE-AMPLITUDE KELVIN WAVES
USING 3D ViC METHOD
Summary. One of the most fascinating and still investigated phenomenon relates to evolution of Kelvin waves (sinusoidal distortions) which appear on circumference of vortex ring. Numerical results reported by Kiknadze and Mamaladze [4] and Barenghi et al. [1] showed that a translational velocity of perturbed ring in superfluid depends on frequency and amplitude of the Kelvin waves. Their numerical simulations refer to the case of inviscid vortex ring with infinitesimal thickness of a core which is realistic for superfluid conditions (zero viscosity and microscopic vortex core thickness). It was interesting to verify if the same behaviour could be observed in the limit of classical inviscid fluid. Numerical results presented in the paper refer to a classical inviscid vortex ring with finite radius of inner core. The dynamics of the ring is modelled using three-dimensional vorticity particle-in-cell method [2].
In spite of the complicated dynamics of rings, it seems that its translational velocities are the same regardless of perturbed or unperturbed structures. In the presented simulations, the Kelvin waves develop along the circumference of the ring and lead to the complicated evolution of vorticity field but do not have the significant influence on the translational velocity.
