1. No category

UNIVERZA V LJUBLJANI
FAKULTETA ZA NARAVOSLOVJE IN TEHNOLOGIJO
ODDELEK ZA FIZIKO
Aleš Berkopec
PROBLEM PROSTE POVRŠINE
V TOKU VISKOZNE TEKOČINE
DIPLOMSKA NALOGA
Mentor:
prof. dr. Alojz Kodre
Ljubljana, februar 1994
POVZETEK
Diplomska naloga je posvečena problemu nestisljive viskozne tekočine s prostim
robom. Pri stacionarnih problemih te vrste je poleg hitrostnega in tlačnega polja
del rešitve tudi oblika prostega roba.
V prvem poglavju smo predstavili spremenljivke, s katerimi opišemo tekočino.
Ob upoštevanju nestisljivosti in konstantne viskoznosti lahko problem skrčimo na
reševanje sistema dveh enačb. V drugem poglavju smo izbrali konkreten problem
ter enačbe oklestili do brezdimenzijske oblike in definirali robne pogoje. Stanje
raziskav na področju numeričnega reševanja problemov s prostim robom je opisano
v uvodnem delu tretjega poglavja, ki je v celoti posvečeno reševanju glavnega primera z numeričnimi algoritmi. Izbrali smo tri načine reševanja: dve sorodni celični
metodi in metodo končnih elementov. Po podrobni obdelavi smo z vsako od metod
rešili preizkusni primer. Metoda sledilnih delcev je kot najstarejša predstavljena
prva, naslednja pa njej sorodna metoda robnih sledilnih delcev. Popolnoma nov
način reševanja temelji na metodi končnih elementov.
Po primerjavi v četrtem poglavju smo izbrali najustreznejšo metodo za pregled prostora parametrov, ki mu je namenjeno peto poglavje. V njem smo razdelili prostor
parametrov na dva dela glede na obstoj stacionarne oblike. Ker nam je na voljo le
del prostora, smo na koncu predzadnjega poglavja pokazali, kje lahko pričakujemo
stacionarne oblike in označili nekaj zanimivih primerov.
V zaključku so podani sklepi in nakazana smer raziskav na področju izpopolnitve
predstavljenih algoritmov.
i
ABSTRACT
The work discusses the incompressible viscous free boundary flow. The steady flow
solution consists not only of velocity and pressure field. The position of the free
boundary is a part of the solution as well.
In Chapter One we introduce the quantities that govern the flow. After the incompressibility and the constant viscosity postulate are applied, the problem can be
reduced to two basic equations. The main problem is presented in Chapter Two,
furthermore the nondimensional form of the equations and boundary conditions
are derived. We begin Chapter Three with general overview of the solution methods for incompressible viscous free surface flow. Three algorithms are described in
detail: benchmark problems and the main problem are solved with each of them.
The Marker-and-Cell (MAC) method, as oldest, is considered first. From MAC
the Surface-Marker (SM) method was derived which is presented next. Completely
new approach is introduced with the Finite Element Method.
After comparison of methods in Chapter Four we choose the most appropriate
method for exploration of the parameter space. Chapter Five takes us to the
examination of the parameter space in order to find the parameter values for which
the steady solution of the main problem can be reached. The parameter values are
limited, therefore only part of the parameter space is actually available. A sketch
at the end of the chapter gives a hint about conditions that lead to the steady free
surface flow.
The final Chapter summarises the work and presents some ideas for future research
in the area of numerical methods.
ii
VSEBINA
POVZETEK
i
ABSTRACT
ii
VSEBINA
iii
1 UVOD
1
2 FORMULACIJA PROBLEMA
2
2.1 Osnovne enačbe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2 Brezdimenzijske količine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.3 Robni pogoji . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5
3 NUMERIČNE METODE
7
3.1 Metoda sledilnih delcev (MAC) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.1.1 Opis metode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.1.2 Preizkusni primer: odtekanje pravokotnega bloka tekočine . . . . . . 11
3.1.3 Rešitve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Metoda robnih sledilnih delcev (SM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2.1 Opis metode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2.2 Preizkusni primer: odtekanje pravokotnega bloka tekočine . . . . . . 15
3.2.3 Rešitve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.3 Metoda končnih elementov (FEM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3.1 Opis metode - Galerkinov pristop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
iii
VSEBINA
iv
3.3.2 Povezava tlaka s hitrostnim poljem (penalty method) . . . . . . . . . . 21
3.3.3 Matrični členi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.3.4 Preizkusni primer 1: dvodimenzionalni vrtinec v votlini . . . . . . . . . 28
3.3.5 Preizkusni primer 2: odtekanje eliptičnega bloka tekočine . . . . . . . 29
3.3.6 Rešitve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4 PRIMERJAVA METOD
32
4.1 Prostorska zahtevnost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.2 Časovna zahtevnost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.3 Ostale prednosti in slabosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5 RAZDELITEV PROSTORA PARAMETROV . . . . . . . . . . . . . . . . . . . . 35
5.1 Posamezni zgledi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5.2 Razdelitev prostora . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6 ZAKLJUČEK
39
ZAHVALA
40
LITERATURA
41
1
UVOD
Gibanje tekočine določimo s hitrostnim poljem v(x, y, z, t) in dvema izmed termodinamskih spremenljivk, recimo tlakom p(x, y, z, t) in gostoto ρ(x, y, z, t). Tlak v
mirujoči tekočini nastane zaradi lastne teže, imenujemo ga hidrostatski tlak. Kadar
se tekočina premika, vplivajo na tlačno polje tudi medsebojne sile delcev in togih
sten. Spreminjanje gostote je odvisno od stisljivosti in temperaturne razteznosti
tekočine. Obravnavali bomo nestisljive tekočine s konstantno viskoznostjo (newtonske). Posebnost naloge je modeliranje prostega roba. Prosti rob je tisti del
tekočine v polju sile teže, ki je v stiku z zunanjim medijem, kamor ne sodijo toge
stene. Lega prostega roba v splošnem ni znana in predstavlja del rešitve problema. Tu veljajo posebni robni pogoji, katerih posledica je popolnoma druga narava
gibanja, saj se delci na prostem robu lahko premikajo v poljubne smeri, tekočina ob
togi neprepustni steni pa ne. Analitični prijemi so za te vrste problemov primerni
le za majhno množico posebnih primerov, a kot omenjata J.M. Floryan in H. Rasmussen v svojem preglednem članku [8], do danes še ni razvita splošna numerična
metoda za reševanje dinamike tekočine s prostim robom. Obstaja vrsta algoritmov, ki prav tako kot analitične rešitve dajejo zadovoljive rezultate le za posebne
primere. Tako smo si tudi v diplomski nalogi kot problem zastavili reševanje toka
viskozne nestisljive tekočine v izbrani geometriji, kjer je bil eden od ciljev razvoj
nove metode za iskanje stacionarnih stanj s pomočjo trikotniške mreže končnih
elementov.
Omejili se bomo na dvodimenzionalni primer - v vogalu se vrti neskončna prizma
tekočine, opazujemo obliko njenega prečnega prereza. Pri računanju ne bomo upoštevali sil površinske napetosti. Njihov vpliv lahko zanemarimo, če je kaplja dovolj
velika. Glavni cilj naloge je raziskati spekter parametrov, ki določajo tok tekočine.
Prostor parametrov želimo razdeliti v dva dela glede na obstoj stacionarne oblike.
Hkrati bomo preizkusili nov algoritem za reševanje danega problema s prostim
robom, ki temelji na relaksaciji trikotne mreže končnih elementov. Primerjali ga
bomo z že znanimi in raziskali njegove prednosti in slabosti.
1
2
FORMULACIJA PROBLEMA
V pravokotnem vogalu, ki ga tvorita spodnja in leva stena, obe neskončno dolgi,
je velika kaplja nestisljive newtonske viskozne tekočine. Zanimajo nas stacionarna
stanja kaplje v težnem polju, če se spodnja plošča premika s konstantno hitrostjo
u0 proti pokončni plošči (sl. 2.1), težni pospešek in materialne lastnosti kaplje pa
so parametri gibanja.
...
...
...
.
........
.....
B
C
g0
.........................
.................. ......
..................... .....
........................ ..............
......................
............................ ............................................................ ...................
................... ................. ........... ............
................... ................. ........... .............
................... ................. ........... .............
................... ................. ........... .............
................... ................. ........... .............
................... ................. ........... ............
. . . . . . . . . . . . . . . . . . . . . . . . . ....... ...
...............................................
u0
A
Sl. 2.1: Skica obravnavanega problema, ki nam je služil kot primer pri modeliranju
nestisljive newtonske viskozne tekočine s prostim robom. Iskali bomo
stacionarno obliko kaplje pri različnih parametrih toka. Robova BC in
AC sta toga: prvi miruje, drugi se giblje v nakazani smeri s konstantno
hitrostjo u0 . Oblika prostega roba AB predstavlja del rešitve problema.
2.1
Osnovne enačbe
Predpostavimo, da je temperatura tekočine konstantna in povsod ista. Ker je
tekočina nestisljiva, je tudi gostota ρ konstantna, zato velja kontinuitetna enačba
v obliki:
∇·v = 0 ,
(2.1)
∂
∂
kjer je vektor ∇ = ( ∂x
, ∂y
), vektor v = (u, v) pa hitrost. Iz drugega Newtonovega
zakona sledi ohranitveni zakon za gibalno količino [6,15]
η
∂v
+ (v · ∇)v = f − ∇p + ∇2 v
∂t
ρ
.
(2.2)
To je Navier-Stokesova enačba, kjer so f volumske sile, p tlak in η viskoznost, ki je
po predpostavki konstantna. Čeprav obravnavamo stacionarni problem, časovnega
odvoda ne bomo opustili, ker bomo pri reševanju uporabili tudi psevdočasovne
metode.
2
2 FORMULACIJA PROBLEMA
2.2
3
Brezdimenzijske količine
Sistem enačb poenostavimo z uvedbo karakterističnih vrednosti in brezdimenzijskih spremenljivk. Vsako spremenljivko nadomestimo s produktom karakteristične
vrednosti in brezdimenzijske spremenljivke. Za primer parametrizirajmo hitrost:
v → u0 · v ,
(2.3)
kjer je u0 karakteristična hitrost, v na desni strani pa brezdimenzijska hitrost. Za
karakteristično hitrost u0 za našo nalogo izberemo hitrost spodnje plošče. Ker
pričakujemo, da bo v stacionarnem stanju u0 največja hitrost v sistemu, smo tako
približno omejili brezdimenzijsko polje absolutnih vrednosti hitrosti | v | na vrednosti med 0 in 1.
Podobno naredimo pri tlaku. Največji tlak pričakujemo v točki C (sl. 2.1), kjer
se tekočina, ki se lepi ob gibajočo se spodnjo ploščo, nariva ob levo steno. Hitrost
spodnje plošče je u0 , prav tako pa tudi tekočine, ki jo plošča nosi s seboj. Ob
navpični levi steni se v vogalni točki izniči vodoravna komponenta hitrosti, tako
da se tekočina ustavi in začne lesti ob navpično steno. Hkrati se v točki C ustvari
zastojni tlak pC = 21 ρu20 , ki bo v vogalni točki največji, saj v nobeni drugi točki ne
pride do tako občutnih sprememb hitrosti. Če vzamemo za karakteristično vrednost
tlaka kar zastojni tlak 12 , je tlačno polje tekočine normirano na vrednost 1 (če ne
upoštevamo hidrostatskega tlaka). Pisanje si poenostavimo, če tlake normiramo
na dvojno vrednost zastojnega tlaka v vogalni točki, P0 = ρu20 . Tako so tlaki v
tekočini navzgor približno omejeni na brezdimenzijsko vrednost 12
p → P0 p
.
(2.4)
Na koncu parametrizirajmo še čas in dolžino:
t → T0 t ,
(2.5)
(dx, dy) → L0 (dx, dy) .
(2.6)
T0 je povezan s karakterističnima hitrostjo in dolžino kot T0 = L0 /u0 . Enačba
nam pove, da je karakteristični čas tisti, v katerem delec tekočine s hitrostjo u0
prepotuje razdaljo L0 . Smiselna vrednost karakteristične dolžine
√ L0 je kvadratni
koren ploskve, ki jo zavzema dvodimenzionalna tekočina L0 = S.
Tako definirane spremenljivke vstavimo v Navier-Stokesov sistem (2.2), upoštevajoč
silo teže kot edino prostorninsko silo. V sistemu enačb brezdimenzijskih količin
∂v
1 2
0
+ (v · ∇)v =
∇ v
(2.7)
− ∇p +
1
− F r2
∂t
Re
sta ostala le še dva parametra. Parameter, v katerem nastopa težni pospešek g0 ,
smo poimenovali F r. Imenujemo ga Froudovo število F r in podaja razmerje med
2 FORMULACIJA PROBLEMA
4
kinetično in potencialno energijo ali kvadrat razmerja med hitrostjo spodnje plošče
in hitrostjo, ki jo doseže tekočina pri prostem padu z višine L0
1
g0 L0
= 2 .
2
Fr
u0
(2.8)
Drugi parameter, ki določa dinamiko kaplje, je Reynoldsovo število Re. Povezuje
karakteristično hitrost u0 in dolžino L0 z gostoto ρ in viskoznostjo sredstva η. Re
podaja razmerje med inercialnimi silami in silami zaradi notranjega trenja
Re =
ρu0 L0
.
η
(2.9)
Pri majhnih Reynoldsovih številih je tok tekočin urejen, saj v enačbi (2.7) prevladuje člen z Laplaceovim operatorjem ∇2 v, katerega rešitve so pohlevne in stabilne. Táko gibanje tekočine imenujemo laminarno. Njegovo nasprotje je turbulentno gibanje, ki ga lahko opazimo pri velikih Re, ko je tok neurejen in kaotičen,
ker ima v Navier-Stokesovi enačbi najmočnejši vpliv nelinearni člen (v · ∇)v. V
Reynoldsovem številu sta skriti tudi viskoznost η in gostota ρ, edini materialni
konstanti snovi, ki ju potrebujemo v izračunu. Njun količnik ima ustaljeno ime
- kinematična viskoznost ν = ηρ . V tabeli 2.1 navajamo vrednosti viskoznosti,
gostote in kinematične viskoznosti za nekaj različnih snovi.
snov
kg
ρ [m
3]
kg
η [ ms
]
voda
zrak
etilni alkohol
glicerin
živo srebro
998
1.29
790
1260
13550
10−3
1.8 · 10−3
1.2 · 10−3
0.85
1.56 · 10−3
ν=
η
ρ
2
[ ms ]
10−6
1.4 · 10−5
1.5 · 10−6
6.7 · 10−4
1.15 · 10−7
Tab. 2.1: Gostota, viskoznost in kinematična viskoznost za različne snovi [14,15].
Po poenostavitvi enačb, ki opisujejo gibanje tekočine v gravitacijskem polju, vidimo, da so stacionarne oblike kaplje odvisne le od dveh parametrov. Dvoparametrični prostor (Re, F1r2 ) bomo v zadnjem poglavju pregledali glede na formiranje
stacionarnih oblik. Preden se tega lotimo, poglejmo, kakšni so robni pogoji, in
izberimo učinkovito računsko metodo.
2 FORMULACIJA PROBLEMA
2.3
5
Robni pogoji
Robni pogoji povezujejo spremenljivke na robu definicijskega območja enačbe. Dani
problem ima tri robove: dve togi steni, od katerih se spodnja giblje s konstantno
hitrostjo u0 v tangencialni smeri, in prosti rob. Zahteve, ki naj bodo izpolnjene na
robovih so lahko posledica postavitve problema, kot je v našem primeru s hitrostmi
na togih stenah, lahko pa jih izpeljemo iz gibalnih enačb, kar bomo storili pri
izpeljavi robnih pogojev za tlak.
Če predpostavimo, da sta steni neporozni in nepremični v normalnih smereh, potem
sta normalni hitrosti ob stenah enaki nič. Tekočina naj se na steno lepi, zato ima
tam enako hitrost kot stena sama, ali
vn = 0 ,
vt = vstena
(2.10)
(2.11)
,
kjer sta n normalni, t pa tangencialni vektor na steni, vn in vt pa normalna in
tangencialna hitrost. Tlak na togih stenah izpeljemo iz Navier-Stokesovih enačb.
Sledi robni pogoj za tlak
∂p
1 ∂ 2v
=
∂n
Re ∂n2
.
(2.12)
Na prostem robu je tangencialna komponenta hitrosti neomejena; če naj na prosti
rob ne deluje nobena sila, pa velja
2 ∂u nx −
Re ∂x
2 ∂v p−
ny −
Re ∂y
p−
1
Re
1
Re
∂u ∂v +
ny = 0 ,
∂y
∂x
∂u ∂v +
nx = 0 ,
∂y
∂x
(2.13)
(2.14)
kjer je n = (nx , ny ). V stacionarnem stanju so hitrosti na prostem robu le tangencialne, zato je
v·n = 0 .
(2.15)
Namesto tlačnih robnih pogojev (2.13) in (2.14), bomo uporabili približek za velika
Reynoldsova števila [10]
Re → ∞ ⇒ p
= 0 .
(2.16)
∂Ω
Prosti rob je v stacionarnem stanju tako izobara. Za to aproksimacijo se odločimo
zaradi delikatnega računanja normale. Zanjo potrebujemo lokalno ukrivljenost
roba, s tem pa je, še posebej pri celičnih metodah, precej težav.
2 FORMULACIJA PROBLEMA
6
v·n=0
u=0
v=0
∂p
∂x
=
1 ∂2u
Re ∂x2
.............................
.................... .......
...............................
.............................. ........
.
.
.
.
.
.
.
.
.
.
.
....................................... ................
............................
............................ ................
............. . . . ..........
............................ ............................................................................ .........................
............................ ........................ ................ ........... ......
............................ ........................ ................ ............. .......
............................ ........................ ................ .............. ....
............................ ........................ ................ ..................
............................ ........................ ................ ..................
............................ ........................ ................ ....................
............................ ........................ ................ ..................
............................ ........................ ................ ..................
............................ ........................ ................ ................
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .......... .....
p=0
u = u0
v=0
∂p
∂y
=
1 ∂2v
Re ∂y 2
Sl.2.2: Robni pogoji.
Začetnim pogojem posebnega poglavja ne bomo posvečali. Iščemo stacionarno
stanje, ki je od začetnih pogojev neodvisno, seveda pa bomo računali z začetnimi
oblikami, ki so podobne pričakovanim končnim, ker tako zmanjšamo računski čas.
3
NUMERIČNE METODE
Ne poznamo splošne metode, s katero bi lahko reševali raznolike probleme s prostimi
robovi. Na voljo so metode z omejenim območjem uporabnosti, ki delujejo sicer
dobro, vendar le pri posebnih pogojih. Prosti rob se lahko zvija, prepogiba, širi
ali še kako drugače deformira, zato mu je težko slediti. Tekočini lahko sledimo po
Eulerjevem ali Lagrangeovem načinu [8]. Pri prvem je računska mreža pritrjena
na laboratorijski sistem, v drugem pa vezana na tekočino. Obstajajo tudi metode,
ki jih lahko uvrstimo v oba razreda.
Pri Eulerjevem načinu moramo vnaprej določiti računsko domeno, kar je prikladno
za prostorsko omejena gibanja (tekočina v posodi, ipd.). Prosti rob je določen le
do mrežne razdalje natančno in zaradi te pomankljivosti naletimo na probleme
pri vklapljanju vplivov, ki so odvisni od oblike površine, na primer površinske
napetosti, ki je sorazmerna ukrivljenosti površine. Te metode so precej potratne,
saj tekočina zapolnjuje le del računske domene.
Pri Lagrangeovem načinu opisovanja naletimo na drugačne probleme. Čeprav s
točkami, vezanimi na tekočino, položaj roba natančno poznamo, smo pri divjem
gibanju (zvijanje površine, lokalna rotacija) v težavah, če opazujemo časovni razvoj
oblike, saj se topologija računske mreže lahko zlomi. Metoda je ekonomična, vendar
zaradi omenjenih težav uporabna za reševanje manjšega dela prostorobnih problemov, predvsem stacionarnih in takšnih, kjer vsaj približno lahko vnaprej uganemo
obliko roba.
Oba zahtevnejša postopka (iskanje roba in premikanje mrežnih točk) sta netrivialni nalogi, najustreznejšo metodo reševanja pa lahko izberemo šele, ko poznamo
problem.
Iz želje po čim boljšem postopku za reševanje diplomskega primera, se je rodila ideja
o novi metodi, ki bi tekočino in njen prosti rob simulirala s prilagodljivo trikotno
mrežo končnih elementov. Izračunov nismo mogli preveriti z eksperimentalnimi
rezultati, saj nam ni znano, da bi bili kjerkoli objavljeni. Zato smo, da bi izvedeli,
kako uspešna je nova metoda, primerjali rezultate z drugimi, že uveljavljenimi
metodami, kot sta metodi sledilnih in robnih sledilnih delcev.
7
3 NUMERIČNE METODE
3.1
Metoda sledilnih delcev (MAC)
3.1.1
Opis metode
8
Metoda MAC je bila razvita v šestdesetih letih tega stoletja in je namenjena opazovanju časovnega razvoja nestisljive viskozne tekočine s prostim robom [9,10].
Na prostor, kjer bomo opazovali gibanje tekočine, najprej položimo pravokotno
mrežo, ki domeno računanja razdeli v kvadrate, ki jim rečemo celice (sl. 3.1).
Vrednosti spremenljivk izračunavamo le v diskretnih mrežnih točkah, vrednosti v
vmesnih točkah prostora pa lahko dobimo z interpolacijo. Potem z označitvijo celic
definiramo del prostora, ki ga zapolnjuje tekočina. Če v dani celici ni tekočine, jo
označimo za prazno. Če je v dani celici tekočina, katerakoli od njenih sosed pa je
prazna, jo označimo za robno. V primeru, ko celica s tekočino nima praznih sosed,
je polna. Toge stene definirajo zidne celice.
Sl. 3.1: Razdelitev računske domene na celice.
V polne in robne celice nato enakomerno posujemo sledilne delce ali markerje. To
so brezmasni delci, vezani na tekočino. Vsak marker potuje s tekočino tako, kot mu
narekuje hitrostno polje na mestu, kjer se trenutno nahaja. Ko markerji potujejo
po domeni, se celicam spreminja status. Prej robno celico lahko zapustijo vsi
markerji in tako postane prazna. Postopek spreminjanja statusa celic je sestavljen
iz malega števila korakov, ki jih bomo kasneje opisali. Status celic potrebujemo, da
določimo območje računanja. Na prostem robu je tlak določen z ustreznim robnim
pogojem, zunaj tekočine pa je enak zunanjemu tlaku. Ker tlačno polje računamo
v diskretnih celičnih točkah, si z robnimi in zidnimi celicami omejimo računsko
domeno za tlak.
Stacionarno stanje in rešitev problema poskušamo doseči s psevdočasovnim postopkom. Polni Navier-Stokesov sistem (2.2) rešujemo kot časovni problem. Če se
3 NUMERIČNE METODE
9
oblika kaplje med dvema časovnima korakoma dovolj malo spremeni, smo dosegli
rešitev. Pri tem postopku v vsakem časovnem koraku izračunamo hitrostno polje
iz tlačnega, nato pa tlačno polje popravimo tako, da ima v naslednjem časovnem
koraku hitrostno polje divergenco nič. Nato popravimo hitrostno polje in premaknemo markerje. Po vsakem premiku markerjev celicam še popravimo status.
Z metodo se je uveljavila uporaba premaknjene mreže (staggered grid). To je
posebna mreža, pri kateri računamo različne spremenljivke v različnih mrežnih
točkah [7]. Numerični postopki, s katerimi rešujemo tekočinske probleme, so dostikrat nestabilni. Z uporabo premaknjene mreže pa dosežemo, da so krajevni odvodi
prvega reda aproksimirani z diferencami prvega reda, s čimer smo se znebili parazitskih rešitev, ki razpihnejo vrednosti in naredijo račun nestabilen. Na več načinov
se lahko odločimo, v katerih točkah bomo izračunavali u, v in p. Naš izbor premaknjene mreže je prikazan na sl. 3.2.
•
•
k+1/2
vj+1
k+1/2
vj
• k
•
uj−1/2 pkj
• k
•
uj+1/2 pkj+1
•
• k
uj+3/2
•
k−1/2
vj
k+1/2
vj+1
Sl. 3.2: Premaknjena mreža (staggered grid).
Zapis unj,k je diskretizirana spremenljivka u v točki (j, k) ob času n · ∆t. V taki
shemi se ux v točki (j, k) zapiše kot
∂u
∂x
n
=
unj+1/2,k − unj−1/2,k
∆x
j,k
+ O(∆x2 )
.
(3.1)
Navier-Stokesovi enačbi in enačbo za nestisljivost izračunamo v različnih točkah:
komponento x Navier-Stokesovega sistema v točkah (j + 1/2, k), komponento y v
točkah (j, k + 1/2) in ∇ · v = 0 v točkah (j, k) [7].
Pri diskretizaciji enačb potrebujemo naslednje diferenčne izraze
∂u
∂t
n
j+1/2,k
=
n
un+1
j+1/2,k − uj+1/2,k
∆t
+ O(∆t)
,
(3.2)
3 NUMERIČNE METODE
∂u2
∂x
n
=
j+1/2,k
n
∂(uv)
∂y j+1/2,k
2 n
∂ u
∂x2 j+1/2,k
2 n
∂ u
∂y 2 j+1/2,k
n
∂p
∂x j+1/2,k
10
(u2 )nj+1,k − (u2 )nj,k
+ O(∆x2 )
∆x
,
(3.3)
=
(uv)nj+1/2,k+1/2 − (uv)nj+1/2,k−1/2
=
unj+3/2,k − 2unj+1/2,k + unj−1/2,k
=
unj+1/2,k−1 − 2unj+1/2,k + unj+1/2,k+1
∆y
∆x2
+ O(∆y 2 )
+ O(∆x2 )
∆y 2
pnj+1,k − pnj,k
+ O(∆x2 )
=
∆x
,
(3.4)
,
+ O(∆y 2 )
(3.5)
,
.
(3.6)
(3.7)
V gornjih enačbah smo uporabili nekaj izrazov, ki na sl. 3.2 niso definirani. Izračunamo jih lahko z linearno interpolacijo, npr.
unj+1,k
=
unj+3/2,k + unj+1/2,k
.
(3.8)
2
Ko osnovno Navier-Stokesovo enačbo diskretiziramo na premaknjeni mreži, lahko
tlak in hitrostno polje izračunavamo izmenično prvega iz drugega in obratno, če
sistem napišemo v implicitni shemi
vn+1 = Fn − ∆t∇d P n+1
,
(3.9)
kjer F ≡ (F, G)T in ∇d pomeni diskretni gradientni operator. Polje F je definirano
z izrazom
uj+3/2,k − 2uj+1/2,k + uj−1/2,k
n
n
Fj+1/2,k = uj+1/2,k +∆t
Re∆x2
uj+1/2,k−1 − 2uj+1/2,k + uj+1/2,k+1
+
Re∆y 2
(3.10)
u2j+1,k − u2j,k
−
∆x
(uv)j+1/2,k+1/2 − (uv)j+1/2,k−1/2 n
.
−
∆y
Podobno definiramo tudi G. Če zahtevamo, da je divergenca hitrosti v n + 1-tem
časovnem koraku enaka nič, ∇d · vn+1 = 0, dobimo iz (3.9)
∇d · vn+1 = ∇d · Fn − ∆t∇2d P n+1 = 0
,
(3.11)
Poissonovo enačbo za tlak
n+1
pj,k−1 − 2pj,k + pj,k+1
pj−1,k − 2pj,k + pj+1,k
+
=
∆x2
∆y 2
Gj,k+1/2 − Gj,k−1/2 n
1 Fj+1/2,k − Fj−1/2,k
+
∆t
∆x
∆y
(3.12)
.
3 NUMERIČNE METODE
11
To rešujemo v n-tem časovnem koraku. Ko izračunamo hitrost vn+1 iz (3.9),
premaknemo markerje za ustrezen vn+1 ·∆t. Poissonovo enačbo za tlak smo reševali
z metodo SOR (simultaneous over relaxation) [17].
Algoritem za reševanje s sledilnimi delci lahko pretvorimo v naslednji enačbi
vn+1 = Fn − ∆t∇d P n+1
1
∇2d P n+1 =
∇d · Fn .
∆t
,
(3.13)
(3.14)
Po vsakem premiku markerjev moramo celicam popraviti status, ker markerji
skoznje potujejo. Začetni status celic je eden od začetnih pogojev, napotek za njegovo spreminjanje lahko združimo v nekaj korakov, s katerimi med premikanjem
posameznih markerjev zagotovimo, da
• če se marker premakne v prazno celico, ta postane robna,
potem, ko smo premaknili vse markerje, naredimo še en pregled celic in takrat
• vse robne celice, ki ne vsebujejo markerjev spremenimo v prazne,
nato pregledamo vse celice še enkrat in
• vse polne celice, ki jim je sosednja katera od praznih, spremenimo v robne,
• vse robne, ki nimajo praznih sosed, spremenimo v polne.
Sprememba statusa celic je zadnji del posameznega računskega cikla.
3.1.2
Preizkusni primer - odtekanje pravokotnega bloka tekočine
Kot testni primer smo pri obeh celičnih metodah uporabili odtekanje pravokotnega
bloka tekočine. To je znani problem podiranja jezu (broken dam problem). Je
poučna simulacija prostega robu v časovni odvisnosti, rezultati pa so dostopni v
literaturi [4]. Problem je enak diplomskemu primeru pri pogoju, da je hitrost spodnje plošče u0 = 0. Na začetku poskusa tekočina miruje med dnom in stranskima
stenama. Ob času t = 0 desno steno odmaknemo in pustimo, da tekočina prosto
izteče. Poskus delamo v večji posodi, od katere se tekočina odbije. Na sl. 3.3 je
prikazan časovni razvoj procesa za brezdimenzijska parametra Re = 2 in F r = 1,
kjer lahko opazimo tudi pojave, ki pričajo o slabosti metode. Ti so posledica končne
dimenzije celic, tako da sledilni delec, ki je za mrežno razdaljo oddaljen od stene,
že čuti njen vpliv. Na sl. 3.3, ki kaže položaj sledilnih delcev po devetih časovnih
enotah, vidimo, da se v spodnjem desnem kotu delci ne dotikajo stene, a se kljub
temu narivajo v navpični smeri. Podobno napako kaže tudi slika po dvanajstih
časovnih enotah. Stena na delce vpliva, čeprav na slikah ne vidimo, da bi se stene
dotaknili. Stik tekočine s steno je namreč določen z dolžino celice. Stene se tekočina
dotika, če je od nje oddaljena manj kot za dolžino celice. Večja natančnost stika
zahteva manjšo dimenzijo celice, kar pa drastično povečuje računski čas. Zato smo
prisiljeni narediti kompromis.
3 NUMERIČNE METODE
12
t=0
t=9
t=3
t = 12
t=6
t = 15
Sl. 3.3: Simulacija podrtja jezu z metodo sledilnih delcev za Re = 2,
1
F r2 =
1.
Pri gostoti mrežnih celic, ki smo jo uporabili, so dobljeni rezultati smiselni in
uporabni za nadaljno obravnavo, tipičen izračun pa na osebnem računalniku vzame
okrog pol ure časa (glej tudi poglavje o primerjavi numeričnih metod). Zaradi
premikanja sledilnih delcev je dodatna slabost metode, da stacionarnega stanja ne
moremo ugotavljati avtomatično. Postopek iteriranja zato ustavljamo interaktivno.
3 NUMERIČNE METODE
3.1.3
13
Rešitve
Rešitve diplomskega problema začnemo z istimi začetnimi pogoji kot testni problem, seveda pa ob času t = 0 začnemo premikati spodnjo ploščo s hitrostjo u0 od
desne proti levi (sl. 1.1). Izbrali smo tri primere (sl. 3.4–3.6).
Sl. 3.4: Stacionarna oblika za parametre toka Re = 0.75,
Sl. 3.5: Stacionarna oblika za parametre toka Re = 1,
1
F r2 =
1
F r2 =
Sl. 3.6: Stacionarna oblika za parametre toka Re = 0.3,
1.
1.
1
F r2 =
11.
3 NUMERIČNE METODE
3.2
Metoda robnih sledilnih delcev (SM)
3.2.1
Opis metode
14
Metoda MAC nudi informacijo tudi o delcih v notranjosti tekočine, ki pa nas ob
iskanju stacionarnih oblik roba ne zanimajo. Za take in podobne probleme so S.
Chen, D.B. Johnson in P.E. Raad [4] metodo MAC predrugačili v metodo robnih
markerjev (SM). Osnovna razlika med metodama je položaj markerjev, ki so pri
metodi SM razpostavljeni le po površini. Tudi robne markerje premikamo enako,
kot smo to počeli z markerji v metodi MAC, le zaradi morebitnega krčenja in
širjenja meje potrebujemo postopek za vstavljanje novih in brisanje odvečnih markerjev. Uporaba metode je smiselna pri problemih prostega roba zaradi manjše
prostorske in predvsem časovne zahtevnosti v primerjavi z metodo MAC. Pri
SM metodi uporabljamo enake gibalne enačbe kot pri metodi MAC, premaknjeno mrežo, rešujemo Poissonovo enačbo za tlak in premikamo markerje, ki so v
tem primeru le na površini tekočine. Ker opazujemo časovni razvoj oblike, robne
markerje razpostavimo po površini na začetku in pustimo, da jih tok tekočine premika, kakor mu narekujejo enačbe. Nobenega zagotovila nimamo, da ne bodo na
začetku robni delci, označeni z robnimi markerji, odplavali v notranjost. Do preklopa površin (sl. 3.8–3.10) pride, ker se tekočina lepi na stene in se vrti. Dimenzije
celice so velike v primerjavi z razdaljo med sosednjima legama poljubnega delca,
zato moramo robu, kot pri metodi MAC, verjeti z določeno toleranco. Zgodi se
lahko, da na slikah dobimo večkratno površino - več preklopljenih krivulj. To ne ustreza resničnemu dogajanju, ker z metodo robnih sledilnih delcev opazujemo, kam
potujejo delci, ki so bili na začetku robni. Ker smo na začetku procesa označili
robne delce, jim med časovnim razvojem sledimo tudi v notranjost, ko s seboj
potegnejo nase pripeti rob.
Bistvena razlika med metodo MAC in metodo SM je v drugačnem položaju markerjev in postopku sprememinjanja statusa celic. Poleg označb iz metode MAC
prazna, polna, robna in zidna, potrebujemo še novo označbo: zadeta, ki jo uporabimo le v postopku spreminjanja statusa celic kot pomožno označbo. Postopek
lahko opišemo v nekaj korakih:
• vse celice, ki vsebujejo (robne) markerje označimo kot zadete.
Potem izvršimo premik in dodamo nove markerje, kjer so razmiki med posameznimi
markerji preveliki - markerji imajo zdaj drugačen položaj kot so ga imeli, ko smo
označili zadete celice
• vse celice, ki zdaj vsebujejo markerje označimo trenutno za robne,
• za ostanek zadetih celic preverimo, če imajo kakšno sosedno celico s statusom
prazna in jih v tem primeru spremenimo v prazne, sicer pa v polne,
• na koncu še pregledamo, če katera od robnih celic nima prazne sosede in ji
spremenimo status v polno.
3 NUMERIČNE METODE
3.2.2
15
Preizkusni primer - odtekanje pravokotnega bloka tekočine
Zaradi primerjave z metodo MAC smo za testni primer uporabili podrtje pravokotnega bloka tekočine z istima parametroma (sl. 3.7).
t=0
t=9
t=3
t = 12
t=6
t = 15
Sl. 3.7: Odtekanje pravokotnega bloka tekočine z metodo robnih sledilnih delcev
za Re = 2, F1r2 = 1.
Metoda ima isto pomankljivost kot MAC, to je, da toge stene vplivajo na dani
delec že, ko so oddaljene za dolžino celice. To lahko opazimo na sliki, ki kaže obliko
tekočine pri t = 15. Zdi se, kot da so v tekočini zajede zunanjega medija, kar pa
ni res. Nastale so zaradi prevelike dolžine celice. Pri manjši mrežni razdalji bi bila
prosta površina na isti višini, zajeda pa manjša. Pri preklopu površin se na ta
način poveča prostornina tekočine zaradi numerične napake.
3 NUMERIČNE METODE
3.2.3
16
Rešitve
Zaradi primerjave z metodo MAC smo pognali simulacije danega problema za iste
tri pare vrednosti parametrov. Pomankljivost prepogibanja površin kot posledica
končno velike dolžine celice je tu še bolj očitna. Tako kot pri MAC tudi tu računanje
ustavljamo interaktivno.
Sl. 3.8: Stacionarna oblika za parametre toka Re = 0.75,
Sl. 3.9: Stacionarna oblika za parametre toka Re = 1,
1
F r2 =
1
F r2 =
Sl. 3.10: Stacionarna oblika za parametre toka Re = 0.3,
1.
1.
1
F r2 =
11.
3 NUMERIČNE METODE
17
3.3 Metoda končnih elementov (FEM)
3.3.1
Opis metode - Galerkinov pristop
Pri celičnih metodah smo izračunali vrednosti količin v danih točkah, medtem ko
pri FEM predpostavimo, da lahko iskano polje U diskretiziramo z
U =
N
X
uj ψj
,
(3.15)
j=1
kjer so uj koeficienti, ψj pa znane analitične funkcije, ki so različne od nič na
omejenem območju okrog točk j [7,19]. Funkcijam ψj pravimo testne funkcije.
Naj bo L(Ũ) = 0 vodilna enačba problema, ki ga rešujemo. Polje Ũ je točna rešitev,
če pa rešitev aproksimiramo s (3.15), potem enačba ne bo točno izpolnjena, ampak
bo
L(U ) = R ,
(3.16)
kjer je R ostanek. Rešitev problema so koeficienti uj , ki jih izračunamo z zahtevo,
da naj bo R minimalen v smislu
ZZ
Rϕi dΩ = 0
,
(3.17)
Ω
zato tudi
ZZ
Ω
L(U )ϕi dΩ = 0 ,
(3.18)
kjer so ϕi utežne funkcije. Te so definirane v točkah z neznanimi vrednostmi
spremenljivk i ∈ [1, M ]. Točke z neznanimi vrednostmi spremenljivk ui imenujemo
aktivne točke. Vseh točk je N , aktivnih pa M ≤ N . Enačba (3.17) pomeni,
da je ostanek R ortogonalen na utežne funkcije. Ker je utežnih funkcij M , kolikor
tudi neznanih koeficientov, dobimo z vstavljanjem vseh ϕi natanko M enačb za M
neznanih ui . Za Galerkinovo formulacijo je značilno, da so utežne funkcije enake
testnim, torej
ϕi ≡ ψi .
(3.19)
Testne funkcije praviloma izbiramo med odsekoma zveznimi polinomi istega nizkega
reda. Tako je ψj tista testna funkcija, ki ima v točki j vrednost 1, v ostalih točkah
0, vmes pa je interpolirana z ustreznim polinomom, ki določa red aproksimacije in
testne funkcije. Kot primer si oglejmo tri sosednje testne funkcije za enodimenzionalni problem, ko so le–te odsekoma zvezni polinomi prvega reda (linearne) (sl.
3.11).
3 NUMERIČNE METODE
1
ψj−1
18
ψj
ψj+1
....
....
.
.........
.. .. .
.. .. .
.... ........
....
....
....
....
....
.
.
.
.
.
.
.
.
.
.....
...
..
....
....
.....
.....
....
....
.... . ..
.
. ........
.
.
.
.
.
.
.
.
.
.
.
.
......
.....
..
....
.....
.
.... ...
.
.
.
.
.
.
.
.
.
.
.
.
.....
.
.....
.....
...
.
..
.
.
.
.....
.
.
..
.
.
.
.
.
.
.
.
.
.....
.
..
.
.....
.
.
.
...
.
.
.
.
.
.....
.. .
.
..
.
.
.
.
.
.
.
.
....
.
.....
....
....
xj−1
xj
xj+1
Sl. 3.11: Slika linearnih testnih funkcij za enodimenzionalni primer.
Analogno lahko definiramo testne funkcije za večdimenzionalne probleme. Med
enostavnimi liki (elementi), s katerimi pokrijemo domeno, so najbolj razširjeni
kvadrati, pravokotniki in trikotniki. Pravokotni elementi so uporabni predvsem pri
pravokotnih geometrijah, zato smo se pri diplomskem primeru odločili za trikotnike,
katerih mrežo lahko tesneje prilagodimo robovom nepravilnih oblik (sl. 3.12). Linearne testne funkcije za dvodimenzionalni problem so piramidalne funkcije, njihova
osnovna ploskev pa je sestavljena iz trikotnikov.
Sl. 3.12: Pri FEM na domeno razpnemo mrežo elementov.
Metoda ima svoje prednosti in slabosti. Med prednostmi je najpomembnejše dejstvo, da poznamo rob v robnih točkah, vmes pa ga lahko izračunamo z interpolacijo.
To je veliko bolj natančno kot pri celičnih metodah, če primerjamo taki simulaciji,
pri katerih sta število mrežnih točk (FEM) ter število celic (MAC, SM), enaki. Ker
imamo pri FEM točke na robu, lahko aproksimiramo tudi normalo in ukrivljenost
roba. Slabe strani se pokažejo pri divjih gibanjih tekočine, če se mreža premočno
deformira, in je za nadaljne račune neuporabna.
Pri celičnih metodah so različne začetne oblike vodile do istih stacionarnih stanj,
pri metodi FEM pa račun konvergira k stacionarni obliki le pri začetnih oblikah,
ki so blizu končni.
3 NUMERIČNE METODE
19
Stacionarno hitrostno polje izračunamo iz divergenčnega stavka (2.1) in stacionarne
Navier-Stokesove enačbe. Na začetku poglavja smo opisali, da diskretiziramo polja
u, v in p s testnimi funkcijami
X
X
X
u=
uj ψj , v =
vj ψj , p =
pj ψj .
(3.20)
j
j
j
Za diskretizacijo vseh polj bi radi uporabili linearne testne funkcije, ker je z njimi
izračun skalarnih produktov (3.18) bistveno enostavnejši. Čeprav večina strokovnjakov [12,20] priporoča uporabo testnih funkcij drugega reda za hitrostno in prvega
reda za tlačno polje, bomo s trikom iz poglavja 3.3.2 tlak približno izrazili s hitrostnim poljem, zato nas zdaj zanima le, kakšne naj bodo testne funkcije za hitrost.
Videli bomo, da so pri danih problemih rezultati, dobljeni z linearnim interpoliranjem vseh polj, popolnoma smiselni in primerljivi z rezultati celičnih metod.
Da bi prišli do Galerkinove formulacije stacionarnega Navier-Stokesovega problema, ukrepamo kot v (3.18). Stacionarno Navier-Stokesovo enačbo pomnožimo z
utežnimi funkcijami ψi in integriramo po Ω. Enako storimo z enačbo (2.1). Sledi
ZZ ux + vy ψi dS = 0 ,
(3.21)
Ω
ZZ 1
2
(u )x + (uv)y + px −
(uxx + uyy ) ψi dS = 0 ,
(3.22)
Re
Ω
ZZ
ZZ 1
1
2
(vxx + vyy ) ψi dS = − 2
ψi dS . (3.23)
(v )y + (uv)x + py −
Re
Fr
Ω
Ω
Rešitev stacionarnega sistema pri dani obliki domene predstavljata vektor hitrosti
v in tlačno polje p.
V stacionarnem stanju veljata na prostem robu pogoja (2.15) in (2.16). Pogoj
(2.16) uporabimo za izračun stacionarnega hitrostnega polja pri trenutni obliki
domene, pogoj (2.15) pa pri spreminjanju njene oblike. Obliko domene spremenimo
s premikanjem prostorobnih točk za
∆rj = (vj ·nj )nj · ∆t ; j ∈ ∂Ω ,
(3.24)
kjer je ∆t parameter, s katerim utežimo premikanje prostorobnih točk. Posebnost pri modeliranju prostega roba, je trojna točka; to je točka, v kateri se stikajo
tekočina, stena in zunanji medij. Gibanje trojne točke v realni tekočini je odvisno
od mnogih parametrov, ki jih v modelu ne moremo upoštevati. V takih primerih
se zatečemo k poenostavitvam. Najpreprostejši možnosti sta: trojnih točk ne premikamo, ali pa jih premikamo tako, da ostane kot med prostim robom in steno
konstanten. Pri konkretnem problemu je bil ustreznejši drugi način, zato je začetni
naklon roba v trojnih točkah enak naklonu v vseh naslednjih korakih.
Premikali smo tudi notranje točke, ker so deformacije prostega roba lahko tako
velike, da bi pri fiksnih pozicijah točk lahko privedle do defektne računske mreže,
3 NUMERIČNE METODE
20
pri kateri ena ali več notranjih točk leži izven domene, omejene z robnimi točkami.
Poleg tega je napaka metode najmanjša, če so trikotniki enakostranični. Zato smo
po vsakem premiku prostorobnih točk premaknili vsako notranjo točko v težišče
njenih sosed.
Iteriranje k stacionarnemu stanju ustavimo, ko je
X
vj ·nj ≤ δ ,
(3.25)
j∈∂Ω
kjer je δ neko majhno, pozitivno število. Ko je zadoščeno pogoju (3.25), so
spremembe roba dovolj majhne. Z zmanjševanjem parametra δ lahko poljubno
povečamo natančnost pri iskanju stacionarnega stanja.
Postopek iteriranja h končni obliki lahko razdelimo na štiri glavne dele:
1. izberemo začetno obliko domene
2. za trenutno obliko izračunamo stacionarno hitrostno polje pri pogoju, da je
tlak na robu enak nič
3. prosti rob premaknemo, kot nam narekuje prostorobni pogoj o tangencialnih
hitrostih
4. 2. in 3. korak ponavljamo, dokler niso premiki prostega roba dovolj majhni.
3 NUMERIČNE METODE
3.3.2
22
Povezava tlaka s hitrostnim poljem (penalty method)
Uporaba penalne metode (penalty nethod) je priljubljena predvsem v francoski
šoli numerične matematike [9,16,19]. Pravi, da lahko tlak približno izrazimo s
hitrostnim poljem kot
1
1
p = − (ux + vy ) = − (∇ · v)
ε
ε
,
(3.26)
kjer je parameter ε majhno, pozitivno število, ali
lim
ε→0
1
− (∇ · v)
ε
= p
.
(3.27)
Za boljšo predstavo zveze (3.26) skiciramo dva primera lokalnega hitrostnega polja.
.....
.....
....
....
....
.....
..... ..
.
....................
...
...
..................
.... ..
.....
.
.
.
...
....
....
.....
.....
.
.....
....
....
.....
.
.
.
.
.....
.. .....
............
..........
•
S
..........
............
.. .....
....
.....
....
....
....
.....
.
Sl. 3.13a: (∇ · v) < 0 ,
S
...........
...........
.. .....
....
.....
.....
....
....
....
..
.
.....
....
....
.....
.
.
.
..
....
.. ......
..........
..........
....
...................
..
.... .
....
.
.
.
...
....
....
.....
....
•
T
....
....
....
....
.....
....
.... .
... .
...................
....
sl. 3.13b: (∇ · v) > 0.
T
Puščice kažejo v smer vektorskega polja. Iz teorije polj sta nam tipa točk S in T
znana kot ponor in izvor, srečamo pa ju pri necevastih poljih. Taka polja imajo
divergenco različno od nič, lokalno pa so v točkah z negativno divergenco ponori,
v tistih s pozitivno pa izvori polja. Pri delu z nestisljivo tekočino lahko povežemo
tlak s hitrostnim poljem, če zahtevi ∇ · v = 0 popustimo in malenkostne odmike
od ničelne vrednosti, povežemo s tlakom, kot smo že pokazali v (3.26). Tako v
točki tipa S (sl. 3.13a) vidimo nadtlak, v točki tipa T (sl. 3.13b) pa podtlak, oba
sorazmerna ∇ · v. Vpliv poudarimo ali zanemarimo z ustrezno izbiro parametra
ε. Prednost metode je v tem, da smo tlak izrazili s hitrostnim poljem, in se tako
znebili ene od neznank. Oznako p bomo še vedno uporabljali, vendar bo od začetka
naslednjega poglavja naprej pomenila isto kot v (3.26).
3.3.3 Matrični členi
3 NUMERIČNE METODE
23
Pri reševanju z metodo končnih elementov je vmesni cilj matrični sistem A · x = b.
V matriki A so skriti deli enačb, ki so odvisni od vrednosti v aktivnih točkah.
V vektor b vpletemo robne pogoje in znane vrednosti, odvisno pač od danega
primera. Vektor x je sestavljen iz vrednosti spremenljivk v aktivnih točkah. Sistem
lahko točno nastavimo samo za linearne probleme, pri nelinearnih si pomagamo z
linearizacijo, v danem primeru v → v0 + v. V vektorju x so tako popravki hitrosti
v ali


u1
 .. 
 . 


u 
x =  N ,
(3.28)
 v1 
 . 
 .. 
vN
če je N število aktivnih točk. S testnimi funkcijami diskretiziramo spremenljivke
(3.20), prav tako pa tudi njihove odvode, npr.
X ∂(uj ψj )
X ∂ψj
X
∂u
=
=
=
uj
uj ψjx
∂x
∂x
∂x
j
j
j
.
(3.29)
Pri linearnih testnih funkcijah je definiran le prvi odvod, višji so enaki nič. Pri
drugih odvodih si pomagamo z Greenovo formulo
ZZ
Z
ZZ
(Px + Qy )ψdS =
(P, Q)·n ψds −
(Px ψx + Qy ψy )dS ,
(3.30)
Ω
∂Ω
Ω
kjer pomeni Ω domeno, ∂Ω njen rob, n = (nx , ny ) pa normalni vektor na ∂Ω v
smeri navzven Ω.
Najprej sistem lineariziramo in pustimo še v osnovni obliki. Iz
0
1 2 0
0
0
∇ (v + v) =
(v + v) · ∇ (v + v) + ∇p −
− F1r2
Re
sledi
1 2
(v · ∇)v0 + (v0 · ∇)v + ∇p −
∇ v =
Re
1
0
− (v0 · ∇)v0 − ∇p0 + ∇2 v0
1
− F r2
Re
,
(3.31)
(3.32)
kjer smo zanemarili člen (v · ∇)v ter na novo označili p= − 1ε ∇·v in p0= − 1ε ∇·v0 .
Prepišimo v Galerkinovo obliko (3.22), (3.23) in sproti razstavimo na komponente
ZZ 1
0
0
0
0
uux + vuy + u ux + v uy + px −
(uxx + uyy ) ψdS =
Re
Ω
ZZ 1 0
0
0 0
0 0
0
(u + uyy ) ψdS
− u ux − v uy − px +
Re xx
Ω
(3.33)
3 NUMERIČNE METODE
ZZ Ω
24
1
+
+ u vx + v vy + pu −
(vxx + vyy ) ψdS =
Re
ZZ 1 0
1
0 0
0 0
0
0
− u vx − v vy − py +
(v + vyy ) ψdS
−
F r2
Re xx
Ω
uvx0
vvy0
0
0
(3.34)
.
Nekajkrat naletimo na drugi odvod: v (3.33) pri uxx + uyy in px , saj p = − 1ε ∇·v.
Ustrezne člene integriramo per partes po (3.30), za primer dva iz (3.33)
ZZ
(uxx + uyy )ψdS =
Ω
Z
∂Ω
ZZ
(ux , uy )·n ψds −
Z
px ψdS =
Ω
∂Ω
pψnx ds −
ZZ
(ux ψx + uy ψy )dS
,
(3.35)
Ω
ZZ
pψx dS
.
(3.36)
Ω
Posvetimo se robnima integraloma v (3.35) in (3.36). Vanju vključimo robne
pogoje, ti pa so odvisni od posameznega primera. Pri prostem robu so aktivne
tudi prostorobne točke, vendar bomo v njih predpisali p = 0, da zadostimo enemu
od robnih pogojev. Tako od robnih integralov ostane v enačbi le tisti iz (3.35).
Izpeljali bomo matrične elemente za prostorobni problem, zato bomo v nadaljevanju izpeljave ustrezne robne integrale upoštevali.
Združimo (3.33)-(3.36):
ZZ ZZ
0
0
0
0
uux + vuy + u ux + v uy ψdS +
pψx dS−
Ω
Ω
Z
ZZ 1
1
(ux , uy )·n ψds +
ux ψx +uy ψy dS
−
Re ∂Ω
Re
Ω
ZZ ZZ
0 0
0 0
= −
u ux + v uy dS −
p0 ψx +
Ω
Ω
Z
ZZ 1
1
0
0
0
0
(u , u )·n ψds −
ux ψx + uy ψy dS
Re ∂Ω x y
Re
Ω
ZZ Ω
ZZ
+
+ u vx + v vy ψdS +
pψy dS−
Ω
Z
ZZ 1
1
(vx , vy )·n ψds +
vx ψx +vy ψy dS
−
Re ∂Ω
Re
Ω
ZZ ZZ
0 0
0 0
= −
u vx + v vy dS −
p0 ψy +
Ω
Ω
Z
ZZ 1
1
0 0
0
0
(v , v )·n ψds −
vx ψx + vy ψy dS
Re ∂Ω x y
Re
Ω
uvx0
vvy0
0
(3.37)
,
0
(3.38)
in za zgled do konca
RR izpeljemo
RRle levo stran
R (3.37), saj so ostali
RRčleni podobni. Štirje
0
značilni členi so
uux ψdS,
pψx dS, (ux , uy )·nψds ter (ux ψx + uy ψy )ψdS.
3 NUMERIČNE METODE
25
Vsakega posebej razvijemo in dobimo
X
ZZ
ZZ X
0
0
uux ψdS =
uj ψj
uk ψkx ψi dS =
Ω
Ω
=
j
=
j
Ω
k
XX
j
ZZ
Ω
pψx dS = −
∂Ω
j
ZZ
Ω
u0k
k
ε
Ω
1X
−
ε j
(ux , uy )·n ψds =
X Z
=
nxj
uj u0k ψj ψkx ψi dS =
ZZ X
1
1X
= −
ε j
Z
k
X X ZZ
ZZ
Ω
ψj ψkx ψi dS uj
(uj ψjx + vj ψjy )ψix dS =
ZZ
ψjx ψi dS uj
ZZ
ψjy ψi dS vj
(i,j)∈∂Ω
Ω
(3.40)
,
(3.41)
ψjy ψi ds uj
,
(ψjx ψix + ψjy ψiy )dS uj
.
j
X ZZ
j
,
j
X y Z
ψjx ψi ds uj +
nj
(ux ψx + uy ψy )dS =
(3.39)
(i,j)∈∂Ω
(3.42)
Zaradi dveh enačb v Navier-Stokesovem sistemu razpadeta vektorja x in b v dva
bloka, matrika A pa v štiri. Medsebojno prepleteni enačbi (3.37) in (3.38) povežemo
v matrični sistem tako, da v zgornjo polovico sistema vpišemo enačbo za komponento x (3.37), v spodnjo pa enačbo za komponento y (3.38). Če je število aktivnih
točk in neznanih vektorjev hitrosti N , je dimenzija vektorjev x in b 2N , dimenzija
matrike A pa 2N × 2N . Matrični sistem ima obliko



 
a1,1
...
a1,N
a1,N +1
...
a1,2N
b1
u1
  ..   .. 
 ..
..
..
..
..
..
 .   . 
 .
.
.
.
.
.




  uN 
 aN,1
...
aN,N
aN,N +1
...
aN,2N
  bN 
·

, (3.43)
=


 aN +1,1 . . . aN +1,N aN +1,N +1 . . . aN +1,2N   v1   bN +1 





  .. 
 .
..
..
..
..
..
  .. 

 ..
.
.
.
.
.
.
.
vN
b2N
a2N,1 . . . a2N,N
a2N,N +1 . . . a2N,2N
kjer so posamezni elementi
X
X
1 i,j
i,j,k 0
i,j,k 0
−
Fy vk + Gx,x
Fx uk +
ai,j = 2
ε
k
k
1
1
x i,j
y i,j i,j
i,j
n Hx + n Hy Gx,x + Gy,y
+
Re
(i,j)∈∂Ω Re
,
(3.44)
3 NUMERIČNE METODE
26
1 i,j
Fyi,j,k u0k + Gy,x
ε
X
ai,N +j =
k
,
1 i,j
,
Fxi,j,k vk0 + Gx,y
ε
k
X
X
1 i,j
i,j,k 0
i,j,k 0
aN +i,N +j = 2
Fy vk +
Fx uk + Gy,y
−
ε
k
k
1
1
x i,j
y i,j i,j
i,j
n Hx + n Hy Gx,x + Gy,y
+
Re
(i,j)∈∂Ω Re
aN +i,j =
bi =
X
j
X
(3.46)
,
(3.47)
( X
X
i,j,k 0
i,j,k 0
0
Fy uk vj0
Fx uk uj −
−
k
k
1
i,j 0
G i,j u0 + Gy,x
vj
ε x,x j
)
1
1
i,j
i,j
nx Hxi,j + ny Hyi,j Gx,x
+ Gy,y
u0j
u0j −
+
Re
Re
(i,j)∈∂Ω
−
1
bN +i = − 2 I i +
Fr
( X
X
X
i,j,k 0
i,j,k 0
0
Fx vk u0j
Fy vk vj −
−
j
(3.45)
k
k
1
i,j 0
i,j 0
− Gx,y uj + Gy,y vj
ε
)
1
1
i,j
i,j
+
nx Hxi,j + ny Hyi,j Gx,x
+ Gy,y
vj0
vj0 −
.
Re
Re
(i,j)∈∂Ω
(3.48)
,
(3.49)
Uvedeni znaki pomenijo:
Fxi,j,k
Fyi,j,k
i,j
Gx,x
i,j
Gy,x
i,j
Gx,y
=
ZZ
ψj ψkx ψi dS
,
(3.50)
ψj ψky ψi dS
,
(3.51)
Ω
=
ZZ
Ω
=
ZZ
ψjx ψix dS
,
(3.52)
ψjy ψix dS
,
(3.53)
ψjx ψiy dS
,
(3.54)
Ω
=
ZZ
Ω
=
ZZ
Ω
3 NUMERIČNE METODE
27
i,j
Gy,y
=
ψjy ψiy dS
,
(3.55)
Ω
Hxi,j
Hyi,j
I
ZZ
i
=
Z
ψjx ψi ds
,
(3.56)
ψjy ψi ds
,
(3.57)
∂Ω
=
Z
∂Ω
=
ZZ
ψi dS
.
(3.58)
Ω
Ostane nam še, da izračunamo skalarne produkte (3.50)–(3.58). Naj bo
X X i,j,k
ψi =
χi
,
j
(3.59)
k
kjer trojica (i, j, k) tvori trikotnik, χi,j,k
pa je ravnina nad trikotnikom (i, j, k) z
i
i,j,k
vrhom v točki i. Enačba za χi
je tako
χi,j,k
=
i
αi + βi x + γi y
2∆
,
(3.60)
kjer so parametri enačbe ravnine
αi = xj yk − xk yj
βi = yj − yk
,
(3.61)
,
(3.62)
γi = xk − xj ,
1 xi yi 1 ∆ =
1
x
y
j
j
2 1 xk yk (3.63)
.
(3.64)
Parameter ∆ je ploščina trikotnika (i, j, k). Pri izračunu si pomagamo s formulo
[18]
ZZ
l
n
m
m!n!l!
· 2∆ .
(3.65)
χk dS =
χj
χi
(m + n + l + 2)!
∆
Naravna števila (z ničlo) m, n in l pomenijo, kolikokrat nastopa ustrezna funkcija v
skalarnem produktu, še vedno pa velja, da (i, j, k) tvorijo trikotnik. Tako so od nič
različni integrali le tisti, katerih točke vrhov so točke trikotnika (lahko tudi enake).
Za zgled izračunajmo enega od integralov:
ZZ
1 1
1
i,j,k
χj χkx χi dS =
Fx
=
∆




∆
βk







ZZ
; i=j
; i=j



βk
βk
1
1
6
12
=
χj χi dS =
·
.
=



2∆
2∆ 
∆


 ∆ ; i 6= j 


 βk ; i 6= j 
12
24
(3.66)
3 NUMERIČNE METODE
28
Ko napolnimo matriko A in vektor b , z ustreznim programom za reševanje matričnega sistema, izračunamo vektor popravkov hitrostnega polja x , ki je vektor
popravkov hitrosti v. Postopek ponavljamo, dokler popravki niso dovolj majhni.
Tako dobimo stacionarno hitrostno polje pri dani obliki domene.
Matrika je redka (ima malo neničelnih elementov), ker ima vsaka od točk največ
šest sosed. Zato je v vsaki vrstici od štirih blokov matrike A šest neničelnih elementov plus diagonalni. To lastnost lahko izkoristimo, da prihranimo računski čas.
Namesto standardnih metod za reševanje matričnih sistemov (npr. Gaussova eliminacija) zato uporabimo program za reševanje redkih matričnih sistemov. Posegli
smo po programu sparse.pas iz programskega paketa Numerical Recipes in Pascal [17], ki uporablja algoritem za minimizacijo skalarne funkcije f (x) = |Ax − b|.
Namesto časovne zahtevnosti O(n3 ) (Gaussova eliminacija), je zahtevnost reda
10n, kjer je n dolžina vektorjev. Občuten prihranek v času pri razsežnih sistemih
nas stane nekaj natančnosti. V danem primeru je matrika A sestavljena iz števil,
ki se razlikujejo za k velikostnih redov, če je parameter penalne metode ε = 10−k .
Zato mora program teči z razširjenim tipom realnih števil, kar pa ni nujno slabost,
saj se zaradi posebnega načina prenosa podatkov na nekaterih računalnikih tako
še celo skrajša čas računanja.
3 NUMERIČNE METODE
3.3.4
29
Preizkusni primer 1: dvodimenzionalni vrtinec v votlini
Standardni in velikokrat predstavljeni [7,19] problem dvodimenzionalnega vrtinca v
votlini (driven cavity) nam je služil kot osnovni preizkus metode. Votlina kvadratne
oblike, v katero je do vrha natočena tekočina, je pokrita s premičnim pokrovom.
Tega vlečemo s hitrostjo u0 ter tako tekočino spravimo v gibanje (sl. 3.14).
u = u0 = 1, v = 0
...
................................................................................
u = 0,
v=0
u = 0,
v=0
tekočina
u = 0, v = 0
Sl. 3.14: Skica modela dvodimenzionalnega vrtinca v votlini z robnimi pogoji.
Votlina nima prostega roba, zato v matričnih členih ni robnega integrala iz (3.35).
Brez težnega pospeška g0 = 0 je tudi F1r2 = 0. V vseh primerih je hitrost zgornje
plošče u0 = 1.
Na slikah 3.15a in 3.15b sta primera z Reynoldsovima številoma Re = 1 in Re =
600. Pri Re > 100 se v spodnjem desnem kotu pojavi še en vrtinec velikosti stotinke
ploščine votline, ki kroži v nasprotni smeri.
|v|max = 0.407
Sl. 3.15a: Re = 1, ε = 0.01,
|v|max = 0.359
sl. 3.15b: Re = 600, ε = 10.
Sliki 3.15a in 3.15b lahko primerjamo z zgledi, ki jih najdemo v [19], kjer so prav
tako uporabljene linearne testne funkcije.
3 NUMERIČNE METODE
3.3.5
30
Preizkusni primer 2: odtekanje eliptičnega bloka tekočine
Podobno kot v poglavjih 3.1.2 in 3.2.2, le z drugačno začetno obliko, pustimo da
tekočina prosto izteče. Poleg dodatnega robnega integrala iz (3.35) se ta primer
od prejšnjega (driven cavity), razlikuje tudi po tem, da potrebujemo postopek za
premikanje roba. Pri dani obliki domene izračunamo stacionarno hitrostno polje s
tlakom, ki ima na robu vrednost nič. Potem prostorobne točke premaknemo, kot
nam narekujejo hitrosti na robu (3.24).
S parametrom ∆t utežimo premikanje roba. Z izbiro prevelikega parametra se
lahko mreža premočno deformira že pri prvem premiku. Če pa je ∆t premajhen,
potrebujemo veliko število iteracij, da izračunamo celoten potek. Preverili smo, da
je za vrednosti ∆t, ki so manjše od neke mejne vrednosti, časovni potek približno
enak.
Po premiku prostorobnih točk premaknemo še notranje točke z zahtevo, da ostanejo
čim bolj enakomerno razmetane po domeni – vsaka točka naj bo v težišču svojih
sosed. Pri odtekanju tekočine stacionarne oblike, kjer velja (3.25), ne bomo dosegli.
Sledili bomo robu, dokler se mreža preveč ne deformira. Na sl. 3.16 so oblika
začetnega stanja, polje hitrosti v začetnem stanju tik pred prvim premikom mreže,
in polje hitrosti tik pred štirimi naslednjimi zaporednimi premiki mreže.
3 NUMERIČNE METODE
Sl. 3.16: Odtekanje eliptičnega bloka tekočine, Re = 2,
31
1
F r2
= 1.
3 NUMERIČNE METODE
3.3.6
32
Rešitve
K diplomskemu primeru preidemo od drugega primera s premikanjem spodnje
plošče in merjenjem napake normalnih hitrosti (2.15). Ko pade |∆r| pod dovolj
majhno mejno vrednost, smo z rešitvijo zadovoljni. Na tok tekočine poleg Re in
1
F r 2 vplivata tudi parameter ε in začetna oblika kaplje, ki je podana z razmerjem polosi elipse. Začetna oblika mora biti dovolj blizu končnemu stanju, da ga
z zaporednimi premiki roba sploh dosežemo. Predstavljamo dva primera, ki sta
narejena tudi z metodo sledilnih delcev.
Sl. 3.17: Oblika po dolgem času za Re = 1,
1
F r2
Sl. 3.18: Oblika po dolgem času za Re = 0.3,
= 1.
1
F r2
= 11.
4
PRIMERJAVA METOD
Vse tri metode smo primerjali po prostoru in času, ki ga vzamejo računalniškemu
pomnilniku. Predstavili smo še glavne prednosti in slabosti ter se potem odločili
za metodo, ki jo bomo uporabili pri pregledu spektra parametrov. Od metode
zato zahtevamo, da deluje čez široko območje različnih vhodnih vrednosti ter da je
časovno čim manj zahtevna.
4.1
Prostorska zahtevnost
Računski mreži sta pri metodah MAC in SM enaki, zato primerjava ni sporna,
pri FEM pa za primerjavo z delčnima metodama (število celic) vzamemo število
mrežnih točk. Oblika roba je določena do celice natančno pri delčnih metodah, ter
do mrežne točke natančno pri FEM metodi.
Pri metodi MAC je na računsko domeno, ki obsega celotno pravokotno posodo,
napeta mreža velikosti JM AX ×KM AX (sl. 3.2). Pri računanju potrebujemo štiri
polja realnih števil: dve hitrostni, tlačno, eno pomožno polje za izračun tlačnega
polja (3.12), ter črkovno polje z označbami statusa celic. V celicah, kjer je na
začetku tekočina, imamo tudi sledilne delce. Simulacije so narejene z d (16 do 36)
delci na celico. Če začetna oblika pokriva približno četrtino celotne domene, je
velikost polja sledilnih delcev
Nsled.
delci
=
1
JM AX × KM AX × d .
4
(4.1)
Eden od trenutno zmogljivejših osebnih računalnikov je zmogel premleti izračune
do velikosti mreže JM AX × KM AX = 30 × 30. Izračun prostorsko bolj omejuje program PASCAL 6.0, kot pa računalniški pomnilnik. Program namreč ni
sposoben izkoristiti vsega spomina, ki je na voljo v računalniku. Če programi
tečejo v razširjenem spominu, smo omejeni šele pri nekajkrat večjih mrežah.
Tudi pri SM metodi potrebujemo štiri polja realnih in eno polje črkovnih vrednosti.
Nekaj pomnilniškega prostora prihranimo zaradi manjšega števila markerjev. Ker
so robni, so zdaj namesto šestnajstih v navadni robni celici le štirje, v ogliščni celici
pa jih je sedem. Število robnih sledilnih delcev v navadni robni celici označimo z d.
33
4 PRIMERJAVA METOD
34
Tako je število markerjev za primer, ko sta stranici pravokotnika približno enako
dolgi
Nrob.
sled. delci
=
1
(JM AX + KM AX) × d + 4 · (2 × d − 1) .
2
(4.2)
Produkt 4 · (2 × d − 1) pomeni število robnih markerjev v štirih ogliščnih celicah.
Razlika (4.1) in (4.2) je ves prostorski prihranek, ki ga pridobimo s prehodom od
MAC k SM metodi.
Pri programu FEM največ pomnilniškega prostora potrebujemo za dve hitrostni
polji dimenzije N , matriko sistema dimenzije 2N × 2N in vektor desnih strani
dimenzije 2N . Prostor jemlje še množica pomožnih polj: polje trikotnikov, aktivnih
točk, itd. Če primerjamo MAC ali SM program z danim številom polnih in robnih
celic ter FEM program z enakim številom aktivnih točk, slednji zahteva okrog 50%
več prostora.
4.2
Časovna zahtevnost
Vrednosti karakterističnih časov so bile ocenjene na osebnem računalniku s procesorjem INTEL 80486, matematičnim koprocesorjem in frekvenco sistemske ure 33
MHz.
Pri metodi MAC porabi program največ časa pri računanju tlačnega polja iz Poissonove enačbe (3.12) s podprogramom SOR. Ker opazujemo časovni razvoj kaplje,
je čas, potreben za izračun stacionarnega stanja, odvisen od parametrov toka, groba
ocena njegove karakteristične vrednosti je okrog 30 minut.
Tudi metoda SM porabi največ časa pri postopku SOR, v primerjavi z MAC pa
prihrani nekaj časa pri premikanju sledilnih delcev. Pri istih pogojih porabi približno 70-80% časa MAC simulacije za dosego stacionarnega stanja.
Pri metodi končnih elementov je časovno najpotratnejša procedura izračuna redkega matričnega sistema sparse.pas, katere časovna zahtevnost močno raste z
dimenzijo sistema. Ker iteriranja ne ustavljamo interaktivno, je primerjava z metodama sledilnih delcev vprašljiva. Groba ocena značilnih časov iteracije je približno
60 minut.
4.3
Ostale prednosti in slabosti
Pri metodi MAC gosto polje sledilnih delcev daje nazoren opis prostega roba,
čeprav ne smemo pozabiti, da bi pravo obliko dobili šele z zmanjšanjem dolžine
4 PRIMERJAVA METOD
35
celice. Metoda deluje za širok prostor parametrov kaplje, kar je ugodno pri pregledovanju spektra stacionarnih stanj, čeprav slednjih ne moremo zaznati z ustreznim
podprogramom in smo prisiljeni program ustaviti interaktivno.
Če bi bile pri metodi SM robne točke ves čas procesa na robu, bi metoda nudila
idealno sliko prostega roba, tako pa nas sili k različnim interpretacijam zajed in
zavihkov, ki nastajajo med časovnim razvojem kaplje. Ostale slabosti so omenjene
že pri metodi MAC.
Prednost metode FEM je v sistematiziranem iskanju stacionarnega stanja, ki je
jasno določeno s prostorobnim pogojem (2.15), vendar ga dosežemo le za ozek
prostor parametrov, kar je glavna slabost metode.
Vse tri metode imajo precejšnjo prostorsko in časovno zahtevnost, tako da so bile
pri izbiri najustreznejše metode za preiskovanje spektra parametrov odločilne ostale
prednosti in slabosti. Odločili smo se za metodo MAC: nazornejša je od metode
SM ter hitrejša, stabilnejša in splošnejša od metode končnih elementov.
5
RAZDELITEV PROSTORA
PARAMETROV
5.1
Posamezni zgledi
Pri numeričnem računanju si ne moremo privoščiti vogala z neskončnima stranicama, zato smo kapljo, kot že v prejšnjih poglavjih, zaprli v oglato posodo.
Re = 0.01
1
F r 2 =0.01
Re = 0.1
1
F r 2 =10
Re = 0.01
1
F r 2 =100
Re = 0.1
1
F r 2 =100
Sl. 5.1: Oblike kaplje po dolgih časih. V levem stolpcu za Re=0.01, v desnem za
Re=0.1. Vse oblike so stacionarne, razen spodnje desne s parametroma
Re=0.1 in F1r2 =100, kjer posamezni markerji kažejo živahno dogajanje
v repu tekočine (primerjava z velikostjo celice na sl. 3.1). Sliki v levem
stolpcu za Re = 0.01 očitno kažeta kvazi-stacionarno stanje kaplje z zelo
veliko viskoznostjo.
36
5 RAZDELITEV PROSTORA PARAMETROV
37
Iz enačbe (2.7) vidimo, da je dinamika kaplje odvisna od parametrov Re in F1r2 , ki
sta določena s karakteristično dolžino L0 , hitrostjo u0 , z materialnimi lastnostmi
preko kinematične viskoznosti ν = ηρ in s poljem sile teže preko g0 .
Re = 5
1
F r 2 =0.5
Re = 5
1
F r 2 =1
Sl. 5.2: Obliki kaplje po dolgem času za Re=5. Leva slika kaže skoraj stacionarno
obliko, pri desni se največje spremembe prostega roba dogajajo v repu.
Re = 1
1
F r 2 =1
Re = 2
1
F r 2 =1
Re = 1
1
F r 2 =10
Re = 2
1
F r 2 =10
Sl. 5.3: Končne oblike kaplje. V levem stolpcu Re=1, v desnem pa Re=2 . Gornji
sliki sta sliki stacionarne oblike.
5 RAZDELITEV PROSTORA PARAMETROV
38
Sistematičnega pregledovanja smo se lotili pri majhnih Re in F1r2 ter nadaljevali
postopoma do velikostnih redov 1000. Zaradi že omenjenih slabosti numerične
metode MAC (interaktivno ustavljanje) bomo končno razdelitev spektra ponudili
le kot namig, v katerem delu lahko pričakujemo stacionarna stanja. Na sl. 5.15.3 smo prikazali nekaj končnih oblik. To so oblike po dovolj dolgih časih, ko
se je tekočina polegla, pri nestacionarnih končnih stanjih pa ponavadi opazujemo
turbuletno gibanje v repu tekočine (glej npr. sl. 5.3, Re = 1, F1r2 = 10).
5.2
Razdelitev prostora
S pomočjo podatkov iz prejšnjega podpoglavja lahko približno nakažemo tiste dele
spektra, v katerih pričakujemo stacionarna stanja kaplje. Poiskali smo tudi nekaj
posebnih primerov gibanja in jih zabeležili v prostoru parametrov. V zgledih
nastopajo različne snovi, vendar vse v primerih, kjer veljajo za nestisljive, sile
površinske napetosti pa so zanemarljive zaradi velikosti sistema. Primere ločimo
med seboj po tekočinah, ki v njih nastopajo, in so z imenom tekočine označeni tudi
na sl. 5.5.
voda: Kapljica vode je na tekočem traku, ki izginja pod steno, kot na sl. 2.1. Njen
−2 kg
presek je S = 1cm2 , hitrost plošče pa u0 = 1 cm
s . Če je η = 10
ms , potem
1
F r 2 = 1000 in Re = 100.
2
glicerin: Enako kot pri vodi, le ν = 6.7 cm
s , zato
1
F r2
= 1000, Re = 0.15.
med: Lahko si mislimo kapljo meda na tekočem traku, kot v prejšnjih primerih, ali
pa, da ga z ostrim predmetom strgamo s podlage. Z istima S in u0 dobimo
1
F r 2 = 1000. Viskoznost medu je v okolici sobne temperature močno odvisna
kg
kg
od temperature . Vzeli smo η ≈ 10 ms
in ρ = 1000 m
3 , zato Re ≈ 0.01.
testo: Predstavljajmo si valjasto kepo testa z osnovno ploskvijo S = 1dm2 , ki jo
valjamo s hitrostjo u0 = 1 dm
s . Valjar je pravokoten in pred sabo nariva testo.
kg
kg
−3
Za viskoznost testa smo vzeli η ≈ 104 ms
in ρ ≈ 1000 m
.
3 . Sledi Re ≈ 10
asfalt: Še preden se asfalt strdi, ga cestni delavci razmažejo s cestnim valjarjem.
Dogajanje je podobno valjanju testa. Karakteristični vrednosti sta L0 = 0.1m
kg
kg
1
in u0 = 0.1 m
s . Za vrednosti η ≈ 100 ms in ρ ≈ 2000 m3 dobimo F r 2 = 100,
Re ≈ 0.2.
sneg: Tudi pluženje snega lahko vzamemo za primer, če si mislimo, da je plug ravne
oblike. Valj preseka S = 1m2 , ki se nabere pred plugom, rinemo s hitrostjo
1
u0 = 1 m
s , zato F r 2 = 10. Poznamo različne vrste snega. Če za viskoznost in
kg
kg
gostoto vzamemo η ≈ 100 ms
in ρ = 200 m
3 , je Re ≈ 2.
5 RAZDELITEV PROSTORA PARAMETROV
39
magma: Zadnji primer smo našli v geologiji. Med živahnim gibanjem tektonskih plasti
pod zemeljskim površjem se ena plast zariva pod drugo. Na prvi plasti je
ogromna kepa vroče magme, ki se dviguje ob navpičnem robu druge plasti.
kg
Gostota magme je približno ρ = 4000 m
3 , viskoznost pa je močno odvisna od
kg
temperature, za naš primer naj bo η ≈ 107 ms
. Tako dobimo F1r2 = 104 in
Re ≈ 10−3 .
4
3
2
log10
1
F r2
1
.
.....
.........
...................
.......
.........................
.........
.........................
......................
...........
...........
.......
..
magma
..
.....
.......
..............
..........
...................
....................
.
...............
.......................
..........................
.............................. .
..................
.................... . .
...............
..............
..............
.............. . . .
.................
........... . . . .
........
.
.
.
.
. ..... . . . . .
..
.
.
.
.
.
.
.
.
.
.
.
....
....
. . . . . ...............
................
.
.
.
.
.
.
.
.
.
.
.
.
. . . . . ..............
......................
. . . . ........................
.........................
. . . ...................
.................
.................
. . ............................
.......
. . . ......................
....................
. . . .................
.........
. . . .........
.........
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . ..........
. . ....................
. .........................
...........................
........
..............................
.........
..............................
.
.
............................ .
................... . .
....................... . . .
............. . . . .
...........
....... . . . . .
..
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
med
voda
testo
nestacionarna
stanja
asfalt
glicerin
0
−1
sneg
stacionarna
stanja
−2
−3
−2
−1
0
1
2
3
4
log10 Re
Sl. 5.5: Razmejitev prostora (Re, F1r2 ) na vrednosti, ki vodijo do stacionarnega
stanja, in na tiste, ki tega ne dosežejo.
6
ZAKLJUČEK
Predstavili smo tri metode za reševanje problemov nestisljive viskozne tekočine s
prostim robom: dve standardni celični in nov algoritem, ki temelji na relaksaciji
trikotne mreže končnih elementov. Metoda končnih elementov se je izkazala kot
uporabna predvsem pri preizkusnih primerih, zaradi občutljivosti na začetne pogoje
pa je slabše služila pri konkretnem primeru iskanja stacionarne oblike. Edina je,
pri kateri smo lahko avtomatično določili, kdaj so hitrosti na robu tangencialne, in
tako izvedeli, kako daleč od stacionarnega stanja smo pri dani iteraciji.
Metoda sledilnih delcev se je med naštetimi izkazala kot najustreznejša za pregled
spektra stanj. V predzadnjem poglavju smo razdelili prostor parametrov in vanj
vključili nekaj zgledov, ki so vsaj približno podobni glavnemu primeru. Med zgledi
je nekaj eksotov, pri katerih je kinematično viskoznost težko najti v tabelah, ali
sploh kako drugače določiti.
Upamo, da bodo dobljeni rezultati potrjeni tudi z eksperimenti. Problem prehoda
v turbulenco je pred nekaj leti spet prišel med popularnejše zaradi prodorov v
eksperimentalni in teoretični fiziki tekočin. Ker scenariji prehoda od urejenega
do turbulentnega gibanja še niso popolnoma raziskani, ostane za nadaljne delo
predvsem izpopolnitev metod, da bomo lahko preverili, če tudi numerični rezultati
vodijo do višjih stabilnih stanj, kot so nihanja tekočine med dvema in več oblikami
roba.
40
ZAHVALA
Doc. dr. B. Šarler, dipl.ing. in A. Košir, dipl.ing., Inštitut Jožef Štefan,
R4, sta mi nudila dostop do strokovne literature in nepogrešljive računalniške opreme. Del naloge je bil po zaslugi prof. dr. V. Valenčiča,
dipl.ing., opravljen tudi na računalnikih Laboratorija za Biokibernetiko
na Fakulteti za Elektrotehniko in Računalništvo v Ljubljani, ki ga vodi
prof. dr. L. Vodovnik, dipl.ing. Vsem se za pomoč in dobro voljo iskreno
zahvaljujem.
41
LITERATURA
[1] Dictionary of Physics, Valerie Illingworth (ed.), Penguin Books, London, 1991.
[2] Batagelj, V., Golli, B., TEX, DFMA, Ljubljana, 1990.
[3] Carey, G.F., Krishnan, R., Penalty Approximation, Iteration, and
Continuation for Navier–Stokes Problems, Finite Elements in Fluids,
Vol. 6, John Wiley & Sons Ltd., 1985.
[4] Chen, S., Johnson, D.B., Raad, P.E., The Surface Marker Method,
Wrobel, L.C., Brebbia, C.A. (eds.), Computational Modelling of Free
and Moving Boundary Problems, Fluid Flow, Vol. 1, Computational
Mechanics Publications and Walter de Gruyter, Southampton and
Berlin, 1991.
[5] De Sampio, P.A.B., A Petrov–Galerkin Formulation for the Incompressible Navier–Stokes Equations Using Equal Order Interpolation
for Velocity and Pressure, Int. J. Numer. Methods Eng., Vol. 31,
1991.
[6] Fetter, A.L., Walecka, J.D., Theoretical Mechanics of Particles and
Continua, McGraw–Hill, New York, 1980.
[7] Fletcher, C.A.J., Computational Techniques for Fluid Dynamics, Springer–
Verlag, Berlin, 1988.
[8] Floryan, J.M., Rasmussen, H., Numerical Methods for Viscous Flows
with Moving Boundaries, Appl. Mech. Rev., Vol. 42, 1989.
[9] Fortin, M., Glowinski, R., Augmented Lagrangian Methods: Applications to the Numerical Solutions of Boundary-Value Problems,
North–Holland, Amsterdam, 1983.
[10] Harlow, F.H., Shannon, J.P., The Splash of a Liquid Drop, J. Appl.
Phys., Vol. 38, 1967.
42
LITERATURA
43
[11] Harlow, F.H., Welch, J.E., Numerical Calculation of Time–Dependent
Viscous Incompressible Flow of Fluid with Free Surface, Physics of
Fluids, Vol. 8, 1965.
[12] Hughes, T.J.R., et al., Finite Element Analysis of Incompressible
Viscous Flow by the Penalty Function Formulation, J. Comp. Phys.,
Vol. 30, 1979.
[13] Hughes, W.F., Brighton, J.A., Fluid Dynamics, Schaumm’s Outline
Series, McGraw–Hill, New York, 1967.
[14] Koškin, N.I., Širkevič, M.G., Priročnik elementarne fizike, Tehniška
založba Slovenije, Ljubljana, 1984.
[15] Landau, L.P., Lifshitz, E.M., Course of Theoretical Physics, Vol. 6,
Pergamon Press, 1987.
[16] Pelettier, D., et al., Are FEM Solutions of Incompressible Flows Really Incompressible? (Or How Simple Problems Can Cause Headaches!),
Int. J. Numer. Methods, Vol. 9, 1989.
[17] Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., Numerical Recipes in Pascal, Cambridge University Press, Cambridge,
1990.
[18] Strang, G., Fix, G.J., An Analysis of the Finite Element Method,
Englewood Cliffs, 1973.
[19] Temam, R., Navier–Stokes Equations, North–Holand, Amsterdam,
1984.
[20] Zienkiewicz, O.C., The Finite Element Method, McGraw–Hill, London, 1977.
[21] Zienkiewicz, O.C., Wu, J., Incompressibility Without Tears – How to
Avoid Restrictions of Mixed Formulation, Int. J. Numer. Methods
Eng., Vol. 32, 1991.