ZESZYTY NATJKOWE POLITECHNIKI ŚLĄSKIEJ Seria: ENERGETYKA z. 87
m i
Nr kol. 806
Henryk KUDELA
Instytut Techniki Cieplnej i Mechaniki Płynów Politechniki Wrocławskiej
BADANIA NAD METOD* DYSKRETNYCH WARSTW WIROWYCH
S t r o s z c z e n i e : V pracy zbadano d oś wi ad c za ln ie wpływ wartości pa
r a m e t r ó w n um er y c z n y c h w y s t ę p u j ą c y c h w me to dz i e w ar st w wirowych na dokładność wyników. Otrzymane rozwiązania porównano ze zna
ny mi ro związaniami dokładnymi. Zbadano w p ły w zjawisko, oderwania się warstwy przyściennej na jakość otrzymyw an yc h wyników.
1 • W stęp
V pracy [3] C h o r i n p o d a ł nu me ry c z n ą metodę, zwaną m e todą kropel wiro
w y c h / pa t rz [7] / y p r ze zn ac z on ą do mo d elowania p r z ep ły wó w z dużymi liczba
mi Rey no ld s a / m a ły mi wa rtościami -współczynniki lepkości/. Istotną własnoś cią m e t od y jest br a k w oblic ze n ia ch siatki numerycznej, '..'prowadza ona zw ykle do obliczeń. tzw, lepkość numeryczną, która moż e być togo samego rzędu co lepkość cieczy. M e to da k r o p e l w i r ow y ch posiada jednak pewne wa
dy, do k t ó r y c h n al e żą długie czasy o bl ic zeniowe /liczba działań w jednym k r o k u czasowym jest p r o p o rc j on al na do k w a dr at u liczby w i r ó w / i powolna zbieżność roz wi ąz a ń w pobl i żu ściany • Próbą przezwy ciężenia ty cii wad jest metoda w a r s t w w i r o w y c h również za proponowana przez Chorina [/ł] . Odbywa się to jednak k o s z t e m u p ro sz cz e ń równ a ń ruchu cieczy. Metoda w a r s t w wi ro w y c h służy do rozwiąz y wa ni a równań wars t wy przyściennej Prantla. Ce le m obecnej pracy by ło do św i ad cz al n e zbadanie wpły wu wartoś
ci p a r a m e t r ó w n u m e r y c z n y c h w y s t ę p u j ą c y c h w met od z ie warstw* wir ow yc h na d okładność o t r z ym yw an y ch wyników, a pona dt o zbadanie jak zjawisko oder
wania w a r s t w y prz y śc ie nn e j wp ływa na otrzymywane wyniki.
2. Krótki opis i wy pr owadzenie letody w a r s t w wi rowych Vi1
R ó w n a n i a wars tw y przyściennej v formie równań transportu wirowości w r a z z r ó w n an ie m ciągłości mają postać:
a t § + ( i i - v) E, = v a y y § (
1
)gdzie
3^u + 8y v = 0, (2)
5 = - 3 y U (3 )
a 3(.) oznacza różnicz ko wa ni e po zmiennej (.) , V - opera to: gradientu, V - k in e ma t y c z n y ws pó łc z y n n i k lepkości, u =(u,v) - wektor prędkości, kto-
286 H. Kudela
ro«;o składowa u jest prędkością styc zn ą do ściany, a v składową normalną, L ' x jest skierowana stycznie do ściany, a oś y p ro st o p a d l e do niej. Zak
łada s j ^ , ¿o ściana znajduje się na dodatniej pó łosi x /y=0/, a ciecz w ypełnia całą półpłaszczyznę y > 0 . W a r u n k i g r a n ic z ne są następujące:
u = 0 dla y = 0 , x > 0 (4)
u = ( U . , 0) dla y = oo (5)
W metodzie wars tw w i r o w y c h w yk o rz y s t u j e się tzw. algo r yt m r o zs zc z ep ie nia, zgodnie z k t ó r y m równanie g łó wn o rozkłada się na sumę d wó ch równań elementarnych. R o zw ią za n ie ot rz y mu je się przez k o l e jn o rozwiązy wa n ie rów
nań elcmentarnych. n a s z y m pr zy pa dk u równanie (1) rozkłada się na rów
nanie Eulera / V =0/ oraz równanie dyfuzji. R ó w n a n i e E u le r a rozwiązuje się w z m ie nnych Eagrange*a, w pr ow a d z a j ą c do o b l i cz eń elementy wirow e zwane dalej warstwami wirowymi. R ó w n a ni e dyfuzji n a to mi as t rozwiązuj©
się metodą przy p ad ko we g o błądzenia, dodając do składowej deterministycz
nej położenia warst w y wirowej, losową o rozkła d zi e Gaussa z odpowiednio dobranymi parametrami. R ównania E u le ra w zm ie n ny ch L a g r a n g e ’a opisujące p rz em ieszczenie infinitezy ma ln e j cząstki cieczy mają postać:
dx
dt = H <61
Jak wiadomo, powy żs z e r ównanie opisuje również ruch w a r s t w wirowych. Je
żeli do zaaproksymowania równań (ó) użyjemy ilorazu różni c ow eg o w przód, to algorytm p r ze mi eszczania w a r s t w w i r o w y c h już z u w z g l ę d ni en ie m skład
nika losowego symulującego lepkość ma postać:
x . n+1 = x. + A t • u (?)
x x v 7
n + 1 n _
y i = y i + *V + ? •
gdzie = x. (n-ńt) , x^ = ( x ^ i y , At - k r o k czasowy, ^ - zmienna lo
sowa o rozkładzie Gaussa z war to śc i ą średnią E [ "?]= 0 i w a r i a nc ja var[ ę] = 2-\>At. Prawą stronę równań (6) można wyrazić poprzez wrirowość (3), p osługując się wzor am i (3 } i (2) następująco:
00
u(x,y) = U«, (x) -
J
§ (x,y) dy (8)y
yf
V (x,y) = - a y | u(x,y) dy ¡9)
V m etodzie w a rs tw w i r o w y c h wir ow oś ć zastępowana jest z b io re m nie obraca
jących się równoległych do osi x odc i nk ów / w a r s t w wir ow y ch / o długości h i intensywności 5 . równej różnicy p r ę d ko śc i n a d w a rs t wą i p o d warst wą
S j = u nad " Upod- T a k w i ? c:
S
(x,y)= Z g j S
(y_ yj) IjU),
(10)gdzie 5 (y) jest deltą Diraca, a “X funkc ją c h ar a kt er ys t yc zn ą równą jedności, gdy x
£.[x^-li/2
,x^+h/2
1 , a poza tym równą zero.Śred n ią p r ę d k o ś ć ( z jaką porusza się i-ta warstwa w kie ru n ku osi X, można ob liczyć w st a wi aj ąc w y r a że ni e (10) do w z o r u (8):
,xifh/2 -
Uj = U „ l x ) - -i- (
Z
l>j ^ l y - y j )X
i l X ) dy ) d xZX:-h/2 Ju- i
B a d a n i a n a d m e t o d ą d y s k r e t n y c h w a r s t w w i r o w y c h 2 g 7
y\ >
-U.(x) - i / 2 Bi - S g j d j . m i
gdzie sumowanie odbywa się po wszystk ic h warstwach, dla k t ór yc h y.^>y^ j oraz O ^ d . ^ 1 , d .= 1 - |x. -x . | /h. Po jawienie się składnika 1/2 Jf. można wytłu-
— «3 I rJ \ f \
raczyć n a s t ę p u j ą c o : w ia do m o, żo J^ó(t)=1, gdy t e(t1 , t^J , na tomias t , gdy t = t 1 a rb it r al ni e można przyjąć, że całka ta jest równa 1/2.Składową aprok- symuje się ilo ra ze m centr al n ym
( ą - O y Ą jdzie I
(1 2)
+ są pcczybliżonyrai wartościami J u (x 1 h / 2
)
dy OI+ = ^ fxi 1 h/2J - £ + 5 j V j * 1
d j = 1 "lx i 1 h/z ' x j 1 / h
y * = min (yi)yj),
a sumowanie Z!+ ( o d p . 2 _ ) odbywa się po w s z y s t k i c h w a r s t w a c h takich, że d j+ ^ 1 (d j_ ^ 1 ^ * W a r u n k i b r z e g o w e u(x,y =00) = Uot>(x) oraz v(x,y = 0 )= 0 spełnione są automatycznie. Dla spełnienia w a ru n ku u-s = 0 , gdzie s jest w e k t o r e m jedn os tk o wy m s t y cz ny m do ściany, na ścianie odbywa się
■proces generacji w a r s t w wirowych. Zakłada się, że pole pr ędkości jest antysymetry czne u (x,-y)= - u(x,y) / s t ą d g (x,-y) =_§(* ,y)/. Jeżeli
u (x >^)~u 0 / 0» w z dł uż ściany tworzy się war st wa wirowa o jednostkowej i ntensywności równej 2 u q , k t ó r ą dziel i się na odcinki o długości h i poz
wala się tym o dcinkom dyfund ow a ć w obszar przepływu. Dla lepszego odtwo- rżenia p r oc e su dyfuzji pr og ra m o w o zapewnia si(4 , aby każda generowana w a r stwa miała intensywność nie więk s zą od z góry zadanej liczby . y
^niax miejsce pojed yn cz ej w a rs t wy tworzy się więc pewną /patr zy s tą / liczbę w a r s t w 2 - 1 tak, aby ¿ 5 % 2 u q , 2 1 = | > u „ ajj + 1 , gdzie [ ] oznacza
i=1 J część całkowitą wyrażenia.
•./a r un ek J" (x, -yJ =^(x,y) uwz gl ęd n ia ny jest poprzez zwracanie do cieczy w ar st w wirowych, k tó re prze k ro cz ył y ścianę i przeszły do dolnej półpłasz- czyzny y < 0 . Dla zmniejszenia w a r i a n c j i otrzymy w an yc h wyni kó w C ho r in zastosował procedurę "koncówkowania" /szczegóły patrz [^3, [ 0 / * 3« Bad a ni a n umeryczne
3*1« Uwagi wstępne
Do obliczeń wykorzystano, po n ie wi e l k i c h zmianach, p r o gr am Chorina zadokumentowany pr z ez Cheera [ 0 • O bliczenia były p rowadzono koni ec z
ności dla skończonego odcinka płytki, tzn. dla 0 ^ x ^ a . Aby uniknąć
283 B. Kudela p rz eciwnych do k i e r u n k u prz ep ły w u ni e po ż ą d a n y c h odd z ia ły wan C h o r i n zap- proponował n as tępującą procedurę*: warstwy wirowe, k tó re przekroczyły x = a f są z obliczeń eliminowane, a* warstwy, k t ó r y c h środki znalazły się w pasku [a-2h, aj,dozna ją ruchu k o n w e kc y jn eg o tylko w k i e r u n k u osi x. W obocnej pracy przyjęto a=1,4. Jałto pi er w s z e rozw ią zy w an o zagadnienie formowania się warstwy p rz yś ciennej nad płask ą płytką, gdy zewnętrzne pole prędkości jest stałe (x ) = U Q • S ta cj o n a r n y m r o z w ią za ni e m takie
go zagadnienia jest dla (1) -(5) dobrz e znany profil Bla si us a [9 , s . 1 3 4 ], N as tę pn i e rozwiązywano' zagadnienie, w k t ó r y m zewnętrzna prędkość zmie
niała się jak Utso(x)= U Q - bx , U ^ O , b > 0 . Z a gadnienie to by ło roz
wiązywano przez llowarthfa [ó] m e to dą rozwinięć asymptotycznych, Do obec
nej pracy rozwiązania, k t ór e bę dziemy uważali za dokładne, zaczerpnięto 2 [7 s. 598] . Dozwiązania tg zależą od dwóch bezwymiarowych zmiennych x * = - ~ oraz 2 = 1/2- y — • W e d ł u g H o w a r t h ’a dla x * = 0,12 za- chodzi ooderwanie warstwy przyściennej. Do o b l i c z e ń p r z y ję t o najpierw rozkład 1- 0,0 7 5 x, dla k t ó re g o pu nk t oderwania p rzypada dla x-1,6 , czyli poza rozp at ry w an ym obszarem, a n a s t ęp ni e U ^ x ) = 1 -0 ,2 3 x, gdzie punkt oderwania przypada dla x = 0 ,5 1 , tzn. wewnątrz rozpatrywane
go obszaru. Celem zmniejszenia b ł ę d u statyst yc z ne go ro zk ła dy prędkości uśredniano po 20 k r o k a c h czasowych. J e ż el i i nt en sywność no w o pow st aj ą cej warstwy była mn iejsza od z góry zadanej liczby £ o , to w a rs t wa ta była pomijana. K obecnej pracy p r z yj ę to = ^^-^niax'
3.2. Wyniki numeryczne
Parametrami, który mi można w p ł yw a ć na p r oces obliczeniowy, są: § a a x t h i A t . uf pracy i 5l Chor in podał, że dla uniknię ci a oscylacji rozwiązań należy przyjmować A t i h, D laej p r z yj ęt o w e w s z y s t k i c h eksperym e nt ac h
A t = h i oznaczono jedną literą r /At=h=r/. D o bó r Ę oraz r jest max
k o m p ro mi so m pomiędzy dokładnością a cz asom trwania obliczeń. Przykładowo^
dla zagadnienia, gdy U OQ(x) = 1 , £ max = 0 , 1 i r=0, 1 , li czba w a r s t w wiro- w yc h po 6 0 krokach czasowych wynosiła 417, a czas trwania obliczeń pola prędkości dla x = 0 ,2 do 1, co 0,1 na mas z yn ie OD RA -1 30 6 w y n o si ł ok. 100 minut, natomiast dla tego samego zagadnienia z ¿fraax = 0,8 i r=0,2 liczba w i r ó w wyno s ił a 3*1, a czas o b l i c ze ń ok. 2,5 minuty.
W pracy [4] eksperymentalnie stwierdzono, że dla A t ¿'0,2 , h < 0 , 2 oraz 5 raax* 0 , 1 domin uj ąc y m b łę d em jest b ł ą d statystyczny, który w r a z ze w zr os t e m liczby w i ró w m al ej e raczej wolno. Na r y s . 1 pr z ed s t a w i o n o roz- kłady p r ę dk oś c i dla Uoc(x)= 1, a na rys. 2, gdy U ^ x ) = 1-0,075x. Na tle rozwiązań dokład n yc h nanie s io no r oz wiązania n u m e r y c z n e : znak A odnosi się d o j§mQX = 0,1, na t omiast * do J n ax = R y s u n k o m a/ odpowiada r = 0 ,1, a rysunkom b/ r=0,2. P o d rysunkami p o da n o równi eż liczbę N warstw wirowych, k tóre brały udział w obliczeniach.
B a d a n i a n a d m e t o d ą d y s k r e t n y c h w a r s t w w i r o w v c h ^
Iiys. 1. Uśrednione rozkłady prędkości dla U x=1, x=0,5, v=10"4 a/ r=,1 , û - ^ max= ° . 1. N=^17, N=85 b/ r=0,2, * - £ „ „ = 0 . 1 , Ne 12, - - £ „ „ = 0 , 8 , N=34
Rys. 2. Uśrednione rozkłady prędkości dla U x=1,075x, x=1, ^ = 10^
• a/ r=0<1 < ‘*§.nax=0-»1 - N=i*'17- ‘- ^ a x = ° > 8 > N = 121 b/ r = 0 ,2, * - f m „ = 0 , 1 , N= 162, « - ^ „ = 0 , 8 , N=<ł7
Przedstawione wyniki otrzymano po 60 krokach czasowych. Przyjęto, że w przybliżeniu jest osiągany wtedy stan stacjonarny. Częściowo potwier
dziły to obserwacje pól prędkości po 40,Ó0 i 80 krokach czasowych. Naj
lepsze rozwiązania otrzymano dla r = 0 ,2 i ^ =0,1. Można zauważyć, że ro
D max 1
wiązania dla r=0, 1 , .5lnax=°j ® wykazują pewne oscylacjo. Zwiększenie r / r= At=h/ powodowało tłumienie tych oscylacji. Rys. 3 'Przedstawia przykładowe usytuowanie warstw wirowych nad płytką po
60
krokach czasowych dla Uoo(x)= 1 oraz r=0,1 i P = 0 , 1 . max 7
Aby zbadać wpływ oderwania warstwy przyściennej na otrzymywane wyniki przeprowadzono obliczenia dla U ^ x ) = 1-0,23x. Nie uzyskano przepływu stacjonarnego, tfarstwa przyścienna gwałtownie się poszerza. Pojawia się przepływ wsteczny dla x ^ 0 , 5 1 * Rozkłady prędkości, mimo uśredniania, są bardzo nieregularne^ a ich ważność jest wątpliwa. Przykładowy rozkład wirów nad płytką dla t^fxjs
1
-0
,23
x, r=0, 1 i £„,„„=0,5 przedstawiono na290 B . Kulisia
0,08-r
ow-
• • • • ••• • • ©•o • • • e o o*
. • . »>■*»»«
o ...
R y s .
3
. W ar st wy wirowe nad płytką: U^,(x)=1, r=0,1 » ^ n a s 2 0 '1 ’ N=417 o- warstua wirowa o intensywności podatniej, • - warstwa wirowa o intensywności
0 intensywności p o d a t n i e j, •- we 1 ujemnej^ ^=10
0,08r
ąot.
■ •« o* *• o o* o o • • • e 0% m o #•*•« O» O O O o • O o • ooł,o»*• o »o*
• *o * • • ••© o •• o O Q OO » <0 »O łOO««Q o«» »O««
o 1
Rys. 4. Warstwy wirowe nad płytką: U*.(a^= 1 -0,2 3 x , r = 0 ,1 , 5 n ax=^ ^ f N=307 o- warstwa o intensywności dodatniej, •- warstwa o intensyimości ujemnejj V*10
4. »'nioski i uwagi końcowe
Na podstawie przedstawionych wyników można wnioskować, że:
t/ przy ograniczonej liczbie kroków czasowych na ogół nie otrzymuje się najlepszych dokładności przy najmniejszych wartościach At , h , 2/ duże wartości ^ ^ ^ = 0 , 8 przy stosunkowo małych wartościach r / r=dt=h=0,y
mogą powodować, nie mające fizycznego znaczenia, oscylacje rozwiązać, 3/ obecna wersja metody warstw wirowych nie nadaje się do modelowania
przepływów z oderwaniem warstwy przyściennej.
Należy zwrócić uwagę, że metoda warstw wirowych /podobnie jak metoda kro
pel wirowych /nie jest zorientowana na otrzymywanie bardzo dokładnych rozwiązań, lecz nacisk położony jest w niej na symulację fizycznej strony przepływu cieczy/ wprowadzenie do obliczeń elementu wirowego, generacja wirów na ścianie/.
Równania Prandtla zostały wyprowadzone na podstawie przyjęcia pewnych hi
potetycznych założeń, z których wynika między innymi v ^ u . C h o r i n w pry
watnej korespondencji zasugerował, że być może należałoby sprawdzać, czy ten warunek jest spełniony dla każdej warstwy. Jeżeli n i e , to warstwę wi
B a d a n i a n a d m e t o d ą d y s k r e t p y o h w a r s t w w i r o w y ob 291
rową należałoby zastąpić kroplą wirową /por. [
2
], [5
j/. V ten sposób zastępowałoby "dopasowywanie*równań PranĄtla do równań Naviera-Stokesa i być może pozwoliłoby to na modelowanie przepływów z oderwaniem warstwy przyściennej. Prace w tym kierunku zostały już pdojętd przez autora.Wyniki zostaną zaprezentowane w terminie późniejszym.
Literatura
Lii Cheer A.Y.: BOUNDL: A Program For Calculating Flow Past a Semi- Infinite Flat Plate Using the Vortex Sheet Method, Lawrence Berkeley
Laboratory University of California, Physics Computer Science and Math. Division, IBL - 66kj Supplement, 1978.
[2] Cheer A.Y.J A Study of Incompressible 2-D Vortex Flov Past a Circular Cylinder, ibid, LBL-9950,1979
[3] Chorin A.J.; Numerical Study of Slightly Viscous Flow, J.Fluid Mech.
¿ 2 s - 785-796, 1973
[4] Chorin A.J.j Vortex Sheet Approximation of Bottdary Layers, J.Comp.
Phys. 22 s * *»28-*»*»2, 1978
[5] Chorin A.J. ; Vortex Models and Boundary Layer Instability, SIAM J.Sci. Stat. i, s. 1-21, 1980
[
6
] Howarth L.; On the solution of the Laminar Boundary Layer Equations, Proc. Roy. Soc. London A919, 169, s. 5*17-579, 1938[
7
] Kudela H.J Metody dyskretnych wirów, w aMetody Analityczne i Kune- ryczne w Mechanice Płynów - wykłady - s. 111-130, Szkoła Letnia Mechaniki Płynów, Mikołajki 198318] Kocin N.E., Kibel I.A., Rozę N.M.ZTeoreticeskaja gidromechanika, OGiZ Moskwa
1963
, t. UL9] ęchłichting H.J Teoria progranl&iowo słoja, Nauka, Moskwa 197*».
HCCJIĘHOBAHKE MESOM BH2P0BUX CliOEB
P e s
x u
e£ paOote HsyveBo skcnepKKSHTaaLKO B jm z s z e pa3Mepa vacxeHBHx napaMeipos weToxa BExpoEHX cjtoBb Ha ToiHoctt noayvaessHx pesynbiaioB. Pe3yxbiaiu cpas~
H e H O c tohhmmh pemeHHSMB. H c c x e x o B a s o B i M B u e io v k h cenapaRHB norpaHHtsoro exon sa noxyvaeHue pe3yxs,tastu.
292
S TUDY OF T H E V O R T E X SHEET M E T H O D
S u m m a r y
The paper experimental? investigates the influence of the numerioal parameters in the vort ex sheet m ethod on the accuracy of the results.
The results were compared with the known aoorate solutions. The influence of the boundary separation point on the results w hi ch were being obtained was also investigated.