Navierových-Stokesových rovnic. Numerical Solution of Navier-Stokes equations using the finite element method

VSˇB – Technicka´ univerzita Ostrava Fakulta elektrotechniky a informatiky Katedra informatiky

Numericke´ rˇesˇenı´ Navierovyćh-Stokesovyćh rovnic pomocı´ metody konecˇnyćh prvku˚ Numerical Solution of Navier-Stokes equations using the finite element method

2011

Martin Hasal

Prohlasˇuji, zˇe jsem tuto bakala´rˇskou praći vypracoval samostatneˇ. Uvedl jsem vsˇechny litera´rnı´ prameny a publikace, ze kteryćh jsem cˇerpal.

V Ostraveˇ 5. kveˇtna 2011

.............................

Za ochotu a trpeˇlivost bych ra´d podeˇkoval vedoucı´mu diplomove´ praće panu Prof. RNDr. Pavlovi Burdovi, CSc. Za trpeˇlivost a poskytovańı´ cennyćh rad deˇkuji svy´m trˇem ucˇitelu˚m doc. RNDr. Jirˇ´ımu Bouchalovi, Ph.D., Prof. RNDr. Radimovi Blahetovi, CSc. a Doc. Ing. Daliborovi Luka´sˇovi, Ph.D. Da´le bych chteˇl podeˇkovat sve´ rodineˇ za umozˇneˇnı´ studia na vysoke´ sˇkole. Me´ prˇ´ıtelkyni a jejı´ rodineˇ za pomoc prˇi realizaci te´to praće. V neposlednı´ rˇadeˇ bych ra´d podeˇkoval cele´ Katedrˇe aplikovane´ matematiky, zvla´sˇt’panu Prof. RNDr. Zdenˇkovi Dosta´lovi, DSc. za umozˇneˇnı´ psanı´ praće na vy´sˇe uvedene´ te´ma.

Abstrakt V praći jsou prˇedstaveny za´klady numericke´ho rˇesˇenı´ Navierovyćh-Stokesovyćh rovnic pomocı´ metody konecˇnyćh prvku˚. K diskretizaci oblasti, v nı´zˇ hleda´me rˇesˇenı´, pouzˇ´ıva´me smı´sˇenou metodu konecˇnyćh prvku˚, a to pomocı´ Q2 −Q1 prvku˚. U teˇchto prvku˚ je oveˇrˇena take´ jejich stabilita. Prˇi numericke´m rˇesˇenı´ Navierovyćh-Stokesovyćh rovnic vycha´zı´me ze Stokesovy u´lohy, ktera´ je jejı´m linearizovany´m prˇ´ıpadem. Na´sledneˇ pomocı´ Picardovy a Newtonovy iteracˇnı´ metody dopocˇ´ıta´me rˇesˇenı´ Navierovyćh-Stokesovyćh rovnic. Klıćˇova´ slova: Stokesovy rovnice, Navierovy-Stokesovy rovnice, smı´sˇena´ metoda konecˇnyćh prvku˚, Newtonova iteracˇnı´ metoda, Picardova iteracˇnı´ metoda, diplomova´ praće

Abstract In the work are introduced the basics of the numerical solution of the Navier-Stokes equations using the finite element methods. For domain discretization is using the mixed finite element methods by means of Q2 − Q1 elements, where the stability was verified. In the numerical solution of the Navier-Stokes equations was used Stokes equations, which is a linearized case. After that we use a Picard and Newton iteration for the solving of the Navier-Stokes equations. Keywords: Stokes equations, Navier-Stokes equations, mixed finite element methods, Newton iteration, Pickard iteration, master thesis

Seznam pouzˇityćh zkratek a symbolu˚ MKP Ω Γ = ∂Ω Γ0,...,n t u p ν φ, ψ (·, ·) a(·, ·), b(·) n = n(x) RN

– – – – – – – – – – – – –

Metoda konecˇnyćh prvku˚ Oblast-omezena´ otevrˇena´ mnozˇina Hranice oblasti Cˇa´st hranice, ∪ni=0 Γi = ∂Ω Cˇas Rychlost Tlak Viskozita (vazkost) Ba´zove´ funkce Skala´rnı´ soucˇin Bilinea´rnı´ forma Jednotkovy´ norma´lovy´ vektor N-rozmeˇrny´ rea´lny´ prostor

1

Obsah 1

´ vod U 1.1 Trochu historie, fyziky a odvozovańı´ . . . . . . . . . . . . . . . . . . . . . .

2

´ vod do metody konecˇnyćh prvku˚ U 2.1 Jednorozmeˇrna´ u´loha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Dvourozmeˇrna´ u´loha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Metoda konecˇnyćh prvku˚ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Galerkinova metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Smı´sˇena´ formulace a smı´sˇena´ metoda konecˇnyćh prvku˚ pro elipticke´ rovnice

9 9 12 15 16 23

3

Stokesova u´loha 3.1 Prˇedstavenı´ Stokesovy rovnice . . ´ lohy . . . . . . . . . . . . . . . . . 3.2 U 3.3 Slaba´ formulace . . . . . . . . . . . 3.4 Smı´sˇena´ metoda konecˇnyćh prvku˚ 3.5 Podmıńky rˇesˇitelnosti . . . . . . . 3.6 Stabilnı´ prvky Q2 − Q1 . . . . . . . ˇ esˇenı´ syste´mu rovnic . . . . . . . 3.7 R

4

3 3

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

32 32 34 36 38 40 42 45

Navierova-Stokesova u´loha 4.1 Rovnice konvekce-difuze . . . . . . . . . . . . . . . 4.2 Bezrozmeˇrny´ tvar . . . . . . . . . . . . . . . . . . . 4.3 Slaba´ formulace a linearizace . . . . . . . . . . . . 4.4 Bifurkacˇnı´ analy´za a teorie stability . . . . . . . . . 4.5 Newtonova a Picardova metoda . . . . . . . . . . 4.6 Aproximace smı´sˇenou metodou konecˇnyćh prvku˚ ´ lohy k rˇesˇenı´ . . . . . . . . . . . . . . . . . . . . . 4.7 U

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

47 47 48 49 50 53 54 56

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

5

Za´veˇr

58

6

Literatura

59

Prˇı´lohy A Apendix A.1 Opera´tory . . . . . . A.2 Vlastnosti opera´toru A.3 Veˇty . . . . . . . . . . A.4 Prostory . . . . . . .

60 . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

61 61 61 62 62

2

Seznam obra´zku˚ 1 2

3 4 5 6 7 8 9 10

11 12 13 14 15 16 17 18

Pu˚sobenı´ sil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Linea´rnı´ rozlozˇenı´ rychlosti v prˇi lamina´rnı´m proudeˇnı´ mezi dveˇma deskami (vlevo) Rychlostnı´ profil a tecˇne´ napeˇtı´ (vpravo) . . . . . . . . . . . . . . . . . . . . Ba´zova´ funkce (φi (xj ) = δij ) . . . . . . . . . . . . . . . . . . . . . . . . . . . Integracˇnı´ oblast, Γ0 Dirichletova hranicˇnı´ podmıńka a Γ1 Neumannova hranicˇnı´ podmıńka . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Q1 ba´zova´ funkce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Transformace sourˇadnic Q1 prvku . . . . . . . . . . . . . . . . . . . . . . . Druhy cˇtyrˇsteˇnovyćh elementu˚ . . . . . . . . . . . . . . . . . . . . . . . . . Ba´zove´ funkce Mh (hornı´ obra´zek), Xh (spodnı´ obra´zek) . . . . . . . . . . Zobrazeni norma´lovyćh smeˇru˚ na troju´helnı´ku . . . . . . . . . . . . . . . . Obecny´ troju´helnı´k T (vlevo), mi uzly ve strˇedech stran, ni jednotkove´ norma´love´ vektory smeˇrˇujıćı´ ven z T . Referencˇnı´ troju´helnı´k Tˆ (vpravo), mi uzly ve strˇedech stran, ni jednotkove´ norma´love´ vektory, i = 1, 2, 3. . . ˇ esˇenı´ rychlosti Stokesovy u´lohy, prˇ´ıklad 1: vtok-vy´tok . . . . . . . . . . . R ˇ esˇenı´ tlaku Stokesovy u´lohy, prˇ´ıklad 1: vtok-vy´tok . . . . . . . . . . . . . R ˇ esˇenı´ Stokesovy u´lohy - proudeˇnı´ v kaviteˇ, rychlost (vlevo), tlak (vpravo) R Q2 − Q1 prvek (• dveˇ slozˇky rychlosti, ◦ tlak) . . . . . . . . . . . . . . . . . Q2 − Q1 makroprvek M = k ∪ m . . . . . . . . . . . . . . . . . . . . . . Q2 − Q1 sı´t’(◦ rychlost, × tlak) . . . . . . . . . . . . . . . . . . . . . . . . . ˇ esˇenı´ Stokesovy u´lohy na jemne´ sı´ti . . . . . . . . . . . . . . . . . . . . . . R ˇ esˇenı´ Navierovy-Stokesovy u´lohy na 256 × 256 sı´ti . . . . . . . . . . . . . R

6

7 12 13 22 22 23 25 29

29 35 36 37 42 44 45 46 57

3

1 1.1

´ vod U Trochu historie, fyziky a odvozovańı´

Uzˇ od prvopocˇa´tku, kdy lidstvo zacˇalo uzˇ´ıvat tekutiny ve svu˚j prospeˇch, ne jen k pit’´ı a dyćhańı´, se snazˇ´ı pochopit za´kladnı´ za´konitosti mechaniky tekutin. Starˇ´ı rˇ´ımane´, kterˇ´ı povazˇovali vodu za jeden ze cˇtyrˇ zˇivlu˚, se naucˇili vyuzˇ´ıvat jejı´ za´kladnı´ vlastnost, zˇe tecˇe z mı´st vy´sˇe polozˇenyćh do mı´st nı´zˇe polozˇenyćh, kterou vyuzˇili prˇi stavbeˇ akvaduktu, ˇ ´ıma. A tak bychom mohli da´le pokracˇovat prˇes strˇetehdy nejdu˚lezˇiteˇjsˇ´ı stavby stare´ho R doveˇky´ rozvoj morˇeplavectvı´ azˇ do doby dnesˇnı´, kdy se bez podrobne´ho modelovańı´ tekutin neobejdeme. Burdzˇ Dubaj (Dubajska´ veˇzˇ), most Pont de Millau (Millau Viaduc), BAC Concorde, program Space Shuttle, zde jen ma´lo prˇ´ıkladu˚, ktere´ by bez modelovanı´ proudeˇnı´ nebylo mozˇne´ realizovat. Vsˇe nenı´ jen o technice, ale take´ o osobnostech:

Archime´des (287-212, prˇ.n.l.) - staroveˇky´ rˇecky´ matematik a fyzik, mechanik, vyna´lezce. Je vyna´lezcem sˇroubove´ho cˇerpadla, jezˇ je vyuzˇ´ıvańo azˇ do dnes, hlavneˇ v cˇistı´rnaćh odpadnıćh vod. Zna´my´ je take´ Archime´du˚v za´kon: „Teˇleso ponorˇene´ do kapaliny je nadlehcˇovańo silou, ktera´ se rovna´ tı´ze kapaliny teˇlesem vytlacˇene´“. Leonardo da Vinci (1452-1519) na´vrh mlyńa s vodnı´m kolem, objevenı´ za´kona rychlosti a pru˚rˇezu v · S = konst , na´vrhy vhodnyćh tvaru˚ pro cˇluny, pada´ku, atd.. Galileo Galilei (1564–1642) - toskańsky´ astronom, filosof a fyzik teˇsneˇ spjaty´ s veˇdeckou revolucı´ oveˇrˇil, zˇe tlak v kapalineˇ za´visı´ na „vy´sˇce“ sloupce kapaliny. Blaise Pascal (1623-1662) - francouzsky´ fyzik a matematik a jeho za´kon: „Tlak vyvolany´ vneˇjsˇ´ı silou v kapalineˇ je ve vsˇech smeˇrech a ve vsˇech mı´stech stejny´“. A mnoho dalsˇ´ıch nemeńeˇ vy´znamnyćh „pru˚kopnı´ku˚“ v oblasti mechaniky kontinua. Doporucˇuji prˇecˇ´ıst [12]. Pro na´s budou v tuto chvı´li nejdu˚lezˇiteˇjsˇ´ı dva pańove´. Tı´m prvnı´m je Navier Claude Louis (1785–1836) - francouzsky´ inzˇeny´r a fyzik, ktery´ v roce 1822 navrhl jistou soustavu parcia´lnıćh diferencia´lnıćh rovnic jako model popisujıćı´ visko´znı´ nestlacˇitelne´ tekutiny. Kdyzˇ se tehdejsˇ´ı fyzikove´, matematici podı´vali na jeho odvozenı´ rovnic, okamzˇiteˇ je smetli ze stolu, fyzika´lnı´ prˇedpoklady byly naprosto nerealisticke´. Pozdeˇji, v roce 1845, Gabriel Stokes 1 odvodil mnohem rigoro´zneˇjsˇ´ım zpu˚sobem model stejne´ rovnice jako prˇed 23-mi lety Navier. Dobre´ je si uveˇdomit rozdı´l mezi tekutinami, tedy kapalinami a plyny a teˇlesy pevne´ho skupenstvı´. Rozdı´l tkvı´ ve znacˇne´ pohyblivosti cˇa´stic. Nemajı´ vlastnı´ tvar a jsou snadno deˇlitelne´. Z toho take´ plyne, zˇe koheznı´ sı´ly jsou uvnitrˇ tekutin mnohem mensˇ´ı nezˇ u pevnyćh teˇles, tekutiny kladou relativneˇ maly´ odpor sila´m pu˚sobıćı´m ve smeˇru vneˇjsˇ´ı norma´ly plochy, ktera´ tekutinu omezuje. Z tohoto du˚vodu je proto vhodneˇjsˇ´ı mluvit u tekutin o tlaku mı´sto o napeˇtı´. 1 (1819-1903), jeden z nejvy´znamneˇjsˇ´ıch irskyćh veˇdcu˚. Je povazˇovań za jednoho ze zakladatelu˚ hydrodynamiky. Urcˇil metodu pro meˇrˇenı´ kinematicke´ viskozity (Stokesu˚v viskozimetr). Uka´zal na prostou linea´rnı´ za´vislost napeˇtı´ na deformacˇnı´ rychlosti, cˇ´ımzˇ zobecnil Newtonu˚v za´kon.

4

Odpor tekutiny proti zmeˇneˇ tvaru se nazy´va´ viskozita (vazkost). Vlastnost, jezˇ se projevuje vzˇdy, kdyzˇ tekutina nenı´ v rovnova´ze a prˇizpu˚sobuje se vneˇjsˇ´ım sila´m na tekutinu pu˚sobıćı´m. Sı´la viskozity ma´ snahu zeslabit rozdı´l rychlostı´ (vza´jemnyćh) v proudıćı´ tekutineˇ, cˇ´ımzˇ prˇipomıńa´ sı´ly vznikajıćı´ prˇi trˇenı´, proto taka´ neˇkdy viskozitu nazy´va´me vnitrˇnı´m trˇenı´m 2 . Vıće informacı´ o tekutinaćh v [20] o rozlisˇenı´ hranice mezi plyny a kapalinami v [7]. V na´sledujıćı´ rˇa´dcıćh si nebudeme rigoro´zneˇ odvozovat Navier-Stokesovu rovnici, obecneˇ lze nale´zt v [20], [7] matematicke´ za´lezˇitosti (funkciona´lnı´ analy´za) lze najı´t v [21], [24], [8], [26], ale pokusı´me se uka´zat za´kladnı´ vlastnosti mechaniky tekutin, ktere´ bud’ prˇ´ımo popisuje Navier-Stokesova rovnice nebo je prosteˇ vhodne´ si uveˇdomit. Budeme vycha´zet z [12], odkud cˇerpa´me obra´zky. Prˇi rˇesˇenı´ u´loh v mechanice tekutin se vycha´zı´ z prˇedstavy tekutiny jako spojite´ho, stejnorode´ho prostrˇedı´. Stejnorodostı´ neboli izotropiı´ rozumı´me stejne´ vlastnosti vsˇech cˇa´stecˇek kapaliny neza´visle´ na jejich poloze a smeˇru pu˚sobenı´ sil. Da´le je vhodne´ zave´st elementa´rnı´ objem tekutiny jedna´ se o objem velmi maly´ vzhledem k rozmeˇru na´doby cˇi proudu, avsˇak dostatecˇneˇ velky´ vzhledem ke strˇednı´ volne´ dra´ze molekuly budeme jej znacˇit dV = dxdydz jako objem hranolu. Sı´ly, ktere´ mohou pu˚sobit na elementa´rnı´ objem: • sı´ly plosˇne´ • sı´ly hmotnostnı´ (objemove´),ktere´ lze deˇlit: – vneˇjsˇ´ı hmotnostnı´ sı´ly (naprˇ. tı´hovou) – sı´ly od pohybu kapalin (naprˇ. odstrˇedivou, setrvacˇnou) kde jsou definovańy jako vektory o trˇech karte´zskyćh sourˇadnicıćh. Da´le vı´me, zˇe proudeˇnı´ skutecˇne´ (vazke´) tekutiny mu˚zˇe by´t dvojı´ho druhu: • lamina´rnı´ proudeˇnı´, kdy se cˇa´stice tekutiny se pohybujı´ v tenkyćh vrstvaćh, anizˇ se prˇemı´st’ujı´ po pru˚rˇezu zkoumane´ oblasti, • turbulentnı´ proudeˇnı´, kdy se cˇa´stice tekutiny majı´ kromeˇ pode´lne´ rychlosti take´ turbulentnı´ (fluktuacˇnı´) rychlost, jı´zˇ se prˇemı´st’ujı´ po pru˚rˇezu. Proudeˇnı´ mu˚zˇe by´t take´ za´visle´ na usporˇa´danı´ v prostoru - v jednorozmeˇrne´m Euklidovske´m prostoru v = (x1 ) atd. Rozlisˇujeme rovneˇzˇ proudeˇnı´ podle za´vislosti na cˇase: ∂ = 0, • usta´lene´ (staciona´rnı´), ktere´ je neza´visle´ na cˇase v 6= v(t); ∂t

• neusta´lene´ (nestaciona´rnı´), ktere´ je za´visle´ na cˇase v = (x1 , t). 2 V dynamice se cˇasto lze setkat s pohybem kdy se viskozita neuvazˇuje nebo je prˇ´ılisˇ mala´, takovou tekutinu nazy´va´me dokonalou tekutinou.

5

Prˇechodu lamina´rnı´ho proudeˇnı´ v turbulentnı´ je pro urcˇite´ potrubı´ a tekutiny dańo kritickou rychlostı´. Z pokusu˚ i teorie podobnosti vyply´va´, zˇe prˇechod lamina´rnı´ho proudeˇnı´ v turbulentnı´ (a naopak) je urcˇen Reynoldsovy´m kriticky´m cˇ´ıslem, ktere´ je definovańo vztahem: v·d Re = , (1.1) ν kde v je strˇednı´ rychlost proudıćı´ tekutiny, d je charakteristicky´ rozmeˇr (naprˇ. prˇi proudeˇnı´ v potrubı´ jeho pru˚meˇr) a ν je kinematicka´ viskozita proudıćı´ tekutiny. Prvnı´ rovnicı´, kterou zacˇneme bude Eulerova rovnice hydrostatiky, ktera´ vyjadrˇuje obecnou podmıńku rovnova´hy sil hmotnostnıćh (F0 ) (resp. vneˇjsˇ´ı objemove´ sı´ly (naprˇ. gravitacˇnı´ nebo odstrˇediva´ sı´la)) a tlakovyćh (Fp ) pu˚sobıćıćh na pro jednotku hmotnosti 1 kg tekutiny v klidu. Euler rˇ´ıkal, zˇe vektorovy´ soucˇet hmotnostnı´ a tlakove´ sı´ly je roven nule F0 + Fp = 0.

(1.2)

Eulerovy rovnice lze prˇepsat do diferencia´lnı´m tvaru ∂p ∂p ∂p aρ − gradp = 0, kde gradp = , , . ∂x ∂y ∂z

(1.3)

Rozsˇ´ırˇenı´m je Eulerova rovnice pro proudeˇnı´ idea´lnı´ tekutiny, ktera´ je rozsˇ´ırˇenı´m Eulerovy rovnice hydrostatiky o setrvacˇne´ sı´ly od vlastnı´ho pohybu cˇa´stic dokonale´ tekutiny. F0 + Fp = Fs .

(1.4)

Opeˇt lze prˇepsat rovnici do diferencia´lnı´ho tvaru3 ∂v 1 Dv = + (v · grad) v = a − gradp. Dt ∂t ρ

(1.5)

Pozna´mka 1.1 

   vx vx ∂ ∂ ∂  vy  = (v · grad) v = (v · grad)  vy  = vx + vy + vz ∂x ∂y ∂z vx vx 

x vx ∂v ∂x

 y =  vx ∂v ∂x z vx ∂v ∂x

x + vy ∂v ∂y

∂v

+ vy ∂yy z + vy ∂v ∂y

x + vz ∂v ∂z



∂v  + vz ∂zy  z + vz ∂v ∂z

Poslednı´m rozsˇ´ırˇenı´m Eulerovy rovnice pro proudeˇnı´ idea´lnı´ tekutiny o trˇecı´ sı´lu dosta´va´me Navier-Stokesovu rovnici pro nestlacˇitelnou tekutinu, ktera´ vyjadrˇuje rovnost setrvacˇne´ sı´ly vu˚cˇi soucˇtu hmotnostnı´, tlakove´ a trˇecı´ sı´ly F0 + Fp + Ft = Fs .

(1.6)

6

Obra´zek 1: Pu˚sobenı´ sil

7

Obra´zek 2: Linea´rnı´ rozlozˇenı´ rychlosti v prˇi lamina´rnı´m proudeˇnı´ mezi dveˇma deskami (vlevo) Rychlostnı´ profil a tecˇne´ napeˇtı´ (vpravo)

Pozna´mka 1.2 Pro matematicke´ vyja´drˇenı´ trˇecıćh sil se vyuzˇ´ıva´ Newtonu˚v za´kon, jehozˇ rigoro´znı´ odvozenı´ vyzˇaduje hlubsˇ´ı znalosti. Pro prˇedstavu uvazˇujme na´sledujıćı´ prˇ´ıklad proudeˇnı´ skutecˇne´ tekutiny mezi dveˇma rovinny´mi deskami. Hornı´ deska se pohybuje rychlostı´ U a druha´ je v klidu, cozˇ take´ znamena´, zˇe na pohybujıćı´ se desce ma´ cˇa´stice rychlost U , kdezˇto na desce v klidu je rychlost cˇa´stice nulova´. Pro ostatnı´ cˇa´stice v mezerˇe mezi deskami jsou rychlosti rozlozˇene´ linea´rneˇ. Pohybujıćı´ se cˇa´stice rychlostı´ strha´va´ sousednı´ cˇa´stice do pohybu na´sledkem visko´znı´ho trˇenı´. Rychlost cˇa´stice ve vzda´lenosti y od desky v klidu bude y v=u . h Smykove´ smykove´ napeˇtı´ od vazkosti je podle Newtona vyja´drˇeno vztahem τ =η

u dv =η , dy h

(1.7)

kde τ je tecˇne´ napeˇtı´, η je dynamicka´ viskozita. Vy´raz vyjadrˇuje vztah mezi tecˇny´m napeˇtı´m a derivacı´ rychlosti podle sourˇadnice kolme´ na smeˇr pohybu. Vra´tı´me-li se zpeˇt k trˇecı´ sı´le, mu˚zˇeme ji vyja´drˇit obdobneˇ jako sı´lu tlakovou Ft = τ S = η

dv dv S = νρ S, dy dy

(1.8)

kde ν je kinematicka´ viskozita. Stanovı´-li se rovnova´ha vsˇech sil pu˚sobıćıćh na elementa´rnı´ objem tı´m, zˇe se drˇ´ıve odvozena´ Eulerova rovnice rozsˇ´ırˇ´ı o zmeˇny norma´lovyćh a tecˇnyćh 3

nazy´va´me substanciona´lnı´ nebo te´zˇ individua´lnı´ derivaci funkce rychlosti podle Kde tota´lnı´ derivaci Dv Dt cˇasu. Levy´ cˇlen v prostrˇednı´ cˇa´sti rovnosti (1.5) nazy´va´me loka´lnı´ derivacı´ a druhy´ cˇlen se nazy´va´ konvektivnı´ (proudova´) derivace funkce podle cˇasu.

8

napeˇtı´, prˇesneˇji rˇecˇeno o jejich derivace podle sourˇadnic, dostane se Navierova - Stokesova rovnice, ktera´ ve vektorove´m za´pise pro nestlacˇitelnou tekutinu v pravou´hle´m sourˇadne´m syste´mu ma´ tvar Dv ∂v 1 = + (v · grad) v = a − gradp + ν∆v, (1.9) Dt ∂t ρ 2 ∂ ∂2 ∂2 kde ∆ = ∇ · ∇ = ∂x + + je Laplaceu˚v opera´tor, v nasˇem prˇ´ıpadeˇ aplikovany´ na 2 ∂y 2 ∂z 2 trˇi slozˇky rychlosti. Tudı´zˇ poslednı´ cˇlen v (1.9) prˇedstavuje sı´lu potrˇebnou k prˇekonańı´ visko´znı´ho trˇenı´ tekutiny. Navier-Stokesou rovnici mu˚zˇeme rozepsat do trˇ´ı smeˇru˚ sourˇadnic x, y, z 2 ∂ v x ∂ 2 v x ∂ 2 vx ∂vx ∂vx ∂vx 1 ∂p ∂vx + vx + vy + vz = ax − +ν + + (1.10) ∂t ∂x ∂y ∂z ρ ∂x ∂x2 ∂y 2 ∂z 2 2 ∂vy ∂ vy ∂vy ∂vy ∂vy ∂ 2 vy 1 ∂p ∂ 2 vx + vx + vy + vz = ay − +ν + + (1.11) ∂t ∂x ∂y ∂z ρ ∂y ∂x2 ∂y 2 ∂z 2 2 ∂ vz ∂vz ∂vz ∂vz 1 ∂p ∂ 2 vz ∂ 2 vz ∂vz + vx + vy + vz = az − +ν + + (1.12) ∂t ∂x ∂y ∂z ρ ∂z ∂x2 ∂y 2 ∂z 2 Dosta´va´me syste´m parcia´lnıćh diferencia´lnıćh rovnic na prvnı´ pohled usporˇa´danyćh s jedinou nelinearitou. Tato jedna nelinearita deˇla´ syste´m o neˇkolik rˇa´du˚ na´rocˇneˇjsˇ´ım. Z hlediska charakterizace nelinearit je pro N = 2 proble´m kriticky´, ale zvla´dnutelny´, zatı´mco pro N = 3 jde o proble´m dosud nezvla´dnuty´. Ne nadarmo nadace Clayova matematicke´ho institutu zarˇadila vyrˇesˇenı´ Navierovy-Stokesovy rovnice na seznam sedmi nejdu˚lezˇiteˇjsˇ´ıch matematickyćh proble´mu, kdy na kazˇdy´ z nich je vypsańa odmeˇna milion dolaru˚[10]. Prvnı´m, kdo se pokusil Navierovy-Stokesovy rovnice serio´zneˇ vyrˇesˇit, byl C.W.Ossen. Na jeho praći nava´zal J. Lerray, ktery´ doka´zal existenci rˇesˇenı´ pro N = 2, ale pro N = 3 neuspeˇl. Doka´zal existenci tzv. „turbulentnı´ho rˇesˇenı´“ (veˇrˇil, zˇe turbulence je zodpoveˇdna´, za prˇ´ıpadne´ singularity). Navrhl take´ mozˇnost, jak uka´zat, zˇe klasicke´ rˇesˇenı´ nemusı´ obecneˇ existovat. Nakonec se uka´zalo, zˇe jeho na´vrh je chybny´, takzˇe NADEˇJE PRO ˇ ESENI´ ZˇIJE !!! NALEZENI´ KLASICKE´HO R

9

´ vod do metody konecˇnyćh prvku˚ U

2

V na´sledujıćıćh podkapitolaćh si uka´zˇeme velmi hruby´ u´vod do variacˇnı´ formulace, metody konecˇnyćh prvku˚ a metody smı´sˇenyćh prvku˚. Budeme vycha´zet [23], [3], [5] a [2], kde lze najı´t podrobneˇjsˇ´ı informace.

2.1

Jednorozmeˇrna´ u´loha

Uvazˇujme na´sledujıćı´ diferencia´lnı´ rovnici sˇ´ırˇenı´ tepla v tenke´m (homogennı´m a izotropnı´m) dra´tu prostrˇednictvı´m kondukce s reaktivnı´m cˇlenem s jednotkovou tepelnou vodivostı´ a s koeficientem µ prˇestupu tepla do vneˇjsˇ´ıho prostrˇedı´ na uzavrˇene´m intervalu rea´lnyćh cˇ´ısel 0 < x < 1, kde u je teplota (formulaci lze najı´t v [2]) Lu ≡ −

du d κ(x) + µ(x)u = f (x), dx dx

(2.1)

hranicˇnı´ podmıńky: u(0) = 0

Dirichletova hranicˇnı´ podmıńka4 ,

(2.2)

du (1) = 0 Neumannova hranicˇnı´ podmıńka.5 (2.3) dx Klıćˇem k metodeˇ konecˇnyćh prvku˚ je zajisˇteˇnı´ ortogonality rezidua (Lu − f ) a vsˇech prvku˚ v z vhodne´ho vektorove´ho prostoru (Lu − f, v) = 0,

∀v ∈ V 6 ,

(2.4)

kde vy´razem (., .) je mysˇlen L2 skala´rnı´ soucˇin prˇes integracˇnı´ oblast Z (u, v) = uv dx Ω = h0, 1i . Ω

Je nutno upravit rovnici (2.4), protozˇe obsahuje druhe´ derivace, jezˇ jsou slozˇiteˇjsˇ´ı na aproximaci v ra´mci metody konecˇnyćh prvku˚. Pomocı´ integrace (2.4) po cˇa´stech vy´jde: (κ

du  du  du dv   , ) + (µu, v) = (f, v) + κ v  − κ v  . dx dx dx 1 dx 0

(2.5)

Pouzˇitı´m Neumannovy hranicˇnı´ podmıńky (2.3) v bodeˇ x = 1 prvnı´ hranicˇnı´ cˇlen vypadne. Druhy´ cˇlen (v x = 0) vypadne, pokud v splnı´ homogennı´ hranicˇnı´ podmıńku jako u v(0) = 0. Cˇasto take´ „esencia´lnı´“ hranicˇnı´ podmıńka. Cˇasto take´ „prˇirozena´“ hranicˇnı´ podmıńka. 6 Myslı´me „pro vsˇechna v lezˇ´ıcı´ ve V.“ 4 5

10

Tı´m dosta´va´me na´sledujıćı´ formulaci pu˚vodnı´ho proble´mu (2.1): hleda´me u takove´ u(0) = 0, ∀v takove´, zˇe v(0) = 0,

a(u, v) = (f, v),

(2.6)

prˇ´ıslusˇna´ bilinea´rnı´ forma: Z

a(u, v) =

0

1

κ

du dv ) + µuv dx. dx dx

(2.7)

Prˇedchozı´ forma je symetricka´ pro jake´koliv u a v. Mu˚zˇeme tedy usoudit, zˇe v bude z podobne´ „trˇ´ıdy“ prvku˚ jako u, jizˇ si pokusı´me definovat. Je trˇeba si uveˇdomneˇt, zˇe proble´m (2.6) mu˚zˇeme formulovat jako optimalizacˇnı´ proble´m (energetickou formulaci), hleda´me u takove´, zˇe u(0) = 0 a J(u) = min J(v), v|v(0)=0

s prˇ´ıslusˇny´m funkciona´lem

1 J(v) = a(v, v) − (f, v). 2 Podı´va´me se na prostor, prˇes ktery´ minimalizujeme funkciona´l J(v). Integra´l obsazˇeny´ v (2.6) a (2.7) musı´ by´t konecˇny´ a prvky u a v musı´ by´t z Hilbertova prostoru, pak jsme schopni nale´zt minimum funkciona´lu J(v). To na´s vede k na´sledujıćı´ definici mnozˇiny. V = mnozˇina takovyćh v pro neˇzˇ platı´ Z 1 |v|2 dx < +∞, 0

Z

0

dv 2 dx < +∞, dx

1

(2.8)

v(0) = 0.

Prˇedpokla´dejme, zˇe κ(x) je veˇtsˇ´ı nezˇ 0 κ(x) ≥ κ0 > 07 . Vsˇimneˇme si, zˇe definice (2.8) prostoru V obsahuje Dirichletovu podmıńku (v(0) = 0), ale neobsahuje Neumannovu podmıńku. Prvnı´ hranicˇnı´ podmıńka ma´ vy´znam v definici prostoru V , zatı´mco Neumannova nema´, du˚vody najdeme v [23]. Dosta´va´me na´sledujıćı´ formulaci proble´mu, hleda´me u ∈ V takove´, zˇe: a(u, v) = (f, v) 7

Pak

p

∀v ∈ V.

a(v, v) definuje normu na V , ekvivalentnı´ klasicke´ Sobolevoveˇ normeˇ:

(2.9) r

R1 0

2 dv |v|2 + dx dx

11

nebo ekvivalentnı´ definice prˇes funkciona´l 1 min J(v) = a(v, v) − (f, v). v∈V 2

(2.10)

Proble´m (2.9) se cˇasto nazy´va´ „variacˇnı´ formulace“ diferencia´lnı´ rovnice (2.1). Budeme se zaby´vat numericky´m rˇesˇenı´m Navier-Stokesovy rovnice, radeˇji zvolı´me termıń „slaba´ formulace“. Metoda konecˇnyćh prvku˚ vycha´zejıćı´ ze slabe´ formulace potrˇebuje ke sve´ konstrukci konecˇneˇ dimenziona´lnı´ podprostor prostoru V : Vh ⊂ V. Pak se zmeˇnı´ formulace (2.9) na formulaci u´lohy hleda´me uh ∈ V takove´, zˇe a(uh , vh ) = (f, vh )

∀vh ∈ Vh .

(2.11)

Popı´sˇeme si ve strucˇnosti nejjednodusˇsˇ´ı po cˇa´stech linea´rnı´ prvek. Nejprve si zavedeme deˇlenı´ intervalu 0 ≤ x ≤ 1: 0 = x0 < x1 < x2 < . . . xN = 1. Funkce vh ∈ Vh musı´ splnˇovat na´sledujıćı´ podmıńky: • vh je spojita´ na (0, 1), • vh je po cˇa´stech linea´rnı´ na intervalech, • vh (0) = 0. Je snadne´ uka´zat, zˇe Vh generuje na´sledujıćı´ mnozˇinu ba´zovyćh 8 funkcı´ viz (Obra´zek 3). x ≤ xi−1 : φi (x) = 0, x − xi−1 , xi−1 ≤ x ≤ xi : φi (x) = xi − xi−1 xi+1 − x , xi ≤ x ≤ xi+1 : φi (x) = xi+1 − xi xi+1 ≤ x : φi (x) = 0. ˇ esˇenı´ hleda´me ve tvaru: Odtud je zrˇejme´, zˇe Vh je dimenze N . R uh (x) =

N X

ui φi (x),

i=1

kde ui je hodnota uh v bodeˇ xi = ux (xi ). Z diskre´tnıćh rovnic (2.11) hleda´me takove´ 8

cˇasty´ nazy´vanyćh „tvarovyćh.“

12

Obra´zek 3: Ba´zova´ funkce (φi (xj ) = δij ) u1 , . . . , uN ktere´ splnˇuje: N X

ui a(φi , φj ) = (f, φj )

∀j = 1, . . . , N,

(2.12)

i=1

(dosta´va´me syste´m N linea´rnıćh rovnic o N nezna´myćh). Pozna´mka 2.1 Matice syste´mu linea´rnıćh rovnic je tridiagona´lnı´, symetricka´ a pozitivneˇ definitnı´. Syste´m, tedy mu˚zˇe by´t rˇesˇen Gaussovou eliminacı´ (pocˇet operacı´ pro vy´pocˇet ≈ O(N )) nebo jinou metodou pro hledańı´ rˇesˇenı´ rˇ´ıdkyćh matic. Pro prˇ´ıpad nehomogennıćh okrajovyćh podmıńek (mı´sto (2.2) ma´me u(0) = u0 ) lze proble´m snadno prˇeve´st na proble´m s homogennı´mi okrajovy´mi podmıńkami vıće informacı´ v [23]. Prˇedpokla´dejme Fourierovu okrajovou podmıńku v x = 1 u(0) = 0, du (1) = q1 , (α > 0, β > 0). dx Z (2.5) dosta´va´me odpovı´dajıćı´ slabou formulaci: Z 1 Z 1 du dv κ(1)α κ(1)α (κ , + µu, v)dx + u(1)v(1) = q1 v(1) f, vdx + dx dx β β 0 0 αu(1) + β

2.2

(2.13)

∀v ∈ V. (2.14)

Dvourozmeˇrna´ u´loha

Bude uvazˇovańa 2-dimenziona´lnı´ analogie rovnice (2.1) −(κij u,j ),j = f

v Ω,

(2.15)

13

Obra´zek 4: Integracˇnı´ oblast, Γ0 Dirichletova hranicˇnı´ podmıńka a Γ1 Neumannova hranicˇnı´ podmıńka u=0

na Γ0 ,

κij u,j nj = 0

Dirichletova hranicˇnı´ podmıńka,

na Γ1 ,

Neumannova hranicˇnı´ podmıńka,

(2.16) (2.17)

Kde - Ω = integracˇnı´ oblast, - Γ = Γ0 ∪ Γ1 = hranice oblasti Ω, viz obra´zek 4, - n = n(x) = jednotkovy´ vektor, norma´la k Ω v bodeˇ xi , smeˇrˇujıćı´ ven z oblasti, se slozˇkami n1 , n2 , - u,j =

∂u ∂xi ,

- κij (x) = tenzor napeˇtı´, splnˇujıćı´ podmıńku κ0 > 0. Musı´ platit podmıńka stejnomeˇrne´ elipticity, jinak viz [23]: κij (x)ξi ξj ≥ κ0 (ξ12 ξ22 )

(2.18)

pro vsˇechna x z Ω, a pro vsˇechny vektory (ξ1 , ξ2 ). Take´ se da´ prˇedpokla´dat, zˇe κij je ohranicˇeny´. Pro dvou rozmeˇrny´ prˇ´ıpad je κij u,j : 2 X i=1

κij u,j = κ1j u,1 + κ2j u,2

14 Prˇı´klad 2.1 κ11 = κ22 = 1, κ12 = κ21 = 0. Z te´to podmıńky se vytvorˇ´ı Poissonova rovnice se dveˇma typy hranicˇnıćh podmıńek: −∆u = f

v Ω,

u = 0 ∂u = 0 ∂n

na Γ0 , na Γ1 .

Homogennı´ Dirichletova podmıńka (2.16) definuje v≡0

na Γ0

(2.19)

Vyna´sobenı´m (2.15) funkcı´ v, na´slednou aplikacı´ integra´lu prˇes Ω s pouzˇitı´m Greenovy veˇty 9 ve vıće rozmeˇrech dosta´va´me Z Z κij u,i v,j dx = f v dx (2.20) Ω

Ω

(Hranicˇnı´ integra´ly jsou identicky rovny 0 dı´ky (2.17) a (2.19)). Potrˇebujeme najı´t takova´ v ∈ V , pro nezˇ by meˇl dany´ integra´l smysl. V = mnozˇina takovy´ v pro nezˇ platı´ Z |v|2 dx < +∞, Z Ω |v,j |2 dx < +∞, Ω v Γ0 = 0.

Pozna´mka 2.2 Jestlizˇe Γ0 bodova´ (ma´ nulovou mı´ru), pak tato poslednı´ podmıńka v (2.21) nema´ smysl a proble´m prˇejde na u´lohu s jedinou Neumannovou hranicˇnı´ podmıńkou, pak u je definovańo azˇ na aditivnı´ konstantu. Slaba´ formulace na´by´va´ tvaru u∈V a(u, v) = (f, v)

∀v ∈ V,

(2.21)

kde a(u, v) = (f, v) =

9

Dvou rozmeˇrna´ integrace po cˇa´stech:

R

Ω

Z

Z

κij u,i v,j dx, Ω

f v dx.

(2.22)

Ω

(φ,i ψ + φψ,i ) dx =

R

Γ

φψni dΓ.

15

Pozna´mka 2.3 Pokud tenzor κij = κji , pak a(u, v) = a(v, u) a (2.21) je ekvivalentnı´ hledańı´ minima energeticke´ho funkciona´lu 12 a(v, v) − (f, v) ve V . Stejnomeˇrna´ elipticita v (2.18), ktera´ zajisˇt’uje existenci a jednoznacˇnost rˇesˇenı´ (viz Lax-Milgaramova veˇta) [5], ukazuje Z Z κij v,j v,j dx ≥ κo v,i v,i dx. Ω

Ω

V symetricke´m prˇ´ıpadeˇ, κij = κji , forma a(v, v)1/2 mu˚zˇe by´t uzˇita, jako norma na V , ktera´ je ekvivalentnı´ obvykle´ Sobolevoveˇ normeˇ.

2.3

Metoda konecˇnyćh prvku˚

Za´kladnı´ za´konitosti metody konecˇnyćh prvku˚10 budou uka´zańy na slabe´ formulaci (2.21), zavedeme uh 11 (diskre´tnı´ rˇesˇenı´): uh =

N X

up φp (x),

(2.23)

p=1

kde φp jsou ba´zove´ funkce a ui jsou stupneˇ volnosti uh . Pak dosta´va´me soustavu rovnic popisujıćı´ diskre´tnı´ rˇesˇenı´: N X

up a(φp , φq ) = (f, φq )

∀q(q = 1, . . . , N ),

(2.24)

p=1

Vh je vektorovy´ prostor generovany´ ba´zovy´mi funkcemi φp (p = 1, . . . , N ), Vh ∈ V . Odtud lze pozorovat, zˇe metoda konecˇnyćh prvku˚ zcela za´visı´ na vy´beˇru prostoru Vh , jiny´mi slovy na vy´beˇru ba´zovyćh funkcı´ φp . Oblast Ω se rozdeˇlı´ na male´ kousıćˇky (prvky, podoblasti, elementy), v jednorozmeˇrne´m Euklidovske´m prostoru na u´secˇky, ve dvourozmeˇrne´m na troju´helnı´ky nebo cˇtyrˇu´helnı´ky a ve trojrozmeˇrne´m na cˇtyrˇsteˇny nebo kva´dry. Ba´zove´ funkce jsou definovańy na kazˇde´ podoblasti. V na´sledujıćı´m prˇehledu budou uka´zańy mnohocˇleny stupneˇ k = 1 nebo 2. Pk = mnozˇina polynomu˚ stupneˇ ≤ k. Odtud pro troju´helnı´ky: • P0 = 1, • P1 = 1, x1 , x2 , • P2 = 1, x1 , x2 , x21 , x1 x2 , x22 , atd.. Prˇi pouzˇitı´ cˇtyrˇu´helnı´ku˚ budou mnozˇiny polynomu˚ = Qk stupneˇ ≤ k: • Q0 = 1, 10 11

Da´le v textu se take´ mu˚zˇe objevit zkratka MKP. „h“ znacˇ´ı rˇa´d sı´teˇ.

16

• Q1 = 1, x1 , x2 , x1 x2 , • Q2 = 1, x1 , x2 , x21 , x1 x2 , x22 , x21 x2 , x1 x22 atd.. Budou uvazˇovańy konecˇne´ elementy Lagrangerova typu, u kteryćh φp je definovańa hodnotami ve specia´lneˇ vybranyćh uzlovyćh bodech. Vy´beˇr uzlovyćh bodu˚ je vedeny´ pozˇadavky na spojitost, aby bylo splneˇno pro kazˇdou φp (a uh ) ∈ V azˇ na nekonformnı´ elementy, viz pozna´mka pod textem, vhodne´ podmıńky jsou: • φp po cˇa´stech polynom, • φp je spojita´ nad kazˇdy´m elementem, • φp = 0 na γ0 . Tyto pozˇadavky na spojitost je jednoduche´ oveˇrˇit, jak bude uka´zańo v pozdeˇjsˇ´ıch kapitolaćh prˇi podrobneˇjsˇ´ım rozboru jednotlivyćh druhu˚ elementu˚. Vıće informacı´ o spojitosti a chybeˇ aproximace v [27]. Pozna´mka 2.4 Budeme se zaby´vat jen konformnı´ metodou konecˇnyćh prvku˚, tj. jedna´ se o situaci, kdy nejsou narusˇeny na´sledujıćı´ variacˇnı´ prˇestupky (variational crimes, viz Fix, Strang) [2] : 1. Ωh 6= Ω oblast je aproximovańa pouze sjednocenı´m konecˇnyćh prvku˚, 2. fh 6= f pouzˇ´ıva´me numerickou integraci, tj. pracujeme s fh mı´sto s f , 3. τbh 6= τb, (b τ prˇedstavuje nehomogennı´ Neumanovu podmıńku naprˇ. κij u,j nj = τb na Γ1 ) pouzˇ´ıva´me numerickou integraci, tj. pracujeme s τbh mı´sto s τb, nebo aproximujeme Γ1 ,

4. u bh 6= u b, (b u prˇedstavuje nehomogennı´ Dirichletovu podmıńku naprˇ. u = u b na Γ0 ) aproximujeme u b nebo Γ0 .

2.4

Galerkinova metoda

Zde budou prˇedstaveny za´klady Galerkinovy metody, vıće v [23] a [2]. Uvazˇujme dvourozmeˇrny´ prˇ´ıpad Poissonovy rovnice na oblasti Ω ⊂ R s Lipschitzovskou hranicı´ [5]. −∆u = −∇2 u = f,

(2.25)

kde ∆ je Laplaceu˚v opera´tor. Poissonova rovnice je tedy parcia´lnı´ diferencia´lnı´ rovnice ˇ esˇenı´ rovnice (2.25) bude take´ vyhovovat hranicˇnı´ podmıńce na ∂Ω, elipticke´ho typu. R uvazˇujme podmıńku na´sledujıćı´ho typu: αu + β

∂u =q ∂n

na ∂Ω,

(2.26)

kde ∂u/∂n znacˇ´ı smeˇrovou derivaci ve smeˇru vneˇjsˇ´ı norma´ly k hranici ∂Ω, α a β jsou rea´lne´ neza´porne´ konstanty. Pokud je β = 0, pak se jedna´ o Dirichletovu hranicˇnı´ podmıńku,

17

v prˇ´ıpadeˇ α = 0 se bude jednat o Neumannovu hranicˇnı´ podmıńku, jinak se jedna´ o Newtonovu hranicˇnı´ podmıńku. Pro vıće informacı´ doporucˇuji prˇecˇ´ıst [23], kapitolu 3, zde si uvedeme jen velmi vzda´leny´ pohled bez vesˇkeryćh podrobnostı´. Pozna´mka 2.5 Nezbytnou podmıńkou pro existenci rˇesˇenı´ Neumannovy u´lohy je podmıńka kompatibility: Z Z g+ f = 0. (2.27) ∂Ω

2.4.1

Ω

Slaba´ formulace

Vezmeme si „rozumnou (testovacı´)“ funkci v : Ω → R a v = 0 na ∂Ω a vyna´sobı´me jı´ Poissonovu rovnicı´ Z ∇2 u + f v = 0. (2.28) Ω

Integra´l je jisteˇ dobrˇe definovany´, ma´ smysl. Pokud zna´me klasicke´ rˇesˇenı´ u, pak zcela jisteˇ musı´ splnˇovat (2.28). Pokud je v dostatecˇneˇ hladka´ funkce pak (2.28), lze upravit pomocı´ Greenovy veˇty. Z Z Z ∂u 2 − v∇ u = ∇u · ∇v − v , ∂Ω ∂n Z Ω ZΩ Z ∂u ∇u · ∇v = (2.29) vf + v . Ω Ω ∂Ω ∂n U poslednı´ rovnice se rˇesˇenı´ u nazy´va´ slabe´ rˇesˇenı´ (vycha´zı´ ze slabe´ formulace u´lohy). Naprˇ. pro Dirichletovu hranicˇnı´ podmıńku (α = 1 a β = 0 v (2.26) budeme hledat slabe´ rˇesˇenı´ u ve tvaru: Z Z Z ∇u · ∇v = vf + vg. Ω

hleda´me tudı´zˇ u takove´, zˇe:

Ω

−∇2 u = f

u = gD na ∂ΩD a

∂Ω

v Ω,

∂u = gN na ∂ΩN . ∂v

(2.30) (2.31)

12 kde ∂ΩD ∪ ∂ΩN = ∂Ω a ∂Ω R D ∩ ∂ΩN ve vsˇech bodech hranice je skoro vsˇude pra´zdny´. Da´le prˇedpokla´dejme, zˇe ∂ΩD ds 6= 0, jiny´mi slovy v (2.31) nebude jen Neumannova podmıńka. Na´sledneˇ definujme prostory pro rˇesˇenı´ a testovacı´ funkce:

HE1 := {u ∈ H 1 (Ω)d |u = gD na ∂ΩD },

(2.32)

HE1 0 := {v ∈ H 1 (Ω)d |u = 0 na ∂ΩD }.

(2.33)

Rozdı´l mezi obeˇma prostory se skry´va´ v Dirichletoveˇ podmıńce v (2.31), jezˇ je prˇ´ımo „zabudovańa“ do definice prostoru HE1 , v neˇmzˇ hleda´me rˇesˇenı´. Oproti tomu funkce 12

Mysˇleno na mnozˇineˇ nulove´ mı´ry.

18

z testovacı´ho prostoru HE1 0 jsou nulove´ na Dirichletoveˇ cˇa´sti hranice. V kontrastu k tomu Neumanovy podmıńky nejsou ani v jednom z prostoru˚ prˇedepsańy na hranici. Z (2.29) je zrˇejme´, zˇe kazˇda´ funkce u, jezˇ splnˇuje (2.30) a (2.31), je take´ rˇesˇenı´m na´sledujıćı´ slabe´ formulace: hleda´me u ∈ HE1 takove´, zˇe: Z Z Z ∇u · ∇v = vf + vgN ∀v ∈ HE1 0 . (2.34) Ω

Ω

∂ΩN

Du˚lezˇity´m poznatkem je, zˇe klasicke´ rˇesˇenı´ Poissonovy u´lohy pozˇaduje aby funkce u byly dvakra´t spojiteˇ diferencovatelne´ v Ω cozˇ je „tvrdsˇ´ı“ pozˇadavek nezˇ konecˇnost integra´lu cˇtverce prvnı´ derivace. Pouzˇitı´ (2.34) opravnˇuje aproximovat rˇesˇenı´ za podmıńek splneˇnı´ urcˇite´ hladkosti a Dirichletovyćh podmıńek zakotvenyćh v (2.32). 2.4.2

Galerkinova metoda konecˇnyćh prvku˚

V te´to cˇa´sti bude odvozena mysˇlenka aproximace u na konecˇneˇ rozmeˇrne´m podprostoru z prostoru HE1 . Budeme vycha´zet ze slabe´ formulace (2.32)-(2.34) z obecne´ho proble´mu (2.30)-(2.31). Da´le prˇedpokla´dejme, zˇe S0h ⊂ HE1 0 je konecˇneˇ n-rozmeˇrny´ vektorovy´ prostor testovacıćh funkcı´, pro ktere´ je {φ1 , φ2 , . . . , φn } vhodna´ ba´ze. Pak je nutne´ zajistit Dirichletovu podmıńku v (2.32) rozsˇ´ırˇenı´ ba´zovyćh funkcı´ o dalsˇ´ı funkce φn+1 ,P φn+2 , . . . , φn+n∂ ∂ a vybrat takove´ pevne´ koeficienty uj , j = n + 1, . . . , n + n∂ , zˇe funkce n+n j=n+1 uj φj inh terpolujı´ hranicˇnı´ funkci gD na ∂ΩD . Na´sledneˇ je konecˇneˇ rozmeˇrna´ aproximace uh ∈ SE T jednoznacˇneˇ prˇirˇazena s vektorem u = (u1 , u2 , . . . , un ) rea´lnyćh koeficientu˚ vy´razem: uh =

n X j=1

uj φj +

n+n X∂

uj φj .

(2.35)

j=n+1

Funkce φj , j = 1, . . . , n v prvnı´ sumeˇ (2.35) definujı´ mnozˇinu tvarovyćh funkcı´. Konstrukce (2.35) velmi zjednodusˇuje charakterizaci diskre´tnı´ho rˇesˇenı´ prˇi pozˇadavku splneˇnı´ Dirichletove´ podmıńky na hranici. h je zajis Konstrukce prostoru SE ˇ teˇna specifickou volbou tvarovyćh funkcı´ v (2.35), jezˇ se shoduje s volbou testovacıćh funkcı´, jezˇ tvorˇ´ı ba´zi prostoru S0h . Obecneˇ se tato metoda jmenuje Galerkinova (prˇesneˇji Bubnov-Galerkinova) aproximacˇnı´ metoda. Vıće obecna´ je metoda, prˇi ktere´ se pouzˇ´ıvajı´ rozdı´lne´ aproximacˇnı´ prostory (2.32)-(2.33), tj. rozdı´lne´ tvarove´ a testovacı´ funkce. Tato alternativa je zna´ma´ jako Petrov-Galerkinova aproximacˇnı´ metoda. Vy´sledek Galerkinovy aproximace je konecˇneˇ dimenziona´lnı´ verze slabe´ formulace: hleda´me uh ∈ S0h takove´ , zˇe Z Z Z 1 ∇uh · ∇vh = vh f + vh g N ∀vh ∈ SE . (2.36) 0 Ω

Ω

∂ΩN

19

Na´sledneˇ je vhodne´ prosadit si (2.36) pro kazˇdou ba´zovou funkci a na´sledny´m dosazenı´m (2.35) vyply´va´, zˇe (2.36) je ekvivalentnı´ hledańı´ uj , j = 1, . . . , n takovyćh, zˇe n X j=1

uj

Z

Z

∇φj · ∇φi = Ω

φi f + Ω

Z

φi gN − ∂ΩN

n+n X∂

Z

uj

j=n+1

∇φj · ∇φi

(2.37)

Ω

pro kazˇde´ i = 1, . . . , n. Prˇedchozı´ formulace mu˚zˇe by´t zapsańa ve formeˇ matice jako syste´m linea´rnıćh rovnic Au = f, (2.38) kde A = [aij ], a f = [fi ],

fi =

Z

aij =

φi f + Ω

Z

Z

∇φj · ∇φi ,

(2.39)

Ω

φ i gN − ∂ΩN

n+n X∂

j=n+1

uj

Z

∇φj · ∇φi .

(2.40)

Ω

Syste´m linea´rnıćh rovnic (2.38) se nazy´va´ Galerkinu˚v syste´m a funkce uh vypocˇ´ıtańa po dosazenı´ rˇesˇenı´ z (2.38) do (2.35) se nazy´va´ Galerkinovo rˇesˇenı´. Matice A se cˇasto nazy´va´ matice tuhosti. Galerkinova matice koeficientu˚ (2.39) je zrˇejmeˇ symetricka´ a take´ je pozitivneˇ definitnı ´ zat. Uvazˇujeme obecny´ vektor v odpovı´dajıćı´ specificke´ funkci vh = Pn ´. Cozˇ lze uka h ˇe j=1 vj φj ∈ S0 tak, z vT Av =

n X n X

vj aij vi =

j=1 i=1

=

n X n X

vj

j=1 i=1

= =

Z

Z

Ω

 

n X j=1

Z

∇φj · ∇φi vi = Ω

 !

vj ∇φj  ·

n X i=1

vi ∇φi

"

=

∇vh · ∇vh = Ω

≥ 0. Odtud lze videˇt, zˇe cˇtvercova´ a symetricka´ matice A je prˇi nejmensˇ´ım semidefinitnı´. vT Av = 0 jen v prˇ´ıpadeˇ, kdy je ∇vh = 0. To nastane pouze je-li vh je konstantnı´ v Ω. Da´le vı´me, zˇe vh ∈ S0h je spojita´ azˇ do hranice a je nulova´ na ∂ΩD , tudı´zˇ pokud ∇vh = 0 musı´ by´t i vh = 0. Testovacı´ funkce tvorˇ´ı ba´zi S0h , jestlizˇe vh = 0, pak v = 0. Neumannu˚v proble´m (2.34) vyzˇaduje hlubsˇ´ı u´vahu, zde budou uka´zańy jen hlavnı´ mysˇlenky. Galerkinova matice je pouze semidefinitnı´, v tomto prˇ´ıpadeˇ a ja´dro vektoru v odpovı´da´ funkci ∇vh = 0. Pro tuto situaci je nutne´ oveˇrˇit, zˇe podprostor S h ⊂ H 1 dany´

20

mnozˇinou tvarovyćh funkcı´ {φj }, j = 1, P . . . , n definuje rozdeˇlenı´ jednotky, tj. kazˇdy´ vektor h s S musı´ by´t prˇirˇazen s vektorem vh = nj=1 vj φj splnˇujıćı´m n X

φj = 1.

(2.41)

j=1

Konstrukce (2.41) rˇ´ıka´, zˇe jestli je vh konstantnı´ funkce, vh ≡ α, pak vh je prˇirˇazeno diskre´tnı´mu vektoru, jezˇ splnˇuje vi = α pro vsˇechny koeficienty. To take´ znamena´, zˇe ja´dro Galerkinovy matice spojene´ s (2.34) ma´ dimenzi rovnou jedne´, sesta´vajıćı´ z vektoru konstantnıćh koeficientu˚. Da´le poznamenejme, zˇe rˇesˇitelnost Diskre´tnı´ho Neumannova syste´mu (analogie (2.38)) vyzˇaduje, aby ja´dro Galerkinovy matice A bylo ortogona´lnı´ k vektoru prave´ strany f, tj. (1, . . . , 1)T f = 0, kde Z Z f = [fi ], fi = φ i f + φi g. (2.42) Ω

∂ΩN

Pozna´mka 2.6 Prˇ´ıpad, kdy α = 0, β = 0 v (2.26), mu˚zˇe by´t u konstantnı´ splnˇujıćı´ homogennı´ proble´m s f = 0, g = 0, pak Poissonova u´loha ma´ rˇesˇenı´ lisˇ´ıcı´ se v aditivnı´ konstantneˇ. Nebo integracı´ (2.25) prˇes Ω uzˇit’´ım Gaussovy veˇty (viz Appendix), dosta´va´me Z Z Z ∂u 2 = − ∇ u = f, (2.43) ∂Ω ∂n Ω Ω cozˇ na´m da´va´ nezbytnou podmıńku pro existenci rˇesˇenı´ Poissonovy u´lohy s Neumannou podmıńkou Z Z g+

∂Ω

f = 0.

(2.44)

Ω

Uzˇitı´ vlastnosti (2.41) ukazuje, zˇe diskre´tnı´ Poissonova u´loha s Neumannou podmıńkou je rˇesˇitelna´ jen pokud tento hranicˇnı´ proble´m je dobrˇe definovań, pokud platı´ (2.44). h a Sh Vrat’me se zpeˇt k obecne´mu prˇ´ıpadu (2.36), chteˇli jsme uka´zat, zˇe vy´beˇr z SE 0 ma´ za´sadnı´ vy´znam, jestli ma´ uh vztah ke slabe´mu rˇesˇenı´ u. Obecny´ pozˇadavek je vybrat h a S h , aby dostatec ˇ na´ aproximace rˇesˇenı´ byla dosazˇena pro dostatecˇneˇ velka´ n, tj. SE 0 pozˇadujeme, aby chyba skutecˇne´ho a aproximovane´ho rˇesˇenı´ ku − uh k se zmensˇovala pro rostoucı´ n. Motivacı´ pro uzˇitı´ MKP je, zˇe hladke´ funkce (naprˇ. rˇesˇenı´) mu˚zˇe by´t aproximovańo po cˇa´stech polynomy (stupneˇ Pn , n ∈ N). Pokud uvazˇujeme Galerkinu˚v syste´m (2.36), mysˇlenkou je vybrat ba´zove´ funkce {φj } v (2.35), jezˇ jsou loka´lneˇ nenulove´ na sı´ti troju´helnı´ku˚ nebo obde´lnı´ku˚ (R2 ), jehlanu˚ nebo cˇtyrˇsteˇnu˚ (R3 ). V te´to praći se budeme zaby´vat vy´hradneˇ cˇtyrˇu´helnı´kovou sı´tı´. 2.4.3

ˇ tyrˇu´helnı´kove´ prvky C

Pro jednoduchost prˇedpokla´dejme, zˇe Ω ⊂ R2 je oblast s Lipschitzovskou hranicı´ [5] a jsme schopni danou oblast rozdeˇlit na mnozˇinu cˇtyrˇu´helnı´ku˚ k , k = 1, . . . , K. Definujme obecne´ deˇlenı´ Th oblasti Ω na cˇtyrˇu´helnı´ky. To znamena´, zˇe hrany sousednıćh obde´lnı´ku˚ sply´vajı´ a zˇe

21

•

S

k

k = Ω,

• p ∩ h = ∅ pro p 6= h, Body, kde se hrany cˇtyrˇu´helnı´ku˚ setka´vajı´ budeme nazy´vat uzly, budou oznacˇeny indexy j = 1, . . . , n. Uzly, lezˇ´ıcı´ na hranici oblasti Ω nazveme hranicˇnı´ uzly a ostatnı´ nazveme vnitrˇnı´mi uzly. Zvolı´me-li neˇktery´ uzel j a sestrojı´me-li jehlan („pyramidu“) s vy´sˇkou rovnou jedne´ v uzlu j a jehozˇ podstava je n-u´helnı´k tvorˇeny´ sousednı´mi uzly uzlu j (v prˇ´ıpadeˇ nasˇ´ı sı´teˇ je podstavou cˇtyrˇu´helnı´k), pak tento jehlan definujeme ba´zovou funkci φj tak, zˇe pro body podstavy jehlanu vytvorˇene´ho v tomto bodeˇ je nenulova´ a vsˇude jinde v Ω nulova´. Jednou z mozˇnostı´ je vzı´t po cˇa´stech linea´rnı´ funkci (o kvadraticke´ funkci pozdeˇji) na kazˇde´m cˇtyrˇu´helnı´ku˚, jezˇ bude rovna jedne´ v uzlu j a nula ve vsˇech ostatnıćh bodech sı´teˇ. Ba´zova´ funkce je zrˇejmeˇ spojita´ na Ω a je dostatecˇneˇ hladka´, tj. φj ∈ H 1 , dı´ky tomu prostor S0h = span(φ1 , φ2 , . . . , φn ) z (2.36). Cˇtyrˇu´helnı´kove´ prvky jsou z podstaty meńeˇ flexibilnı´ nezˇ troju´helnı´kove´, je le´pe uvazˇovat deˇlenı´ oblasti na obde´lnı´ky cˇi cˇtverce, nezˇ obecne´ cˇtyrˇu´helnı´ky, zejmeńa pro oblasti s „linea´rnı´ hranicı´“ (obde´lnı´kove´, L-oblasti a jejich odnozˇe), cˇasto je vhodne´ volit kombinaci s troju´helnı´kovy´mi prvky, ktere´ je lepsˇ´ı vyuzˇ´ıt pro aproximaci zakrˇivenyćh hranic. Nejjednodusˇsˇ´ı cˇtyrˇu´helnı´kovy´ element je Q1 . Pro jednoduchost budeme uvazˇovat pravou´hly´ cˇtyrˇu´helnı´k. Kazˇda´ funkce mu˚zˇe by´t zrˇejmeˇ vyja´drˇena ve tvaru (ax+b)(cy +d). Kazˇdou ze cˇtyrˇ ba´zovyćh funkcı´ φj na elementu definujeme standarteˇ tak, zˇe je rovna jedne´ v uzlu j a nula v ostatnıćh uzlech. Naprˇ. na elementu x ∈ h0, hi, y ∈ h0, hi pak ba´zove´ funkce na elementu, v protismeˇru hodinovyćh rucˇicˇek jsou (1 − x/h)(1 − y/h),

x/h(1 − y/h),

xy/h2 ,

y/h(1 − x/h)

Globa´lnı´ ba´zova´ funkce na podstaveˇ tvorˇene´ cˇtyrˇmi elementy je zobrazena na (5). V prˇ´ıpadeˇ obecne´ho prvku (nepravidelne´ho, s rozdı´lny´mi vnitrˇnı´mi u´hly) se mu˚zˇe libovolneˇ meˇnit poloha bodu˚ a nenı´-li mezi nimi pravidelnost zvolı´me si izoparametrickou transformaci podle nı´zˇ vsˇechny elementy „sjednotı´me“. Hlavnı´ mysˇlenka je v definovańı´ funkcı´: χ1 (ξ, η) = (ξ − 1)(η − 1)/4, χ2 (ξ, η) = −(ξ + 1)(η − 1)/4, χ3 (ξ, η) = (ξ + 1)(η + 1)/4, χ4 (ξ, η) = −(ξ − 1)(η + 1)/4,

(2.45)

na referencˇnı´m prvku ∗ , ξ ∈ h−1, 1i, η ∈ h−1, 1i. Pak mu˚zˇeme prˇeve´st obecny´ cˇtyrˇsteˇn na referencˇnı´ element. To bude provedeno tak, zˇe sourˇadnice vrcholu˚ obecne´ho cˇtyrˇsteˇnu (xv , yv ), v = 1, 2, 3, 4 prˇevedeme na referencˇnı´ element pomocı´ transformace sourˇadnic

22

Obra´zek 5: Q1 ba´zova´ funkce

Obra´zek 6: Transformace sourˇadnic Q1 prvku (6): x(ξ, η) = y(ξ, η) =

4 X

v=1 4 X

xv χv (ξ, η), yv χv (ξ, η).

(2.46)

v=1

Cˇtverka v sumeˇ vyjadrˇuje pocˇet uzlu˚ pro jednoznacˇne´ urcˇenı´ linea´rnı´ho elementu, pro prˇ´ıpad kvadraticke´ho elementu je trˇeba 9 uzlu˚ a pro prˇ´ıpad kubicke´ho elementu˚ 16 uzlu˚, (7). Jak bylo pouka´zańo mu˚zˇeme definovat kvadraticky´ element, prˇidańı´m dalsˇ´ıch uzlu˚ do strˇedu stran a jednoho do strˇedu Q1 elementu. V tomto prˇ´ıpadeˇ budou tvorˇit ba´zi cˇtyrˇi funkce pro vrcholy, cˇtyrˇi pro hrany a jedna ve strˇedu. Aproximaci budeme hledat ve tvaru (ax2 + bx + c)(dy 2 + ey + f ) na kazˇde´m pravidelne´m cˇtyrˇu´helnı´ku, cozˇ je kombinace sˇesti

23

Obra´zek 7: Druhy cˇtyrˇsteˇnovyćh elementu˚ kvadratickyćh funkcı´ 1, x, y, xy, x2 , y 2 , dvou kubickyćh x2 y, y 2 x a jedne´ bikvadraticke´ x2 x2 .Ba´zova´ funkce vybudovańa nad Q2 prvkem je zrˇejmeˇ spojita´. 2.4.4 Sestavenı´ matic pro cˇtyrˇrozmeˇrne´ prvky Budeme zde uvazˇovat tzv. matici tuhosti Ak , jezˇ se vypocˇ´ıta´va´, jak ukazuje (7), z referencˇnı´ho prvku ∗ na dane´m prvku k a na´sledny´m uzˇitı´m urcˇite´ho kvadraturnı´ho vzorce. Pro cˇtyrˇsteˇny (pro jednoduchost uvazˇujme nejdrˇ´ıve Q1 prvek) je transformace definovańa pro vsˇechny body (x, y) ∈ k dana´ vzorci (2.45) a (2.46). Transformace z referencˇnı´ho prvku na k je diferencovatelna´. Tudı´zˇ diferencovatelnou funkci ϕ(ξ, η), mu˚zˇeme transformovat jako # " #" # " ∂ϕ ∂ξ ∂ϕ ∂η

=

∂x ∂ξ ∂x ∂η

∂y ∂ξ ∂y ∂η

∂ϕ ∂x ∂ϕ ∂y

.

Jacobiho matice v (2.47) je tvorˇena linea´rnı´mi funkcemi o sourˇadnicıćh (ξ, η) " P # P4 ∂χj ∂χj 4 y x ∂(x, y) j j ∂ξ ∂ξ Jk = . Pj=1 Pv=1 ∂χj ∂χj 4 4 ∂(ξ, η) x y j j j=1 v=1 ∂η ∂η

(2.47)

(2.48)

Loka´lnı´ matice tuhosti bude vypocˇ´ıtańa na´sledujıćı´m vztahem, odvozeny´m v [13] m X m X ∂ψ∗,i ∂ψ∗,j ∂ψ∗,i ∂ψ∗,j (k) + , wst |Jk (ξ2 , ηt )| aij = ∂x ∂x ∂y ∂y s=1 t=1

kde (ξ2 , ηt ) jsou body pro vy´pocˇet integra´lu pomocı´ Gaussovy kvadratury.

2.5

Smı´sˇena´ formulace a smı´sˇena´ metoda konecˇnyćh prvku˚ pro elipticke´ rovnice

Ve strucˇnosti rˇesˇeno smı´sˇena´ znamena´, zˇe gradient z rˇesˇenı´ je cha´pań jako pomocna´ neza´visla´ promeˇnna´, tj. mı´sı´me (s urcˇity´mi omezenı´mi) dva typy aproximacı´ pro u a pro gradient u. Smı´sˇena´ metoda konecˇnyćh prvku˚ se velmi cˇasto uzˇ´ıva´ prˇi aproximaci rˇesˇenı´ Stokesovy a Navier-Stokesovy rovnice (pro formulaci rychlost-tlak nebo proudova´ funkce-vı´rˇivost). Jedina´ cena, kterou platı´me, je znacˇny´ na´ru˚st pocˇtu rovnic. V te´to kapitole bude vysveˇtlena podstata smı´sˇene´ metody konecˇnyćh prvku˚ a nezbytnyćh podmıńek pro jejı´ konvergenci. Vıće v [14] a [27].

24

2.5.1

Jednorozmeˇrny´ proble´m

Pro lepsˇ´ı prˇedstavu o smı´sˇene´ metodeˇ konecˇnyćh prvku˚ odvodı´me postupy v jednorozmeˇrne´m „umeˇle´m“ prˇ´ıpadeˇ z knihy [27]. Uvazˇujme: d du − κ = f v h0, 1i , dx dx u(0) = 0, (2.49) du (1) = 0. dx (2.49) mu˚zˇeme take´ prˇepsat ve formeˇ: dp = f v h0, 1i , dx du = 0 v h0, 1i , p−κ dx u(0) = 0, −

(2.50)

p(1) = 0. Slaba´ formulace odvozena´ od (2.50): Z 1 dp + f v dx = 0, dx 0 Z 1 dq −1 κ pq + u dx = 0, dx 0

∀v ∈ L2 (Ω), ∀q ∈ X,

(2.51)

kde X je mnozˇina takovyćh funkcı´ q na uzavrˇene´m intervalu 0 < x < 1, pro ktere´ platı´: Z Z

0

1

|q|2 dx < +∞,

0

1

dq 2 dx < +∞, dx q(1) = 0.

(p ∈ X to naznacˇuje Neumannovu hranicˇnı´ podmıńku v (2.50)(p(1) = 0).) Pozna´mka 2.7 Druha´ rovnice ve (2.51) prˇedstavuje slabou formulaci p − κ(du/dx) = 0 s hranicˇnı´ podmıńkou u(0) = 0. Da´le pak nenalezneme zˇa´dne´ derivace u a v ve (2.51), dı´ky tomu mu˚zˇeme uvazˇovat nespojitou aproximaci u. Pro popis smı´sˇene´ metody diferencia´lnı´ rovnice nizˇsˇ´ıho rˇa´du potrˇebujeme: • po cˇa´stech konstantnı´ aproximaci u, • po cˇa´stech linea´rnı´ aproximaci q,

25

Obra´zek 8: Ba´zove´ funkce Mh (hornı´ obra´zek), Xh (spodnı´ obra´zek) • deˇlenı´ intervalu h0, 1i, na N intervalu˚ (u´secˇky v jednorozmeˇrne´ u´loze). Zavedeme si mnozˇinu Mh = mnozˇina po cˇa´stech konstantnıćh funkcı´. Mnozˇina ba´zovyćh funkcı´ z Mh je dańa funkcı´ wi : 1 na h0, 1i , wi = 0 jinde . A mnozˇinu Xh = mnozˇina po cˇa´stech linea´rnıćh a spojityćh funkcı´. Nulova´ v bodeˇ x = 1. Mnozˇina ba´zovyćh funkcı´ z Xh je dańa „strˇ´ısˇkovou“ funkcı´ φi = (viz Obra´zek 8) Pozna´mka 2.8 Xh ⊂ X, Mh ⊂ L2 (0, 1). Pak z (2.51) dosta´va´me rovnice popisujıćı´ diskre´tnı´ rˇesˇenı´: Z 1 dph + f v dx = 0, ∀v ∈ Mh , dx 0 Z 1 dq −1 κ p h q + uh dx = 0, ∀q ∈ Xh , dx 0 (poslednı´ rovnice definuje diskre´tnı´ derivaci uh ), ekvivalentneˇ: Z xi+1 dph + f dx = 0, i = 0, 1, . . . , N − 1, dx xi Z 1 dqi κ−1 ph qi + uh dx = 0, i = 0, 1, . . . , N − 1. dx 0

(2.52)

(2.53)

26

Prˇedpokla´dejme, zˇe κ ≡ κ0 > 0, tı´mto prˇedpokladem nezmeˇnı´me strukturu syste´mu rovnic, jedna´ se cˇisteˇ jen o zjednodusˇenı´. Stanovı´me: xi+1 − xi , 3κ0 = ph (xi ),

λi+1/2 = pi

ui+1/2 = uh na hxi , xi+1 i, Z xi+1 fi+1/2 = (xi+1 − xi ) fdx. xi

(2.54)

Pak rovnice popisujıćı´ diskre´tnı´ rˇesˇenı´ mohou by´t explicitneˇ napsańy: pi+1 − pi = fi+1/2 , pN

i = 0, 1, . . . , N − 1,

= 0,

1 λ1/2 p0 + λ1/2 p1 − u1/2 = 0, 2

(i = 0),

1 1 λi−1/2 pi−1 + (λi−1/2 + λi+1/2 )pi + λi+1/2 pi+1 + ui−1/2 − ui+1/2 = 0 (0 < i < N − 1) 2 2 1 λ pN −2 + (λN −3/2 + λN −1/2 )pN −1 + uN −3/2 − uN −1/2 = 0 (i = N − 1) 2 N −3/2 Ve formeˇ matice:

kde 

       A=       

A BT B 0

p0 u1/2

1 λ1/2 2 λ1/2 1 1 λ (λ 1/2 + λ3/2 ) 2 1/2 2 λ3/2 .. .. .. . . . .. .. . . 0 1 λ (λ + λi+1/2 ) i−1/2 2 i−1/2 .. . .. .

=

0 f1/2

  B=  

−1

1 −1 1

,

(2.56)



0 ..

.

1 2 λi+1/2

..

.

..

..

.

1 2 λN −3/2

1 2 λN −3/2



.

(λN −3/2 + λN −1/2 )

        (2.57)       



    1  −1

(Skutecˇnost, zˇe rozmeˇry A a B se rovnajı´ je zrˇejma´ z formulace jednorozmeˇrne´ u´lohy). Neˇktere´ skutecˇnosti by meˇly by´t zdu˚razneˇny:

(2.55)

27

• neˇktere´ z diagona´lnıćh elementu˚ v matici (2.56) jsou nulove´ a rˇesˇicˇ potrˇebuje tuto informaci k vy´pocˇtu 13 . • Kdyzˇ vyrˇesˇ´ıme (2.55) pro pi , zı´ska´me: pi = ph (x)i = −

Z

1

f (x)dx, xi

cozˇ je prˇesne´ rˇesˇenı´ v uzlech. Naproti tomu: xi 1 uh (0) = u1/2 = p0 + pi = O(h). 3κ0 2

(2.58)

Odtud vidı´me, zˇe Dirichletova podmıńka je splneˇna pouze s urcˇitou chybou. Vıće v [26]. Pro konvergenci te´to metody (nebo metod vysˇsˇ´ıho rˇa´du) mu˚zˇeme vyslovit na´sledujıćı´ tvrzenı´: ! !1/2 2 Z 1 d |ph − p|2 + (ph − p) + |uh − u|2 dx = O(h). dx 0 Podstata du˚kazu se opı´ra´ o na´sledujıćı´ du˚lezˇite´ tvrzenı´, zna´me´ jako LadysˇenskajaBabusˇka-Brezzi podmıńka: R1

dq v dx dx ∀v ∈ Mh , sup 1/2 ≥ C q∈Xh R 1 dq 2 2 q + dx dx 0 0

Z

1

2

|v| dx 0

2

,

(2.59)

kde C je konstanta neza´visla´ na h. Je velmi du˚lezˇite´ porozumeˇt vy´znamu L.B.B. podmıńce. Jestlizˇe leva´ strana nerovnosti bude vyhovat v = uh , pak 1

dq dx = 0, ∀q ∈ Xh . dx 0 R1 Odtud z druhe´ rovnice ve (2.52) vyply´va´ 0 κ−1 p2h dx = 0 a p ≡ 0, da´le pak (2.59) implikuje uh ≡ 0, tj. implikacı´ ze (2.59) plyne, je-li diskre´tnı´ derivace uh nulova´, pak uh ≡ konst. (= 0 z hranicˇnıćh podmıńek). Jiny´mi slovy (2.59) testuje 14 konzistenci aproximacı´ derivacı´, tj. interpolace pro u a p = κdu/dx nemu˚zˇe by´t vybrańa jedna neza´visle na druhe´.V nasˇem specificke´m prˇ´ıkladeˇ R1 je jednoduche´ oveˇrˇit (2.59) (pro neˇjake´ v z Mh , bude q(x) = 0 v(ξ) dξ, pak q ∈ Xh a dq/dx = v). Mnohem teˇzˇsˇ´ı je po technicke´ strańce oveˇrˇit podmıńku (2.59) v dvou a trˇ´ı rozmeˇrne´m prˇ´ıpadeˇ. Z

13

uh

Proble´m je zrˇejmy´, z prˇ´ıme´ho rˇesˇenı´ (2.55), v prvnı´m kroku vyrˇesˇ´ıme pi a pak dosazenı´m do trˇetı´ azˇ pa´te´ rovnice dosta´va´me ui+1/2 14 Take´ z jine´ho pohledu, (2.59) zajisˇt’uje, zˇe syste´m rovnic ma´ jednoznacˇne´ rˇesˇenı´.

28

2.5.2

Dvourozmeˇrny´ prˇ´ıpad

Uvazˇujme dvourozmeˇrny´ prˇ´ıpad: −(κij u,j ),j

v Ω ⊂ R2 ,

= f

u = 0 κij u,j nj

na Γ0 ,

= 0

na Γ1 .

(2.60)

Prˇ´ıpad mu˚zˇe by´t rozdeˇlen pomocı´ zavedenı´ kogradientu p = (p1 , p2 ) pj,j + f

= 0

v Ω,

pj − κij u,j

= 0

v Ω,

u = 0

na Γ0 ,

pj n j ≡ p · n = 0

na Γ1 .

(2.61)

Druhou rovnici mu˚zˇeme prˇepsat jako: λij pi = u,j , platı´: [λij ] = [κij ]−1 .

(2.62)

Z (2.61) a (2.62) lze odvodit slabou formulaci, kterou mu˚zˇeme v obecne´ formeˇ pro p ∈ X a u ∈ M zapsat na´sledujıćı´m zpu˚sobem: ∀q ∈ X,

a(p, q) + b(q, v) = 0 Z b(p, v) + f v dx = 0

∀v ∈ M,

(2.63)

Ω

kde a(p, q) =

Z

λij pi qj , dx, Z Z b(q, v) = qj,j v dx = div qv dx, Ω

Ω

(2.64)

Ω

Z Z 2 2 X = q(x) : |q| dx < +∞, |div q| dx < +∞, p • n = 0 na Γ1 , Ω

Ω

Z 2 2 M = L (Ω) = v : v dx < +∞, .

(2.65)

Ω

Z (2.65) vyply´va´ vhodnost aproximace u ∈ M nespojity´mi funkcemi, naprˇ. po cˇa´stech polynomy stupneˇ k. Hlavnı´ proble´m se ukry´va´ v konstrukci vhodne´ aproximace pro p ∈ X. Omezı´me se nynı´ na troju´helnı´kovou sı´t’, pak prostor vhodnyćh funkcı´ pro u bude: Mh = {v takove´, zˇe na kazˇde´m troju´helnı´ku T, v = konst. }. Konstrukce takove´to aproximace pro p se opı´ra´ o na´sledujıćı´ tvrzenı´:

29

Obra´zek 9: Zobrazeni norma´lovyćh smeˇru˚ na troju´helnı´ku

Obra´zek 10: Obecny´ troju´helnı´k T (vlevo), mi uzly ve strˇedech stran, ni jednotkove´ norma´love´ vektory smeˇrˇujıćı´ ven z T . Referencˇnı´ troju´helnı´k Tˆ (vpravo), mi uzly ve strˇedech stran, ni jednotkove´ norma´love´ vektory, i = 1, 2, 3. 15 Tvrzenı R ´ 2.1 2Jestlizˇe q je hladke´ na kazˇde´m troju´helnı´ku, pak postacˇujıćı´ a nezbytna´ podmıńka pro Ω |div q| dx < +∞ je spojitost norma´lovyćh slozˇek q na spojnici troju´helnı´ku (tj. podmıńka toku), du˚kaz v [26]:

q|T .nT + q|T ′ .nT ′

= 0

nT + nT ′

= 0.

na ∂T ∩ ∂T ′ , (2.66)

Viz Obra´zek 9. Z prˇedchozı´ho vyply´va´, zˇe stupneˇ volnosti z q ∈ Xh by meˇly by´t norma´love´ slozˇky q v bodech hranice troju´helnı´ku. Je vhodne´ pouzˇ´ıt referencˇnı´ troju´helnı´k Tˆ k definovańı´ ba´zovyćh funkcı´ na T , viz Obra´zek 10. Nejprve definujeme ba´zove´ funkce na Tˆ tak, zˇe φî (m ˆ j ).ˆ njT = δij : 15

q|T ∈ H 1 (T ).

30

φˆ1 (ξ) = φˆ3 (ξ) = ˆ φ3 (ξ) =

ξ1 −1 + ξ2 −2 + ξ2 = φˆ31 ξ2 = φˆ32 −2 + ξ2 = φˆ31 ξ2 = φˆ32

Pak definujeme ba´zove´ funkce na T prˇedpisem: φi (x) =

∂FT ˆ +1 φir , 2. obsah (T ) ∂φr

i = 1, 2, 3.

(2.67)

Prostor funkcı´ Xh je generovań ba´zovy´mi funkcemi φi , jezˇ jsou prˇirˇazeny vsˇem uzlu˚m v troju´helnı´kove´ sı´ti 16 . Pak syste´m rovnic popisujıćı´ diskre´tnı´ rˇesˇenı´ naby´va´ tvaru z (2.63): a(ph , q) + b(q, uh ) = 0 Z b(ph , v) + f v dx = 0

∀q ∈ Xh , ∀v ∈ Mh ⊂ L2 (Ω),

(2.68)

Ω

Mu˚zˇeme nynı´ vyslovit dostatecˇne´ podmıńky (platne´ pro jaky´koliv syste´m formy (2.68)) konvergence: 1. Bilinea´rnı´ forma a(., .) je pozitivneˇ definitnı´ na L2 : Z a(q, q) ≥ C1 |q|2 dx.

(2.69)

Ω

2. L.B.B. podmıńka: ∀v ∈ Mh , sup

q∈Xh

s kqk2X

=

Z

Ω

$

b(q, v) ≥ C2 kqkX

Z

v dx Ω

1/2

,

(2.70)

|q|2 + (div q)2 dx.

Poznamenejme, zˇe konstanty C2 , C1 jsou neza´visle´ na velikosti sı´teˇ. Vy´znam (2.69) a (2.70) je stejny´ jako v jednorozmeˇrne´m prˇ´ıpadeˇ: 1. pro dane´ uh , prvnı´ rovnice v (2.68) jednoznacˇneˇ urcˇuje diskre´tnı´ gradient ph , 2. ph prˇirozeny´ gradient uh , tj. ph = 0 spolecˇneˇ s hranicˇnı´mi podmıńkami ma´ za na´sledek, zˇe uh = 0. Du˚kaz tvrzeni v [14]. 16

Pro kazˇdy´ prvek z q ∈ Xh , se da´ oveˇrˇit, zˇe div q je po cˇa´stech konstantnı´ funkce.

31

Pozna´mka 2.9 1. Jestlizˇe je bilinea´rnı´ forma a(., .) symetricka´, tj. κij = κji , pak forma (2.68) mu˚zˇe by´t prˇevedena na hledańı´ sedlove´ho bodu Lagrangeovy funkce: L(p, v) ≤ L(p, u) ≤ L(q, u), kde

∀q ∈ Xh , v ∈ Mh ,

1 L(q, v) = a(q, q) + b(q, v) − (f, v). 2

2. Hlavnı´ proble´m spocˇ´ıva´ v sestavenı´ sı´teˇ, na nı´zˇ u´lohu rˇesˇ´ıme a take´ v nalezenı´ vhodne´ metody pro rˇesˇenı´ linea´rnı´ho syste´mu.

32

3

Stokesova u´loha

3.1

Prˇedstavenı´ Stokesovy rovnice

Uvazˇujme staciona´rnı´ rovnici (prˇi uzˇitı´ znacˇenı´ z prˇedchozıćh kapitoly) σij,j + fi = 0, σij

= −pδij + ν(ui,j + uj,i ),

ui,i = 0,

(3.1)

kde f = (f1 , f2 ) nebo f = (f1 , f2 , f3 ) - vneˇjsˇ´ı sı´la, u = (u1 , u2 ) nebo u = (u1 , u2 , u3 ) - rychlost, p - tlak, ν - viskozita = konstanta17 , σij - tenzor napeˇtı´. Pouzˇitı´m podmıńky nestlacˇitelnosti div u ≡ ui,i = 0 a prˇedpokladu f = 0, Stokesova rovnice mu˚zˇe by´t napsańa ve tvaru −∆u + grad p = −∇2 u + ∇p = 0, div u = 0.

(3.2)

Vy´sˇe uvedeny´ syste´m prˇedstavuje za´kladnı´ model vazke´ tekutiny. Aplikacı´ Laplaceova opera´toru na vektorovou funkci u dostaneme soucˇet parcia´lnıćh derivacı´ v kazˇde´m z karte´zskyćh smeˇru˚. Prvnı´ rovnice (3.2) ve sve´ podstateˇ prˇedstavuje za´kon zachovanı´ hybnosti (tj. rovnost tlakove´ a trˇecı´ sı´ly) vice v [12] a [9]. Druha´ rovnice reprezentuje podmıńku nestlacˇitelnosti tekutiny (tj. vynucuje za´kon zachovańi hmoty). Za´kladnı´ prˇedpoklad pro modelovańı´ kapalin pomocı´ Stokesovy rovnice je velky´ pomeˇr vizkoznıćh sil oproti advekci vnitrˇnıćh tlakovyćh sil. Z teorie mechaniky tekutin mu˚zˇeme prˇevzı´t pojem Reylnoldsova cˇ´ısla [12] a [20]. Dı´ky tomu mu˚zˇeme tvrdit, zˇe Stokesova rovnice popisuje proudeˇnı´ kapalin s nı´zky´m Reylnoldsovy´m cˇ´ıslem (Re ≪ 1), cozˇ popisuje prˇ´ıpad kdy je rychlost kapaliny velmi mala´ a viskozita velmi vysoka´ 18 . Pouzˇ´ıva´ se pro studium pohybu ru˚znyćh mazadel, lubrikantu˚ a lidske´ krve. Da´le pro popis tekutin s pevny´mi cˇa´stmi, jako jsou cˇa´stecˇky prachu v motorove´m oleji, pigmenty v barveˇ, mikroorganismu a ejakula´tu. V neposlednı´ rˇadeˇ se noveˇ zacˇala pouzˇ´ıvat pro popis tekoucı´ la´vy. Cı´lem te´to kapitoly bude uka´zat si za´kladnı´ souvislosti pro modelovanı´ tekutiny pomocı´ Stokesovyćh rovnic, abychom pozdeˇji mohli prˇejı´t na Navier-Stokerovy rovnice, pro neˇzˇ je Stokesovo proudeˇnı´ limitnı´m prˇ´ıpadem (cˇasteˇji se pouzˇ´ıva´ pojem linearizovany´m 17 18

V dalsˇ´ım textu budeme pro jednoduchost uvazˇovat ν = 1. Odtud se neˇkdy rovnice nazy´va´ rovnice „medu“

33

prˇ´ıpadem).Podkladem pro tuto podkapitolu jsou praće [26], [27], [22], [13], [7] a [14]. Pro obeˇ rovnice Stokesovu i Navier-Stokesovu bude platit podmıńka nestlacˇitelnosti tekutiny. Skutecˇnost, zˇe tlak nenı´ obsazˇen v rovnici kontinuity (druha´ rovnice v (3.2)) znamena´ obtı´zˇ pro modelovańı´ pomocı´ metody konecˇnyćh prvku˚. Jednodusˇe rˇecˇeno, konecˇneˇ prvkove´ prostory jezˇ pouzˇ´ıva´me pro aproximaci rychlosti a tlaku nemohou by´t vybrańy bez za´vislosti jednoho na druhe´m. Abychom splnili urcˇitou „kompabilitu“ mezi prostory, jevı´ se le´pe pouzˇ´ıt smı´sˇenou metodu konecˇnyćh prvku˚. Hlavnı´m cı´lem te´to kapitoly bude analyzovat teorii okolo tvorby vhodnyćh prostoru˚ a prvku˚, jezˇ jsou kompatibilnı´ a budou vhodny´mi kandida´ty pro modelovanı´ Navierovi-Stokesovy rovnice. Formulaci vyja´drˇenou rovnicemi (3.2) uvazˇujeme ve dvou a trˇ´ı rozmeˇrne´ oblasti Ω, spolecˇneˇ s okrajovy´mi podmıńkami na ∂Ω = ∂ΩD ∪ ∂ΩN dane´ u = w na ∂ΩD ,

∂u − np = s na ∂ΩN , ∂n

(3.3)

kde n je norma´lovy´ vektor smeˇrˇujıćı´ ven z oblasti a ∂u/∂n vyjadrˇuje derivaci ve smeˇru R norma´ly.V prˇ´ıpadeˇ, zˇe pro Dirichletovu okrajovou podmıńku platı´ ∂ΩD ds 6= 0, pak ma´me zajisˇteˇnou jednoznacˇnost rˇesˇenı´. Da´le pak Dirichletova hranice mu˚zˇe by´t deˇlena skrze svu˚j vztah s norma´lovou slozˇkou prˇedepsane´ rychlosti ∂Ω+ = {x ∈ ∂Ω|w · n > 0},

vy´tokova´ hranice,

∂Ω0 = {x ∈ ∂Ω|w · n = 0},

charakteristicka´ hranice,

∂Ω− = {x ∈ ∂Ω|w · n < 0},

vtokova´ hranice.

Jestlizˇe je rychlost specifikovańa ve vsˇech bodech hranice, tj. ∂ΩD ≡ ∂Ω, pak rˇesˇenı´ Stokesovy u´lohy (3.2) a (3.3), pro tlak je jednoznacˇneˇ definovańo azˇ na konstantu19 . Navıć integrovańı´m podmıńky nestlacˇitelnosti prˇes Ω a aplikacı´ Gauss-Ostogradske´ho veˇty dosta´va´me Z Z Z w · n. (3.4) u·n= 0= ∇u = Ω

∂Ω

∂Ω

Tudı´zˇ nejednoznacˇnost tlaku je spojena s podmıńkou kompatibility hranicˇnıćh podmıńek Z Z w·n+ w · n = 0. (3.5) ∂Ω−

∂Ω+

Jednodusˇe rˇecˇeno, objem kapaliny, ktery´ vstupuje do oblasti se musı´ rovnat objemu kapaliny, jezˇ z oblasti vycha´zı´. Naproti tomu jinou trˇ´ıdu u´loh prˇedstavuje takzvane´ uzavrˇene´ proudeˇnı´20 , jezˇ je charakterizovańo w · n = 0 ve vsˇech bodech hranice ∂Ω (podmıńka (3.5) je automaticky splneˇna pro tento druh proudeˇnı´). Pokud je na cele´ hranici tecˇna´ slozˇka zrychlenı´ nulova´, pak je tekutina v klidu a nedocha´zı´ k zˇa´dne´mu proudeˇnı´. Pokud alesponˇ na male´ cˇa´sti hranice je tecˇna´ slozˇka zrychlenı´ nenulova´, pak docha´zı´ k pohybu kapaliny neboli k proudeˇnı´. 19 20

V teorii mechaniky tekutin tato konstanta prˇedstavuje tlakove´ hradiny. Typicky´m prˇ´ıkladem je uzavrˇene´ proudeˇnı´ v kaviteˇ.

34

R Na druhou stranu pokud proudeˇnı´ nenı´ uzavrˇene´ (pokud ∂Ω− w · n 6= 0), pak u´loha prˇecha´zı´ na tvar vtok-vy´tok, kde je trˇeba dodrzˇet podmıńku (3.5) nebo (3.2) a (3.3) nebude mı´t rˇesˇenı´. Jednou z mozˇnostı´, jak obejı´t podmıńku (3.5), je zavedenı´ „prˇirozene´ho vy´toku“, cozˇ spocˇ´ıva´ v nahrazenı´ „tvrde´“ Dirichletovy podmı R ńky na vy´toku Neumannovou podmıńkou z (3.3), typicky polozˇenı´m s = 0. Jestlizˇe ∂ΩN ds 6= 0, pak u·n automaticky prˇizpu˚sobı´ vy´tokovou hranici k zajisˇteˇnı´ podmıńky zachovanı´ hmoty ve smyslu (3.4). Navıć rˇesˇenı´ tlaku v (3.2) a (3.3) je vzˇdy jednoznacˇne´.

3.2

´ lohy U

V te´to podkapitole si prˇedstavı´me dveˇ u´lohy, na nichzˇ budeme testovat numerickou implementaci Stokesovy u´lohy. 3.2.1

Prˇ´ıklad 1: vtok-vy´tok

Uvazˇujeme cˇtvercovou oblast Ω = (−1, 1)2 , parabolickou podmıńku na vtoku, prˇirozenou podmıńku na vy´toku a pro ilustraci budeme zna´t analyticke´ rˇesˇenı´. Funkce ux = 1 − y 2 , uy = 0, p = −2x + konst. (3.6) prˇedstavujı´ klasicke´ rˇesˇenı´ Stokesovy u´lohy (3.2), reprezentujıćı´ usta´lene´ horizonta´lnı´ proudeˇnı´ v kana´le pohańeˇne´ tlakovou diferencı´ mezi dveˇma konci. Vsˇeobecneˇ (3.6) je take´ rˇesˇenı´ Navier-Stokesovy rovnice, nazy´vane´ jako Poiseuillovo proudeˇnı´21 . Rˇesˇenı´ Poiseuillova proudeˇnı´ na Ω mu˚zˇe by´t vypocˇ´ıtańo numericky, za pouzˇitı´ rˇesˇenı´ pro rychlost v (3.6), jı´mzˇ definujeme Dirichletovu podmıńku na vtokove´ hranici (x = −1) a charakteristickou hranici (y = −1 a y = 1, pro obeˇ platı´ ux = uy = 0). Vy´tokova´ hranice (x = 1, −1 < y < 1) bude splnˇovat Neumanovu hranicˇnı´ podmıńku (ve smyslu prˇirozene´ho vy´toku): ∂ux − p = 0, ∂x ∂uy = 0. ∂x

(3.7)

Poznamenejme, zˇe Poiseuillovo proudeˇnı´ splnˇuje (3.7) prˇi nulove´m tlaku na vy´toku. Pokud ma´me numericke´ rˇesˇenı´ vypocˇ´ıtane´ pomocı´ MKP pro Stokesovu rovnici (3.2) za´visle´ na (3.3) odpovı´dajıćı´ rˇesˇenı´ zobrazı´ tzv. proudnice22 . 21 Po J. L. M. Poiseuille (1797−1869) - francouzsky´ fyzik.Vyvinul metodu meˇrˇenı´ krevnı´ho tlaku. Zkoumal take´ proudeˇnı´ tekutin trubicemi. Spolecˇneˇ s G.H.L. Hagenem (1797-1884) se zaby´vali lamina´rnı´m proudeˇnı´m. V roce 1838 zjistili, zˇe rychlost proudeˇnı´ za´visı´ na pru˚meˇru a de´lce trubice a na rozdı´lu tlaku mezi konci trubice a formulovali Hagen-Poiseuilleu˚v za´kon. 22 Proudnice (te´zˇ proudova´ cˇa´ra) je trajektorie pohybu jednotlivyćh cˇa´stic prˇi proudeˇnı´ kapalin. Rychlost cˇa´stice v libovolne´m mı´steˇ proudu je tecˇnou k proudnici. Kazˇdy´m bodem proudıćı´ kapaliny procha´zı´ pra´veˇ jedna proudnice. Proudnice se nemohou vza´jemneˇ protıńat. Proudnice lze vyuzˇ´ıt ke graficke´mu zobrazenı´ proudeˇnı´. Jsou-li proudnice rovnobeˇzˇne´, jedna´ se o lamina´rnı´ proudeˇnı´, jsou-li proudnice ru˚znobeˇzˇne´ a ru˚zneˇ stocˇene´, jedna´ se o turbulentnı´ proudeˇnı´.

35

Proudnice: uniformni

ˇ esˇenı´ rychlosti Stokesovy u´lohy, prˇ´ıklad 1: vtok-vy´tok Obra´zek 11: R Definice 3.1 Proudova´ funkce pro dvourozmeˇrne´ Stokesovo proudeˇnı´ ψ : Ω→R, je jednoznacˇneˇ definovańa azˇ na konstantu vztahem ux =

∂ψ , ∂y

uy = −

∂ψ . ∂x

(3.8)

Obecneˇ, okrajove´ podmıńky potrˇebne´ pro proudovou funkci ψ dosta´va´me integracı´ norma´love´ slozˇky rychlosti pode´l hranice. Jako naprˇ´ıklad v tomto prˇ´ıpadeˇ platı´, zˇe norma´lova´ slozˇka rychlosti na charakteristicke´ hranici je nulova´. Integracı´ prˇes vtokovou hranici (x = −1) dosta´va´me Dirichletovu podmıńku ve tvaru: ψ(−1, y) =

Z

y

ux (−1, y) ds = −1

y3 2 +y− . 3 3

Tato podmıńka na´m rˇ´ıka´, zˇe ψ = 0 na spodnı´ hranici a ψ = 4/3 na hornı´ hranici. Norma´lova´ slozˇka z vy´tokove´ podmıńky (3.7) integrovana´ prˇes celou vy´tokovou hranici na´m da´va´ Neumannovu hranicˇnı´ podmıńku Z y ∂ψ (1, y) = p(1, s) ds, ∂n −1 kde p(1, s) je tlak vypocˇ´ıtany´ metodou konecˇnyćh prvku˚ na vy´toku. Z vy´sˇe uvedene´ho je videˇt, zˇe proudova´ funkce na´m da´va´ uzˇitecˇny´ na´stroj pro vykreslenı´ proudnic. Na obra´zku 11 mu˚zˇeme videˇt vykreslenı´ proudnic a pokles tlaku.Jedna´ se o Poiseuillovo proudeˇnı´, jehozˇ rˇesˇenı´ je vypocˇ´ıtańo pomocı´ metody konecˇnyćh prvku˚. Jedna´ se o linea´rnı´ proudeˇnı´, objem kapaliny jezˇ vstupuje do oblasti se nemeˇnı´ prˇi pru˚chodu oblastı´, kapalina jezˇ do oblasti vstupuje pod urcˇity´m tlakem, pak na´sledneˇ vystupuje na vy´toku pod nulovy´m tlakem.

36

tlakove pole

16

15

14

13

12

11 1

0.5

0 1 0.5 −0.5

0 −0.5 −1

−1

Obra´zek 12: Rˇesˇenı´ tlaku Stokesovy u´lohy, prˇ´ıklad 1: vtok-vy´tok 3.2.2

Prˇ´ıklad 2: Uzavrˇene´ proudeˇnı´

Uvazˇujme cˇtvercovou oblast Ω = (−1, 1)2 a uzavrˇene´ proudeˇnı´ uvnitrˇ oblasti23 . Hornı´ strana cˇtvercove´ oblasti se pohybuje urcˇitou rychlostı´ cˇi s urcˇity´m zrychlenı´m vpravo nebo vlevo. Ru˚zny´mi vy´beˇry rychlostı´ hornı´ hranice Ω dosta´va´me ru˚zne´ modely, budeme uvazˇovat model {y = 1; −1 ≤ x ≤ 1|ux } = 1. Numericky vypocˇ´ıtane´ rˇesˇenı´ s Dirichletovou podmıńkou nulove´ rychlosti na ostatnıćh hranicıćh a vykreslenı´ proudnic lze videˇt na obra´zku 13 .

3.3

Slaba´ formulace

Slaba´ formulace Stokesova systemu vycha´zı´ z rovnic (3.2) Z v · (−∇2 u + ∇p) = 0, Ω Z q∇ · u = 0,

(3.9) (3.10)

Ω

pro vsˇechna v a q z vhodneˇ zvolene´ho prostoru testovacıćh funkcı´. Da´le pak dosta´va´me: Z Z Z v · ∇p = − p∇v + ∇ · (pv) = Ω Ω ZΩ Z = − p∇ · v + pn · v, (3.11) Ω

23

∂Ω

Cˇasto se pro tento prˇ´ıpad nazy´va´ uzavrˇene´ proudeˇnı´ v kaviteˇ.

37

tlakove pole Proudnice: uniformni 1 20 0.5 10 0

0

−0.5

−10 1 1 0

−1 −1

−0.5

0

0.5

0

1

−1 −1

Obra´zek 13: Rˇesˇenı´ Stokesovy u´lohy - proudeˇnı´ v kaviteˇ, rychlost (vlevo), tlak (vpravo) −

Z

v · ∇2 u = Ω

=

Z

∇u : ∇v −

ZΩ

∇u : ∇v +

Ω

Z

∇ · (∇u · v) =

ZΩ

(n∇u)v.

(3.12)

∂Ω

Kde ∇u : ∇v oznacˇuje skala´rnı´ soucˇin po slozˇkaćh, naprˇ. ve dvou rozmeˇrech ∇ux · ∇vx + ∇uy · ∇vy . Kombinace (3.9), (3.11) a (3.12) da´va´ Z Z Z Z ∂u 2 − pn · v, (3.13) v · (−∇ u + ∇p) = ∇u : ∇v − p∇v − Ω Ω Ω ∂Ω ∂n pro vsˇechna v z vhodneˇ zvolene´ho prostoru testovacıćh funkcı´.Pozitivnı´ faktem je, zˇe hranicˇnı´ integra´l v (3.13) odpovı´da´ Neumannoveˇ podmıńce v (3.3)24 . A druhy´ integra´l v (3.13) „odpovı´da´“ leve´ straneˇ (3.10). Rovnice (3.13) a (3.10) spolecˇneˇ s hranicˇnı´ podmıńkou (3.3) vedou na prostory, ve kteryćh hleda´me rˇesˇenı´ pro rychlost i vybı´ra´me testovacı´ funkce H1E H1E0

:= {u ∈ H 1 (Ω)d |u = w na ∂ΩD }, 1

d

:= {v ∈ H (Ω) |v = 0 na ∂ΩD },

(3.14) (3.15)

kde d = 2 nebo d = 3 je dimenze prostoru. Skutecˇnost, zˇe v (3.13) nenı´ zˇa´dna´ derivace tlaku, vede k mysˇlence, zˇe L2 (Ω) je vhodny´ prostor pro aproximaci tlaku p. Navıć, vy´beˇr 24

V anglickyćh knihaćh se tato podmıńka nazy´va´ „Do nothing“, jedna´ se prˇirozenou vy´tokovou podmıńku.

38

testovacıćh funkcı´ q pro tlak z prostoru L2 (Ω) zajistı´, zˇe integra´l na prave´ straneˇ rovnice (3.10) je konecˇny´. Nakonec dosta´va´me slabou formulaci Stokesovy u´lohy (3.2) a (3.3) hleda´me u ∈ H1E a p ∈ L2 (Ω) takove´, zˇe Z Z Z s·v ∀v ∈ H1E0 , (3.16) ∇u : ∇v − p∇v = Ω ∂ΩN Z Ω q∇ · u = 0 ∀q ∈ L2 (Ω). (3.17) Ω

Z konstrukce slabe´ formulace vyply´va´, zˇe existuje - li neˇjake´ rˇesˇenı´ (3.2) a (3.3), pak musı´ splnˇovat (3.16) a (3.17). Rovnice (3.16) a (3.17) mohou by´t prˇeformulovańy na u´lohu optimalizace, hledańı´ sedlove´ho bodu funkce Z Z Z 2 sup inf |∇v| − q∇ · v − s · v, (3.18) v∈HE0 q∈L2 (Ω) Ω

Ω

∂ΩN

Vıće v [27]. Proble´m je v proka´zańı´ jednoznacˇne´ho rˇesˇenı´ slabe´ formulace (azˇ na aditivnı´ konstantu v prˇ´ıpadeˇ ∂Ω = ∂ΩD ). Pro tento ućˇel uvazˇujme na chvı´li homogennı´ verzi (3.16) a (3.17): hleda´me u ∈ H1E a p ∈ L2 (Ω) takove´, zˇe Z Z ∇u : ∇v − p∇v = 0 ∀u ∈ H1E0 , (3.19) Ω Ω Z q∇ · u = 0 ∀ ∈ L2 (Ω). (3.20) Ω

Pokusı´me se najı´t prˇ´ıpad kdy,Ru = 0 a p je konstanta. Nejdrˇ´ıve zvolı´me u = v v (3.19) a q = p v (3.20). ZR toho plyne Ω ∇u : ∇u = 0 toto spolecˇneˇ se skutecˇnostı´, zˇe bilinea´rnı´ forma a(u, v) := Ω ∇u : ∇v je pozitivneˇ definitnı´ na prostotou H1E0 (viz definice A.3 v apendixu), znamena ´ , zˇe rˇesˇenı´ u = 0 je potrˇebne´. Na´slednou substitucı´ u = 0 do (3.19) R vyplyne, zˇe Ω p∇v = 0 pro vsˇechna v ∈ H1E0 . Jiny´mi slovy je trˇeba uka´zat: Z p = konst. pokud ∂Ω = ∂ΩD 1 (3.21) p∇v = 0 pro vsˇechna v ∈ HE0 =⇒ p=0 jinak. Ω Mnohem vıće o podmıńkaćh existence rˇesˇenı´ a stability u´lohy v [13], [14] a [27].

3.4

Smı´sˇena´ metoda konecˇnyćh prvku˚

Budeme vycha´zet z poznatku˚, ktere´ jsme si uvedli v prvnı´ kapitole, ktere´ v te´to podkapitole aplikujeme na Stokesovu u´lohu. Nahradı´me prostory H1E0 a L2 (Ω) konecˇneˇ dimenziona´lnı´mi podprostory Xh0 ⊂ H1E0 a M h ⊂ L2 (Ω). Tak, zˇe ze slabe´ formulace (3.16) a (3.17) se stane soustava linea´rnıćh rovnic. Hleda´me uh ∈ XhE a p ∈ M h takove´, zˇe Z Z Z ∇uh : ∇vh − ph ∇vh = (3.22) s · vh ∀vh ∈ Xh0 , Ω Ω ∂ΩN Z qh ∇ · uh = 0 ∀qh ∈ M h . (3.23) Ω

39 ˇ esˇenı´ pro vektor rychlosti budeme hledat ve tvaru R uh =

nu X

~j + uj φ

j=1

nu +n∂ X

~j , uj φ

(3.24)

j=nu +1

~ j ∈ Xh . Fixujeme koeficienty ~ j } je vektor ba´zovyćh funkcı´ (rychlostnıćh) a Pnu uj φ kde {φ 0 j=1 uj : j = nu +1, . . . , nu +n∂ tak, zˇe druhy´ vy´raz v (3.24) prˇedepisuje hodnoty nehomogennı´ Dirichletovy podmıńky na ∂ΩD . Tlak vyja´drˇ´ıme mnozˇinou skala´rnıćh ba´zovyćh funkcı´ ~j } na´sledujıćı´m zpu˚sobem: {ψ np X pk ψk , (3.25) ph = k=1

Jak bylo prˇedeslańo formulace (3.22) a (3.23) lze vyja´drˇit jako soustavu linea´rnıćh rovnic A BT u f = . (3.26) B 0 p g

Jedna´ se o blokovou matici, kde jednotlive´ bloky jsou dańy: Z ~ i : ∇φ ~j , A = [aij ], aij = ∇φ

(3.27)

Ω

bkj = −

B = [bkj ],

Z

~j , ψk ∇ · φ

pro i, j = 1, . . . , nu a k = 1, . . . , np . Vektory na prave´ straneˇ jsou vyja´drˇeny: Z Z nu +n∂ X ~i − ~ i : ∇φ ~j , s·φ f = [fi ], fi = − u j ∇φ ∂ΩN

g = [gk ],

gk =

j=nu +1

nu +n∂ X

j=nu +1

(3.28)

Ω

uj

Z

(3.29)

Ω

~j , ψk ∇ · φ

(3.30)

Ω

ˇ esˇenı´ uh a ph zı´skane´ dosazenı´m vektoru˚ u ∈ Rnu a ph ∈ Rnp do (3.24) a (3.25) je rˇesˇenı´ R k neˇmuzˇ potrˇebujeme metodu smı´sˇenyćh prvku˚. Jak bylo zmıńeˇno d prˇedstavuje rozmeˇr Euklidovske´ho prostoru a s tı´m souvisejıćı´ pocˇet slozˇek vektoru rychlosti. Od nyneˇjsˇka budeme uvazˇovat, zˇe d = 2. Da´le pak uvazˇujme standardnı´ prostor ba´zovyćh funkcı´ {φj }nj=1 jezˇ vznikajı´ prˇi metodeˇ konecˇnyćh prvku˚, zvolı´me nu = d · n = 2n a definujme mnozˇinu ba´zovyćh funkci pro rychlost: n o ~1, . . . , φ ~ 2n } := (φ1 , 0)T , . . . , (φn , 0)T , (0, φ1 )T , . . . , (0, φn )T {φ (3.31) Da´ se uka´zat indukcı´, zˇe toto rozdeˇlenı´ po slozˇkaćh je prˇirozene´ rozdeˇlenı´ bloku˚ GalerkiT novi soustavy (3.26). Pro u := |ux |1 , . . . , |ux |n , |uy |1 , . . . , |uy |n , na´sledneˇ pak lze (3.32) prˇepsat do tvaru:      A 0 BxT ux fx  0 A ByT   uy  =  fy  . (3.32) p g B x By 0

40

kde matice A je Laplaceova matice rˇa´du n×n a Bx , By reprezentujı´ slabe´ derivace v smeˇru x, y a jsou rˇa´du np × n: Z A = [aij ],

∇φi · ∇φj ,

aij =

Ω

bx,ki = −

Bx = [bx,ki ],

by,kj = −

By = [by,kj ],

3.5

Z

ψk

∂φi , ∂x

(3.34)

ψk

∂φj , ∂x

(3.35)

Ω

Z

Ω

(3.33)

Podmıńky rˇesˇitelnosti

V te´to cˇa´sti bude uka´zańo, za jakyćh podmıńek je diskre´tnı´ Stokesova u´loha (3.26)-(3.30) rˇesˇitelna´ a jestli odpovı´da´ spojite´ u´loze (3.2),(3.3). Uvedeme si za´kladnı´ poznatky, hlubsˇ´ım studiem se zaby´vajı´ knihy [26], [15], [14] a [13], z nichzˇ jsme cˇerpali. Prvnı´ bude zodpoveˇzena druha´ cˇa´st ota´zky. Podı´va´me se na dimenzi matice (3.26), je zrˇejme´, zˇe je-li np > nu , pak ma´me nedostatecˇnou hodnost matice, prˇi nejmensˇ´ım np − nu , tj. vysoka´ dimenze pro aproximaci tlaku v porovnańı´ k rychlosti. Z Frobeniovy veˇty, vyply´va´, zˇe soustava by meˇla nekonecˇneˇ mnoho rˇesˇenı´. Podobneˇ jak byla rˇesˇena jednoznacˇnost rˇesˇeni v prˇedcha´zejıćı´ kapitole. I zde si uka´zˇeme kdy je matice (3.26) jednoznacˇneˇ rˇesˇitelna´. Nejdrˇ´ıve bude rozebrań homogennı´ syste´m: Au + B T p = 0,

(3.36)

Bu = 0.

(3.37)

Nejprve si prˇena´sobı´me prvnı´ rovnici cˇlenem uT a druhou rovnici cˇlenem pT . Jednoduchy´mi u´pravami dosta´va´me uT Au = 0 nebo prˇi pouzˇitı´ (3.32) dosta´va´me uTx Aux + uTy Auy = 0. Ale A je pozitivneˇ definitnı´ matice, tedy u = 0 cozˇ implikuje jednoznacˇnou rˇesˇitelnost se zrˇetelem na rychlost. Na druhe´ straneˇ, jednoznacˇna´ rˇesˇitelnost s ohledem na tlak je vıće problematicka´. Substitucı´ u = 0 do (3.37) dosta´va´me B T p = 0, ale to take´ rˇ´ıka´, zˇe rˇesˇenı´ tlaku je jednoznacˇneˇ urcˇeno ja´drem matice B T . Pointa je, zˇe pokud (3.26) rˇa´dneˇ reprezentuje spojitou Stokesovu u´lohu, pak prostory, ve kteryćh hleda´me smı´sˇenou aproximaci, musı´ by´t vybrańy zvla´sˇteˇ opatrneˇ. Specia´lneˇ, musı´me zajistit, zˇe null(B T ) = {1} pro prˇ´ıpad uzavrˇene´ho proudeˇnı´ a null(B T ) = {0}, jinak. Podmıńky na matici B T souvisı´ s podmıńkou (3.21) z prˇedchozı´ kapitoly. K vysveˇtlı´ te´to analogie analyzujeme na´rocˇneˇjsˇ´ı proble´m bez prˇirozene´ vy´tokove´ hranice ∂Ω = ∂ΩD a vyuzˇijeme skutecˇnosti, zˇe ja´dro matice B T je izomorfnı´ k omezene´mu prostoru Z  Nh := qh ∈ M h  qh ∇ · vh = 0, ∀vh ∈ Xh0 . (3.38) Ω

Dany´ prostor lze prˇepsat do vektorove´ formy: 

o n  N h := q ∈ Rm  B T q, v = 0, ∀vh ∈ Rnu .

(3.39)

41 Je du˚lezˇite´ si uveˇdomit, zˇe ph ∈ Nh pra´veˇ kdyzˇ podmıńka (3.21) platı´ v diskre´tnı´m prˇ´ıpadeˇ, tj. Z Ω

ph ∇ · vh = 0, ∀vh ∈ Xh0 =⇒ ph = konst.

Pokud prˇedchozı´ tvrzenı´ spojı´me s (3.39) a s vektorem p z (3.25) dosta´va´me

T B p, v = 0, ∀v ∈ Rnu =⇒ p = konst.

(3.40)

(3.41)

cozˇ je ekvivalentnı´ pozˇadavku null(B T ) = {1}. V dalsˇ´ı kapitole se budeme zaby´vat smı´sˇeny´mi prvky, ktere´ jsou stabilnı´ (kompatibilnı´) ve smyslu splneˇnı´ podmıńky (3.40) na dane´ sı´ti. Pozna´mka 3.1 Pokud je rychlost prˇedepsańa na cele´ hranici , pak 1 ∈ null(B T ), tj. diskre´tnı´ Stokesova matice (3.26) je singula´rnı´, a tudı´zˇ vektor na prave´ straneˇ (3.26) musı´ splnˇovat (2.27) v prˇ´ıpadeˇ, aby existovalo rˇesˇenı´. Uvazˇujeme, zˇe ph = p∗h = konst. a poznamenejme Z Z Z ∗ ∗ ph ∇ · vh = ph ∇ · vh = ph vh · n = 0, ∀vh ∈ Xh0 . Ω

Ω

∂Ω

To znamena´, zˇe p∗h ∈ Nh , z cˇehozˇ vyply´va´ 1 ∈ null(B T ). Pak nutna´ podmıńka pro rˇesˇenı´ syste´mu (3.26) je h1, gi = 0, tj. vektor g na prave´ straneˇ (3.26) musı´ v soucˇtu da´vat nulu. K hlubsˇ´ımu pochopenı´ se podı´va´me na definici (3.30) a vyuzˇijeme vlastnosti Neumannovy podmıńky (vycha´zı´ z tzv. Lemmatu o rozdeˇlenı´ jednicˇky), zˇe suma ba´zovyćh funkcı´ pro tlak je rovna jedne´: np X

gk =

np nu +n∂ X X

uj

k=1 j=1+nu

k=1

=

nu +n∂ X

uj

=

uj

j=1+nu

=

nu +n∂ X

uj

j=1+nu

=

Z

Z

Z

nu +n∂ X

∂Ω j=1+n

|

~j ψk ∇ · φ Ω

Z !X np Ω

j=1+nu nu +n∂ X

Z

ψk

k=1

"

~j ∇·φ

~j ∇·φ Ω

~j · n φ ∂Ω

~j · n uj φ

u

{z

=:w·n

}

Cozˇ ukazuje, zˇe podmıńka h1, gi = 0 je diskre´tnı´ reprezentace za´kona zachovańı´ hmoty (3.5). Jednodusˇe rˇecˇeno pro hranicˇnı´ data reprezentujıćı´ rychlost musı´ platit, zˇe hodnoty uj na vtokove´ a odtokove´ hranici musı´ by´t vybrańa. Tak aby objem kapaliny, jezˇ vstupuje

42

Obra´zek 14: Q2 − Q1 prvek (• dveˇ slozˇky rychlosti, ◦ tlak) do dane´ oblasti se rovnal objemu kapaliny vystupujıćı´mu z dane´ oblasti, tj. pro prˇ´ıtok a vy´tok musı´ platit, zˇe Z

nu +n∂ X

∂Ω− j=1+n

~j · n + uj φ

u

Z

nu +n∂ X

∂Ω+ j=1+n

~ j · n = 0, uj φ

(3.42)

u

jinak diskre´tnı´ Stokesova u´loha nenı´ dobrˇe podmıńeˇna´. Cozˇ take´ da´va´ uka´zku toho, procˇ Dirichletova podmıńka na vy´tokove´ hranici mu˚zˇe zpu˚sobit nemale´ potı´zˇe.

3.6

Stabilnı´ prvky Q2 − Q1

Hned na zacˇa´tku je nutno podotknout, zˇe existujı´ i nestabilnı´ prvky neboli takove´, ktere´ nesplnˇujı´ podmıńku diskre´tnı´ inf-sub podmıńku25 . Nestabilnı´ elementy jsou rˇesˇeny v [13], [14] a [27]. Uvazˇujme naprˇ. prˇ´ıklad uzavrˇene´ho proudeˇnı´ popsane´ho Stokesovy´mi rovnicemi. Hlavnı´ mysˇlenka je rozlozˇenı´ oblasti Ω na mensˇ´ı podoblasti M, jezˇ tvorˇ´ı male´ elementy usporˇa´dane´ v neˇjake´ obecne´ topologii. Metoda postavena´ na Q2 − Q1 elementech se cˇasto v literaturˇe oznacˇuje jako Taylor-Hoodova metoda. Pozice uzlu˚ jsou zna´zorneˇny na obra´zku (14). Jak bylo prˇedeslańo budeme se zaby´vat pouze Q2 − Q1 prvky, (14). Soustavu vıće prvku˚ nazveme jako makroprvek a budeme se zajı´mat, jestli je dany´ makroprvek stabilnı´. Oznacˇme M = k kvadraticky´ prvek z obra´zku 14. Uvazˇujme uzavrˇene´ proudeˇnı´ prˇes dany´ prvek M, pak vnitrˇnı´ uzel cˇ.9 bude uzlem pro rychlost a uzly 1,2,3,4 budou tlakove´ uzly. Pro vysˇetrˇenı´ stability na takove´m elementu si sestrojı´me matici B = [Bx , By ] pomocı´ (3.34) a (3.35) a na´sledneˇ oveˇrˇ´ıme, zda-li splnˇuje podmıńku 1 = null(B T ) ve smyslu (3.40) a (3.41). Uzˇ ted’ mu˚zˇeme tvrdit, zˇe syste´m nenı´ urcˇiteˇ stabilnı´, ale velmi hrubeˇ si uka´zˇeme uka´zku, procˇ nenı´ a na konci si rˇekneˇme, procˇ jsme schopni rozhodnout uzˇ v te´to chvı´li. Rychlost v v (3.41) ma´ dveˇ slozˇky rychlosti definovane´ v uzlu 9, zatı´mco tlak p ma´ rozmeˇr 1, ale je obsazˇen 4x na dane´m prvku. Tudı´zˇ Bx , By sloupcove´ vektory, Z ∂φ9 BxT = [bx,j ], bx,j = − ψj , (3.43) ∂x k 25

Cˇasty´ na´zev je Brezzi Babusˇka podmıńka.

43 Z

∂φ9 , (3.44) ∂x k kde φ9 je bikvadraticka´ funkce z podkapitoly 2.4. Tato funkce naby´va´ jednicˇky v uzlu 9 a je nulova´ v ostatnıćh 8 uzlech. Jednodusˇeji ψj , j = 1, 2, 3, 4 jsou standardnı´ bilinea´rnı´ funkce definovana´ na k . Mu˚zˇeme vycˇ´ıslit (3.43) a (3.44) integracı´, uzˇitı´m Simpsonova pravidla: Z |k | (f1 + f2 + f3 + f4 + 4f5 + 4f6 + 4f7 + 4f8 + 16f9 ), f (x, y) = 36 k ByT

= [by,j ],

by,j = −

ψj

kde fi odpovı´da´ f (xi ). To na´m da´va´ Z ∂ψj ∂ψj 4 φ9 = |k | (x9 ), bx,j = ∂x 9 ∂x k Z ∂ψj ∂ψj 4 = |k | (x9 ), φ9 by,j = ∂y 9 ∂y k Vycˇ´ıslenı´m linea´rnıćh funkcı´

∂ψj ∂ψj ∂x , ∂y , j

= 1, 2, 3, 4 pro centra´lnı´ uzel, naprˇ. pro j = 1:

hy ∂ψj (x9 ) = − , ∂x 2|k | cozˇ na´m da´va´ pro vycˇ´ıslenı´ pro vsˇechna j: T 2 −hy hy hy −hy Bx T B = = . ByT 9 −hx −hx hx hx

(3.45)

Pozna´mka 3.2 Diskre´tnı´ opera´tor BxT odpovı´da´ centra´lnı´ diferenci aproximace tlakove´ derivace ve smeˇru x v uzlu (x9 ): 2 hy (−p1 + p2 + p3 − p4 ) 0 = 9 (p2 + p3 )/2 − (p1 + p4 )/2 4 4 ∂p = |k | ≈ |k | (x9 ) (3.46) 9 hx 9 ∂x Z (3.46) vzply´va´, zˇe linea´rnı´ syste´m B T p = 0 je ekvivalentnı´ s −p1 + p2 + p3 − p4 = 0, −p1 − p2 + p3 + p4 = 0, cozˇ je splneˇno vzˇdy, kdyzˇ p1 = p3 a p2 = p4 . Z toho vyply´va´, zˇe B T ma´ dimenzi ja´dra rovnou 2 a na za´kladeˇ tohoto tvrzenı´ mu˚zˇeme tvrdit,zˇe makroelement M prˇedstavujıćı´ jeden Q2 − Q1 prvek nenı´ stabilnı´. Jak jsme prˇedeslali v pocˇa´tku, zˇe mu˚zˇeme rozhodnout zda-li je makroprvek M stabilnı´ a to jednoduchou u´vahou, zˇe syste´m B T p = 0 obsahuje dveˇ rovnice o cˇtyrˇech nezna´myćh. Odtud dimenze ja´dra B T musı´ by´t prˇi nejmensˇ´ım dveˇ. Prˇedchozı´ analy´za stability uka´zala, jak poznat jestli je makroprvek stabilnı´. Zjistili jsme, zˇe jeden Q2 − Q1 prvek nenı´26 , pokusı´me se tedy najı´t vhodneˇjsˇ´ı makroprvek. Na´sledujıćı´m prˇirozeneˇ slozˇiteˇjsˇ´ım makroprvkem je M = k ∪ m se spolecˇnou hranicı´, viz obra´zek (15) Opeˇt uvazˇujeme-li uzavrˇene´ proudeˇnı´ prˇes M, ma´me trˇi vnitrˇnı´ uzly pro 26

Cozˇ, ale vu˚bec neznamena´, zˇe Q2 − Q1 prvek je sˇpatna´ volba na diskretizaci Stokesovy u´lohy.

44

Obra´zek 15: Q2 − Q1 makroprvek M = k ∪ m rychlost 7, 8, 9 a sˇest tlakovyćh uzlu˚ 1 − 6. V tomto prˇ´ıpadeˇ je matice B T ∈ R6×6 , cozˇ by meˇlo poskytnout dostatecˇnou stabilitu dane´ho makroprvku, jezˇ si pokusı´me uka´zat. Zacˇneme od matice (3.45), jizˇ poskytl cˇtyrˇi rˇa´dky o stabiliteˇ syste´mu B T p = 0. Stejny´m postupem jako u prˇedchozı´ho makroprvku na´m vznikne matice B T pro trˇi vnitrˇnı´ rychlostnı´ uzly s tı´m, zˇe pro uzel 7 resp. 9 dosta´va´me podmıńky p1 = p5 , p2 = p6 resp. p2 = p4 , p3 = p5 .Kombinace teˇchto cˇtyrˇech podmıńek na´m da´va´ p1 = p3 = p5 = ξ, p2 = p4 = p6 = η. K potvrzenı´ stability dane´ho makroprvku musı´me uka´zat, zˇe ξ = η, za pouzˇitı´ dvou rovnic, jezˇ platı´ pro rychlostnı´ uzel 8. Vezmeme si nejdrˇ´ıve vertika´lnı´ slozˇku rychlosti, ∂p cozˇ se vyjadrˇuje diskre´tnı´ aproximaci ∂y (x8 ) = 0. Je zrˇejme´, zˇe rˇa´dky odpovı´dajıćı´ uzlu˚m 1, 3, 4 a 6 jsou nulove´. Uvazˇujme naprˇ. koeficient p1 , pak Z ∂φ8 T , By,k = − ψ1 1 ∂y M Z ∂φ8 = − ψ1 , (ψ1 = 0 v m ) ∂y k Z ∂ψ1 , (protozˇe φ8 ny = 0na hranici) = − φ8 ∂y k ∂ψ1 1 |k | (x8 ), (za pouzˇitı´ Simpsonova pravidla) = 9 ∂y = 0, nebot’ ψ1 = 0 na svisle´ hrane´ mezi uzly 2 a 5. Stejny´m postupem pro dalsˇ´ı koeficienty dosta´va´me 2 0 hx 0 0 hx 0 (ByT )k = (3.47) 9 Vra´tı´me-li se zpeˇt k syste´mu B T p = 0 z rovnice (3.47) vyply´va´ p2 = p5 , cozˇ take´ znamena´, zˇe ξ = η, z cˇehozˇ da´le plyne, zˇe 1 = null(B T ), cozˇ je chteˇny´ vy´sledek.

45

Q2−Q1 sit 1

0.8

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.8

−1 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Obra´zek 16: Q2 − Q1 sı´t’(◦ rychlost, × tlak) Prˇedchozı´ analy´za uka´zala, zˇe makroprvek M = k ∪ m slozˇeny´ ze dvou Q2 − Q1 prvku˚ je stabilnı´. Da´le bylo uka´zańo, zˇe stabilita je neza´visla´ na velikosti dvou sousednıćh prvku˚ ani na poloze, otocˇ´ıme-li prvek o 90o . Tato geometricka´ invariance je v podstateˇ klıćˇem k ustanovenı´ stability smı´sˇenyćh prvku˚ na obecnyćh sı´tıćh. Jakmile jednou potvrdı´me stabilitu na dane´m referencˇnı´m makroprvku, pak platı´ pro vsˇechny prvky stejne´ topologie. Navıć mu˚zˇeme tvrdit, obsahuje-li sı´t’sudy´ pocˇet prvku˚, pak je stabilnı´, protozˇe je rozlozˇitelna´ na makroprvek, jezˇ byl uveden v prˇedchozı´m prˇ´ıpadeˇ. V prˇ´ıpadeˇ liche´ho pocˇtu prvku˚, tedy prˇ´ıpad, kdy makroprvek bude M = k ∪ m ∪ l , se da´ doka´zat stabilita obdobneˇ. Jediny´ prˇ´ıpad liche´ho pocˇtu prvku˚, ktery´ nenı´ stabilnı´ je ojedineˇly´ Q2 − Q1 prvek.

3.7

ˇ esˇenı´ syste´mu rovnic R

V te´to podkapitole si uka´zˇeme rˇesˇenı´ syste´mu linea´rnıćh rovnic, jezˇ vznika´ prˇi rˇesˇenı´ Stokesovy rovnice. Tedy sedlo-bodovy´ syste´m matic, ktera´ vznika´ ze Stokesovy u´lohy naby´va´ tvaru A BT . (3.48) B -C Kde matice A = AT ∈ Rn×n je pozitivneˇ definitnı´ (A > 0), B ∈ Rm×n ma´ hodnost r ≤ m ≤ n a C = C T ∈ Rm×m27 je pozitivneˇ semidefinitnı´ (C ≥ 0). Indefinitnost syste´mu 27

V te´to praći se uzˇ´ıva´ stabilnıćh elementu˚ pro ktere´ je matice C = 0.

46

tlakove pole

Proudnice 1 100 0.5

0 −100

0

−200 −0.5

−300 1 1

−1 −1

0 −0.5

0

0.5

1

0 −1 −1

ˇ esˇenı´ Stokesovy u´lohy na jemne´ sı´ti Obra´zek 17: R na´m velmi zpomaluje iteracˇnı´ metody vyuzˇ´ıvajıćı´ Krylovu˚v prostor. Inspiraci pro rˇesˇenı´ (3.48) mu˚zˇeme hledat v cˇlańku [17], jedna´ se o zajı´mavy´ zpu˚sob prˇedpodmıńeˇnı´ (3.48) a na´sledne´ aplikace upravene´ metody sdruzˇenyćh gradientu˚. Nejedna´ se o „standardnı´ “ rˇesˇenı´, jaky´m je naprˇ´ıklad MINRES, ale o zajı´mavy´ pocˇin. Zejmeńa v knize [1] lze najı´t jak MINRES, tak vhodne´ prˇedpomıńeˇni pro syste´m (3.48). Protozˇe u´lohy pro testovanı´ z podkapitoly 3.2, nejsou nikterak slozˇite´, uzˇ´ıva´me pro hledańı´ rˇesˇenı´ „backslash“ opera´tor ´ lohy byly vybrańy za´meˇrneˇ tak, aby byly otesimplementovany´ v syste´mu MATLAB. U tovańy proble´my, jezˇ mohou vznikat u proudeˇnı´ (naprˇ. singularity v rozıćh, prˇirozena´ vy´tokova´ podmıńka atd.) a lze oveˇrˇ´ıt spra´vnost vy´sledku˚ dı´ky porovnańı´ s vy´sledky v [13] a [22]. 3.7.1

Numericke´ vy´pocˇty

V te´to cˇa´sti uka´zˇeme rˇesˇenı´ testovacıćh prˇ´ıkladu˚, jezˇ jsme rˇesˇili numericky. Budeme tedy rˇesˇit linea´rnı´ syste´m (3.26), poprˇ. (3.32), ktery´ vznika´ prˇi rˇesˇenı´ Stokesovy u´lohy. Prvnı´ jsme rˇesˇili Prˇ´ıklad 1: vtok-vy´tok z podkapitoly 3.2. A na obra´zku 11 a 12 mu˚zˇeme videˇt rˇesˇenı´ Stokesovy u´lohy. Uzˇ od pohledu lamina´rnı´ proudeˇnı´ je rˇesˇenı´m u´lohy, cozˇ se na´m take´ potvrdilo prˇi srovnańı´ s vy´sˇe uvedenou literaturou. Poznamenejme, zˇe u´lohu rˇesˇ´ıme na sı´ti z obra´zku 16. ´ lohu opeˇt rˇesˇ´ıme na stejne´ sı´ti Druhou u´lohou k rˇesˇenı´ bylo proudeˇnı´ v kaviteˇ. U z obra´zku 16 (4 × 4) a vy´sledek vidı´me na obra´zku 13. Na obra´zku 17 lze videˇt rˇesˇenı´ na jemneˇjsˇ´ı sı´ti (256 × 256), kde lze pozorovat singularity v rozıćh (mı´sta, kde tlak roste do nekonecˇna).

47

4

Navierova-Stokesova u´loha

V te´to kapitole se budeme zaby´vat Navierovy´mi-Stokesovy´mi rovnicemi28 a jejich matematicky´m modelem. Prˇi tvorbeˇ te´to kapitoly vycha´zı´me z [2], [13], [24], [21], [26] , [22], [25], [4], [11], [27], [18], [15] a [16]. Jak bylo prˇedeslańo v prvnı´ kapitole, Navier-Stokesovy rovnice popisujı´ rovnost setrvacˇne´ sı´la vu˚cˇi soucˇtu hmotnostnı´, tlakove´ a trˇecı´ sı´ly, tzv. Newtonskyćh kapalin (naprˇ. voda a vzduch). Vıće v [12] a [20]. Take´ mohou vycha´zet ze Stokesovy rovnice prˇidańı´m nelinea´rnı´ho konvektivnı´ho cˇlenu u · ∇u := (u · ∇) u (je du˚lezˇite´ si uveˇdomit, zˇe dı´ky nelinea´rnı´mu cˇlenu mu˚zˇe mı´t NS jako okrajova´ u´loha vıće nezˇ jedno stabilnı´ rˇesˇenı´!) a sı´ly f . Budeme se zaby´vat staciona´rnı´ Navierovou-Stokesovu rovnicı´ na oblasti Ω ⊂ Rn −ν∇2 u + (u · ∇) u + ∇p = f, ∇u = 0.

(4.1)

Kde u je rychlost (vektorova´ velicˇina), p je tlak (skala´r) a ν > 0 je kinematicka´ viskozita (konstanta). Prˇi pouzˇitı´ smı´sˇene´ metody konecˇnyćh prvku˚ rovnice (4.1) vedou ne nelinea´rnı´ syste´m algebraickyćh rovnic, ktere´ mohou by´t vyrˇesˇeny jen za pomocı´ iteracˇnıćh metod. Na´s bude zajı´mat hlavneˇ • Picardova iteracˇnı´ metoda, • Newtonova iteracˇnı´ metoda. Podobneˇ jako u Stokesovy rovnice uvazˇujme okrajove´ podmıńky na dvou a trˇ´ı rozmeˇrne´ oblasti Ω spolecˇneˇ s okrajovy´mi podmıńkami na ∂Ω = ∂ΩD ∪ ∂ΩN dane´ u = w na ∂ΩD ,

ν

∂u − np = s na ∂ΩN , ∂n

(4.2)

kde n je norma´lovy´ vektor smeˇrˇujıćı´ ven z oblasti a ∂u/∂n vyjadrˇuje derivaci ve smeˇru norma´ly. Pokud ∂ΩD ≡ ∂Ω, tj. rychlost, je prˇedepsańa na cele´ hranici, pak rˇesˇenı´ tlaku pro (4.1) a (4.2) je jednoznacˇneˇ urcˇeno azˇ na konstantu.

4.1

Rovnice konvekce-difuze

Budeme se jı´ veˇnovat, protozˇe s jejı´ pomocı´ ne jen vıće pochopı´me NS rovnice, ale take´ prˇedstavujı´ jeden z prostrˇedku˚, jak numericky pocˇ´ıtat NS rovnici. Rovnice ve sve´ obecne´ podobeˇ prˇedstavujı´ koncentraci u znecˇisteˇnı´ neˇjakou la´tkou ve smeˇru proudu o dane´ rychlosti w. Teplotu kapaliny pohybujıćı´ se pode´l teple´ steˇny nebo koncentraci elektronu˚ v polovodicˇovyćh soucˇa´stkaćh. Neˇkdy se take´ vyuzˇ´ıva´ advekce-difuze viz [25]. Vice o konvekci-difuzi v [16]. Uvazˇujme tedy rovnice konvekce-difuze w · ∇u − ǫ∇2 u = f, 28

V pozdeˇjsˇ´ım textu se take´ mu˚zˇe objevit zkratka NS.

(4.3)

48

kde ǫ > 0 je difuznı´ konstanta, pro NS rovnici prˇedstavuje viskozitu a u je vektorova´ funkce reprezentujıćı´ rychlost proudeˇnı´. Obecneˇ platı´, zˇe difuze je meńeˇ vy´znamny´ efekt v proudeˇnı´ kapalin nezˇ konvekce. Prˇ´ıkladem na´m mu˚zˇe by´t komıń vypousˇteˇjıćı´ prachove´ cˇa´stice o dane´ koncentraci. Pokud je bezveˇtrˇ´ı, jde cˇisteˇ jen o difuznı´ proble´m, ale pokud zacˇne foukat, bude dominantnı´m jevem konvekce. Zajı´ma´ na´s proudeˇnı´ kapalin ne jejich vlastnosti v klidu, bude tudı´zˇ i pro na´s konvekce du˚lezˇiteˇjsˇ´ı, tedy pokud ǫ ≪ |w|. Uvazˇujeme stejne´ hranicˇnı´ podmıńky jako v NS rovnici (4.2) a stejne´ deˇlenı´ hranici na vy´tokovou, charakteristickou a vtokovou hranici jako u Stokesovy u´lohy. Zameˇrˇ´ıme se nynı´ na konvekcˇnı´ cˇlen w · ∇u. Prˇedpokla´da´me, zˇe pomeˇr ǫ/(|w|L) je maly´, kde L oznacˇuje charakteristicke´ meˇrˇ´ıtko de´lky, jezˇ bude vysveˇtleno pozdeˇji. Rˇesˇenı´ (4.3) inklinuje na veˇtsˇineˇ oblastı´ k rˇesˇenı´ u b hyperbolicke´ rovnice w · ∇b u = f,

(4.4)

Rovnici nazveme redukovany´ proble´m. Diferencia´lnı´ opera´tor (4.4) je vysˇsˇ´ıho stupneˇ nezˇ (4.3) a u b obecneˇ nesplnˇuje okrajove´ podmıńky pro u. Podrobneˇjsˇ´ı zaby´vańı´ se rovnicı´ konvekce-difuze prˇesahuje svy´m rozsahem tuto praći, vıće lze nale´zt v [13], [16] a [22].

4.2

Bezrozmeˇrny´ tvar

Podobneˇ jako v prˇedchozı´ podkapitole je dobre´ mı´t neˇjake´ meˇrˇ´ıtko relativnıćh prˇ´ıspeˇvku˚ konvekce cˇ´ı difuze. Toho mu˚zˇeme dosa´hnout tak, zˇe „normalizujeme“ syste´m (4.1) s ohledem na velikost oblasti a rychlosti jako dane´ velicˇiny. K tomu si zavedeme dane´ L jako charakteristickou de´lku, naprˇ. jako konstantu z Poincare´ho nerovnosti. Vıće v [13]. Jestlizˇe body z Ω oznacˇ´ıme vektorem x , pak k = x/L znacˇ´ı body normalizovane´ oblasti. Podobneˇ pro rychlost bude platit u = u∗ U , kde U je referencˇnı´ rychlost, naprˇ. maxima´lnı´ rychlost na vtoku. Jestlizˇe tlak take´ opatrˇ´ıme meˇrˇ´ıtkem p(Lk) = U 2 p∗ (k) na normalizovane´ oblasti, pak rovnici (4.1), mu˚zˇeme prˇepsat do tvaru L 1 − ∇2 u∗ + (u∗ · ∇) u∗ + ∇p∗ = 2 f, R U kde relativnı´ prˇ´ıspeˇvky konvekce a difuze definujeme Reynoldsovy´m cˇ´ıslem

(4.5)

UL . (4.6) ν Na volbeˇ U a L za´visı´, zda je proble´m difuzneˇ dominantnı´ (R < 1) (pozn. pak lze uka´zat, zˇe rˇesˇenı´ je jednoznacˇneˇ definovańo) nebo konvektivneˇ dominantnı´ (R ≫ 1). Pokud ma´me velmi visko´znı´ kapalinu, vytva´rˇ´ı se v kapalineˇ tzv. smykove´ vrstvy, ktere´ jsou ve sve´ podstateˇ rˇesˇenı´m (4.1) a (4.2). R :=

Pozna´mka 4.1 Lze uka´zat [14], zˇe bezrozmeˇrna´ verze rovnice (4.1), prˇi pouzˇitı´ sˇka´lovańı´ tlaku p(Lk) = (νU/L)p∗ (k), da´va´ na´sledujıćı´ alternativu (4.5) −∇2 u∗ + Ru∗ · ∇u∗ + ∇pu∗ = f∗ A da´ se uka´zat, zˇe vezmeme-li limitu ν → ∞, dosta´va´me Stokesovu rovnici (3.2)

49

4.3

Slaba´ formulace a linearizace

Slaba´ formulace NS rovnice (4.1) a (4.2) je velmi podobna´ rovnicı´m, ktere´ jsme odvodili u slabe´ formulace Stokesovy u´lohy. Definujeme stejne´ prostory pro rˇesˇenı´ a testovacı´ funkce H1E H1E0

:= {u ∈ H 1 (Ω)d |u = w na ∂ΩD }, 1

d

:= {v ∈ H (Ω) |v = 0 na ∂ΩD },

(4.7) (4.8)

Slaba´ formulace, odvozeni v [26] a du˚kaz existence slabe´ho rˇesˇenı´ pro dvourozmeˇrnou NS rovnici v [21]. hleda´me u ∈ H1E a p ∈ L2 (Ω) takove´, zˇe Z Z Z Z ν ∇u : ∇v + (u · ∇u)v − p (∇ · v) = f·v ∀v ∈ H1E0 , (4.9) Ω Ω Ω Ω Z q (∇ · u) = 0 ∀q ∈ L2 (Ω). (4.10) Ω

Z identity (3.13) je zrˇejme´, zˇe kazˇde´ rˇesˇenı´ (4.1) a (4.2) splnˇuje (4.9) a (4.10). Jak bylo prˇedeslańo, hlavnı´ rozdı´l mezi Stokesovou a NS rovnicı´ je ve konvektivnı´m cˇlenu. Trochu prˇedbeˇhneme a prozradı´me, zˇe numericky budeme NS rovnici rˇesˇit tak, zˇe nejdrˇ´ıve vypocˇ´ıta´me jejı´ linearizovany´ prˇ´ıpad, neboli Stokesovu rovnici a na´sledneˇ pak dopocˇ´ıta´me „diference“, tj. rozdı´l v tlacıćh a rychlostech na jednotlivyćh elementech. Tento rozdı´l je ve trˇ´ı-linea´rnı´m funkciona´lu c : H1E0 × H1E0 × H1E0 → R definovane´m jako Z c(z, u, v) =

(z · ∇u) · v

(4.11)

Ω

Pro spojenı´ mezi podkapitolou konvekce difuze si zavedeme podprostor s nulovou divergencı´ jako (4.12) VE0 := {z ∈ H1E0 |∇ · z = 0 v Ω}, a definujeme si bilinea´rnı´ formu az (·, ·) : VE0 × VE0 → R, jako Z az (u, v) := ν ∇u : ∇v + c(z, u, v),

(4.13)

Ω

kde z je dana´ funkce z prostoru VE0 . Jako v podkapitole konvekce difuze potrˇebujeme spojitost a omezenost bilinea´rnı´ formy v normeˇ k∇vk. Za´kladem je, zˇe konvekce je antisymetricka´: c(z, u, v) = −c(z, v, u) prˇes prostor VE0 , cozˇ znamena´ c(z, u, u) = 0 a z toho plynoucı´ elipticita (pozitivnı´ definitnost, viz definice A.3) az (u, u) = νk∇vk2

∀u ∈ VE0 .

(4.14)

c(z, u, u) ≤ Γk∇zkk∇ukk∇vk

(4.15)

a na´sledna´ omezenost formy (viz definice A.1)

50

du˚kaz v [14]. Nynı´ definujeme u := u1 − u2 a p := p1 − p2 , kde (u1 , p1 ) je odlisˇne´ rˇesˇenı´ nezˇ (u2 , p2 ). Ale pozor platı´, totizˇ c(u1 , u1 , v) − c(u2 , u2 , v) 6= c(u1 − u2 , u1 − u2 , v) To znamena´, zˇe nejsme schopni pouzˇ´ıt homogennı´ho proble´mu k urcˇenı´ jednoznacˇnosti rˇesˇenı´ jako v prˇ´ıpadeˇ (3.19) a (3.20). Mı´sto toho, vzhledem k existenci alesponˇ jednoho slabe´ho rˇesˇenı´ splnˇujıćı´ (4.9) a (4.10), je mozˇne´ zadat podmıńky pro silovou funkci f a hranicˇnı´ data w, ktere´ zajisˇt’ujı´ jednoznacˇnost rˇesˇenı´ (azˇ na konstantu tlaku v prˇ´ıpadeˇ ∂Ω = ∂ΩD ). Studiem jednoznacˇnosti se zaby´vajı´ praće [21] a [26]. Ota´zky (ne)existence a (ne)jednoznacˇnost rˇesˇenı´ staciona´rnıćh stavu NS rovnic u´zce souvisı´ s klasickou koncepcı´ hydrodynamicke´ stability a bifurkacˇnı´ teorie. V na´sledujıćı´ podkapitole uvedeme neˇktere´ za´kladnı´ mysˇlenky. Take´ zavedeme pojem linearizace, ktery´ souvisı´ s touto podkapitolou a bude prˇedmeˇtem studia v na´sledne´ podkapitole.

4.4

Bifurkacˇnı´ analy´za a teorie stability

Bifurkace obecneˇ body, v nichzˇ docha´zı´ k nejednoznacˇnosti rˇesˇenı´. Za´klady te´to teorie najdeme v [4], [11] a [18]. Prˇedpokla´dejme v te´to podkapitole, zˇe slabe´ rˇesˇenı´ u(x) a p(x) splnˇuje staciona´rnı´ formulaci (4.9) a (4.10). Budeme se odkazovat na toto rˇesˇenı´ jako na „pevny´ bod“. Mysˇlenka teorie stability, motivovańa prakticky´mi experimenty, je studium vy´voje odchylek od tohoto pevne´ho bodu. Pokud odchylky od stabilnı´ho rˇesˇenı´ jsou zesilovańy s cˇasem, pak pevny´ bod je nestabilnı´ s cˇasem a vy´voj je prˇirozeneˇ nesta´ly´. Naproti tomu pokud se odchylky ustalujı´ s rostoucı´m cˇasem, tı´m pa´dem pevny´ bod je je stabilnı´ a vypocˇ´ıtane´ staciona´rnı´ rˇesˇenı´ je kandida´tem na cˇasoveˇ neza´visle´ rˇesˇenı´ proudıćı´ kapaliny (pozn. tato teorie by byla stejna´ i pro prˇ´ıpad nestaciona´rnı´ho proudeˇnı´). Definujme cˇasoveˇ za´visle´ rˇesˇenı´ proudeˇnı´ tak, zˇe v(x, t) = u(x) + δu(x, t), q(x, t) = p(x) + δp(x, t).

(4.16)

Prˇedpokla´dejme, zˇe rychlost v(x, t) splnˇuje Dirichletovy hranicˇnı´ podmıńky (4.2) tak, zˇe odchylka δu(·, t) je nulova´ na cele´ hranici. Jestlizˇe toto rˇesˇenı´ dosadı´me do nestaciona´rnı´ u´lohy (cˇasoveˇ za´visle´), dosta´va´me analogii slabe´ formulace. Pak odchylky δu(x) a p(x) vyhovujı´ pocˇa´tecˇnı´mu hranicˇnı´mu proble´mu: hleda´me δu(x, t) a p(x, t) takove´, zˇe Z Z Z ∂δu · v + D(u, δu, v) + ν ∇δu : ∇v − δp (∇ · v) = 0 (4.17) Ω Ω Ω ∂t Z q (∇ · δu) = 0. (4.18) Ω

51

pro vsˇechna v ∈ HE01 a q ∈ L2 (Ω), kde D(u, δu, v) je diference nelinea´rnıćh vy´razu˚, to jest Z (δu + u) · ∇(δu + u) · v − (u · ∇u) · v Ω ZΩ Z Z = (δu · ∇δu) · v + (δu · ∇u) · v + (u · ∇δu) · v

D(u, δu, v) =

Z

Ω

Ω

Ω

= c(δu, δu, v) + c(δu, u, v) + c(u, δu, v)

(4.19)

Cı´lem teorie stability je zı´skat informaci o vy´voji odchylek anizˇ by bylo trˇeba vyrˇesˇit (4.17)-(4.19). Klasicka´ „linea´rnı´ stabilita“ se da´ studovat, pokud odstranı´me kvadraticky´ cˇlen c(δu, δu, v) z (4.19). Za prˇedpokladu, zˇe odchylky rostou nebo rozpadajı´ se v cˇase, pak prˇedefinujeme δu(x, t) ← eλt δu(x) a δp(x, t) ← eλt δp(x). To vede na proble´m hledańı´ vlastnıćh cˇ´ısel a odpovı´dajıćıćh vlastnıćh vekoru˚ δu, δp, takovyćh, zˇe Z Z Z c(δu, u, v) + c(u, δu, v) + ν ∇δu : ∇v − δp (∇ · v) = λ δu · v Ω Ω ZΩ q (∇ · δu) = 0. (4.20) Ω

Da´ se uka´zat, pokud vlastnı´ cˇ´ısla jsou za´porna´, pak rˇesˇenı´ u, p pu˚vodnı´ho proudeˇnı´ je (linea´rneˇ) stabilnı´. Naprˇ´ıklad Poiseulovo proudeˇnı´ z nasˇeho prˇ´ıkladu staciona´rnı´ho proudeˇnı´ je nestabilnı´ pro R = 1/ν > 5772. Vıće v [7] a [24]. Linea´rnı´ stabilita je jedina´ nezbytna´ podmıńka pro celou nelinea´rnı´ stabilitu, tj. pokud u´loha nenı´ linea´rneˇ stabilnı´, pak ani nelinea´rnı´ u´loha nenı´ stabilnı´. Vhodna´ podmıńka linea´rnı´ stability se nazy´va´ energeticka´ analy´za, kterou si naznacˇ´ıme da´le. Vezmeme nelinea´rnı´ formulaci (4.17)-(4.18) s Dirichletovou hranicˇnı´ podmıńkou, pak testovacı´ funkce v jsou nulove´ na hranici. A fixujeme cˇas t. Za prˇedpokladu, zˇe δu ∈ H1E0 , mu˚zˇeme vzı´t v = δu. Jestlizˇe zvolı´me, za q = δp, pak poslednı´ cˇlen na leve´ straneˇ rovnosti (4.17) vypadne, tı´m pa´dem dosta´va´me Z ∂δu · δu + D(u, δu, δu) + νk∇δuk2 = 0. (4.21) ∂t Ω Vsˇimneˇme si, zˇe

Z

d ∂δu · δu = dt Ω ∂t

1 kδuk2 2

(4.22)

Pro zjednodusˇenı´ druhe´ho cˇlenu v (4.21), zvolı´me q = ∇·δu cˇ´ımzˇ ∇·δu = 0 v Ω. Pak uzˇitı´m antisymetrie lze uka´zat, zˇe c(δu, δu, u) = 0. Jednodusˇeji, pokud ∇ · u = 0 v Ω dosta´va´me c(u, δu, δu) = 0 a uzˇitı´m (4.19), lze odvodit, zˇe D(u, δu, δu) = c(δu, u, δu) = −c(δu, δu, u). Dosazenı´m do (4.21) a pouzˇitı´m (4.22) na´s vede k na´sledujıćı´mu odhadu pro rozpad odchylky rychlosti v L2 normeˇ. d 1 2 kδuk = c(δu, δu, u) − νk∇δuk2 . (4.23) dt 2

52

V tuto chvı´li upustı´me od bezrozmeˇrne´ho nelinea´rnı´ formulace (4.17)-(4.19). Specia´lneˇ pro bezrozmeˇrne´ velicˇiny u∗ = u/U,

p∗ = pL/νU,

t∗ = tν/L2 ,

k = x/L,

kde L je Poincare´ho konstanta (viz Poincare–Friedrichsova nerovnost v apendixu). Prˇepsańı´m (4.23) a uzˇitı´m pozna´mky (4.1), dosta´va´me tvar za´visly´ na jedine´ promeˇnne´ d 1 kδuk2 = Rc(δu, δu, u) − k∇δuk2 . (4.24) dt 2 Dosta´va´me se do momentu, kdy R je jediny´m parametrem. Poznamenejme, zˇe je-li R nulove´, pak prava´ strana (4.24) je za´porna´ a vesˇkere´ odchylky stejnomeˇrneˇ uby´vajı´ s rostoucı´m cˇasem. Zde je patrny´ vliv Reynoldsova cˇ´ısla na stabilitu rˇesˇenı´. Poznamenejme, zˇe bezrozmeˇrnost rychlosti ma´ takove´ meˇrˇ´ıtko, aby kuk∞ = 1. Pouzˇijeme CauchySchvarzovu nerovnost (viz [3]) na cˇlen konvekce obdrzˇ´ıme Z 1 1 2 2 c(δu, δu, u) = δu · ∇δu · u ≤ kδukk∇δukkuk ≤ Rkδuk + k∇δuk . 2 R Ω Pouzˇijeme-li skutecˇnosti, zˇe prˇepocˇet pomocı´ Poincare´ho lemmatu [musis to poresit] (kδuk ≤ k∇δuk) s ohledem na L, zjistı´me pak d 1 1 1 R2 2 kδuk ≤ kδuk2 − k∇δuk2 ≤ (R2 − 1) kδuk2 . dt 2 2 2 2 Odtud, definujeme-li si K(t) := kδu(·, t)k2 jako energii odchylky, dosta´va´me linea´rnı´ diferencia´lnı´ nerovnost dK + (1 − R2 )K ≤ 0. (4.25) dt K vyrˇesˇenı´, potrˇebujeme kladny´ integracˇnı´ faktor e(1−R To na´m da´ hranici u´bytku energie

2 )t

a integrujeme od cˇasu t = 0.

2

K(t) ≤ K(0)e−t(1−R ) ,

(4.26)

odtud lze pozorovat, zˇe dostacˇujıćı´ podmıńkou pro u´bytek vesˇkeryćh odchylek je R < 1. V tomto prˇ´ıpadeˇ rˇekneme, zˇe je proudeˇnı´ nelinea´rneˇ stabilnı´. V proudeˇnı´ se vyuzˇ´ıva´ tzv. kriticka´ hodnota Reynoldsova cˇ´ısla. Tvorˇ´ı hranici, kdy je proudeˇnı´ (prˇi urcˇite´ rychlosti, vazkosti a tvaru oblasti) lamina´rnı´ a kdy prˇejde v turbulentnı´. Matematicky se tento jev prˇechodu modeluje bifurkacemi,[4] a [18]. Bifurkace jsou mı´sta s nejednoznacˇnostı´ rˇesˇenı´. Bifurkacı´ je mnoho typu˚ u NS rovnic je du˚lezˇita´ Hopfova bifurkace a bifurkace vidle29 . Prvnı´ Hopfova se vyuzˇ´ıva´ v mı´stech, kde vznika´ stabilnı´ periodicke´ rˇesˇenı´. A druha´ v mı´stech, kde docha´zı´ k veˇtvenı´ rˇesˇenı´, naprˇ. proudeˇnı´ na L oblasti na rohu nebo u proudeˇnı´ v trubici, ktera´ se ve smeˇru proudu rozsˇirˇuje. Proble´m na´m zpu˚sobujı´ zejmeńa z hlediska aproximace, take´ majı´ vliv na konvergenci nelinea´rnıćh iteracˇnıćh metod jezˇ budou uka´zańy pozdeˇji. 29

Anglicky pitchfork jedna´ se o prˇeklad autora DP, nenalezen cˇesky´ ekvivalent.

53

4.5

Newtonova a Picardova metoda

NS rovnice budeme rˇesˇit tak, zˇe si nejdrˇ´ıve vyrˇesˇ´ıme linearizovany´ proble´m v nasˇem prˇ´ıpadeˇ Stokesovu rovnici, kdy rˇesˇenı´ bude (u0 , p0 ) ∈ HE × L2 (Ω) a na´sledneˇ vypocˇ´ıta´me sekvencı´ nelinea´rnıćh iteracı´ (u0 , p0 ), (u1 , p1 ), . . . ∈ HE × L2 (Ω), ktere´ „snad“ budou konvergovat ke slabe´mu rˇesˇenı´. Prˇi tvorbeˇ te´to podkapitoly jsme vycha´zeli [13], [22], [14], [27] a [16]. Newtonova iterace, jako jedna z klasickyćh metod, na´m da´va´ rˇesˇenı´ (uk , pk ). Zacˇ´ına´ vy´pocˇtem nelinea´rnı´ho rezidua ze slabe´ formulace (4.9)-(4.10). Dosta´va´me Rk (v) a rk (q) splnˇujıćı´ Z Z Z Rk = f · v − c(uk , uk , v) − ν ∇uk : ∇v + p (∇ · v), (4.27) ΩZ Ω Ω rk = − q (∇ · u), (4.28) Ω

pro vsˇechna v ∈ H1E0 a q ∈ L2 (Ω). Kde u = uk + δuk a p = pk + δpk je rˇesˇenı´ (4.9)-(4.10). Je snadne´ nahle´dnout, zˇe odchylky δuk ∈ H1E0 a δpk ∈ L2 (Ω) splnˇujı´ Z Z Rk (v) = D(uk , δuk , v) + ν ∇δuk : ∇v − δpk (∇ · v), (4.29) Ω Ω Z rk (q) = q (∇ · ∇uk ), (4.30) Ω

pro vsˇechna v ∈ H1E0 a q ∈ L2 (Ω), kde D(uk , δuk , v) je diference nelinea´rnı´ho vy´razu, jako v (4.19). Rozepsańı´m cˇlenu D(uk , δuk , v) a odstraneˇnı´m kvadraticke´ho cˇlenu c(δuk , δuk , v), prˇesneˇ jako v odvozenı´ (4.20), dosta´va´me linea´rnı´ u´lohu pro vcˇechna v ∈ H1E0 a q ∈ L2 (Ω), hleda´me δuk ∈ H1E0 a δpk ∈ L2 (Ω) splnˇujıćı´ Z Z Rk (v) = c(uk , δuk , v) + c(δuk , uk , v) + ν ∇δuk : ∇v − δpk (∇ · v), Ω Ω Z rk (q) = q (∇ · ∇uk ). (4.31) Ω

ˇ esˇenı´ (4.31) se nazy´va´ Newnonova korekce. Aktualizacı´ prˇedchozı´ iterace pomocı´ uk+1 = R uk + δuk a pk+1 = pk + δpk definuje dalsˇ´ı iteracˇnı´ krok. V dalsˇ´ı podkapitole se zameˇrˇ´ıme na konecˇnneˇ dimenziona´lnı´ analogii Newtonovy iteracˇnı´ metody. Ale uzˇ ted’ je jasne´, zˇe Jacobiho matice je regula´rnı´ matice v pevne´m bodeˇ, pak tedy loka´lneˇ Newtonova iteracˇnı´ metoda konverguje kvadraticky. Nejdulezˇitejsˇ´ı podmıńkou je, aby vstupnı´ vektor byl dostatecˇneˇ blı´zko dane´mu rˇesˇenı´. Dalsˇ´ı metodou pro linearizaci je Picardova iteracˇnı´ metoda. Ve vy´razech (4.27) a (4.28) kvadraticky´ vy´raz c(δuk , δuk , v) klesa´ spolu s linea´rnı´m vy´razem c(δuk , uk , v) . Odtud mı´sto (4.31) ma´me na´sledujıćı´ linea´rnı´ proble´m pro vcˇechna v ∈ H1E0 a q ∈ L2 (Ω), hleda´me

54

δuk ∈ H1E0 a δpk ∈ L2 (Ω) splnˇujıćı´ Rk (v) = c(uk , δuk , v) + ν Z rk (q) = q (∇ · ∇uk ).

Z

∇δuk : ∇v − Ω

Z

δpk (∇ · v), Ω

(4.32)

Ω

ˇ esˇenı´ (4.32) se nazy´va´ Picardova korekce. Aktualizacı´ prˇedchozı´ iterace pomocı´ uk+1 = R uk + δuk a pk+1 = pk + δpk definuje dalsˇ´ı iteracˇnı´ krok pro vsˇechna v ∈ H1E0 a q ∈ L2 (Ω), hleda´me uk+1 ∈ H1E0 a pk+1 ∈ L2 (Ω) splnˇujıćı´ Z Z Z c(uk , uk+1 , v) + ν ∇uk+1 : ∇v − pk+1 (∇ · v) = f · v, (4.33) Ω Ω Z Ω q (∇ · ∇uk+1 ) = 0. (4.34) Ω

ˇ a´d konvergence Picardovi Formulace (4.33) a (4.34) se obecneˇ nazy´va´ Oseenu˚v syste´m. R iteracˇnı´ metody je „pouze“ linea´rnı´. Poznamenejme, zˇe vstupnı´ tlak p0 nenı´ potrˇeba, jestlizˇe je iteracˇnı´ krok je zapsań v alternativnı´ podobeˇ (4.33) a (4.34). Prˇi srovnańı´ na´m hlavneˇ vadı´, zˇe polomeˇr koule konvergence Newtonovy metody je typicky u´meˇrny´ viskoziteˇ. Pokud se zveˇtsˇuje Reynoldsovo cˇ´ıslo, potrˇebujeme sta´le lepsˇ´ı odhady vstupnıćh hodnot, aby Newtonova metoda konvergovala. Naproti tomu vy´hoda Picardova iteracˇnı´ metoda je velky´ polomeˇr konvergence v pomeˇru k Newtonoveˇ iteracˇnı´ metodeˇ.

4.6

Aproximace smı´sˇenou metodou konecˇnyćh prvku˚

V te´to podkapitole se budeme zaby´vat sestavenı´m Galerkinovy matice. Budeme se opı´rat o Stokesovu kapitolu, jen s tı´m rozdı´lem, zˇe zde ma´me „problematicky´“ cˇlen konvekce, ktery´ je trˇeba linearizovat, jezˇ byl du˚sledneˇji rozebrań v podkapitole konvekce a difuze. Nahradı´me prostory H1E0 a L2 (Ω) konecˇneˇ dimenziona´lnı´mi podprostory Xh0 ⊂ HE0 a M h ⊂ L2 (Ω). Tak, zˇe ze slabe´ formulace (4.9) a (4.10) se stane soustava linea´rnıćh rovnic. Hleda´me uh ∈ XhE a p ∈ M h takove´, zˇe Z Z Z Z ν ∇uh : ∇vh + (uh · ∇uh )vh − ph (∇ · vh ) = f · vh ∀vh ∈ Xh0 , (4.35) Ω Ω Ω Ω Z qh (∇ · uh ) = 0 ∀qh ∈ L2 (Ω). (4.36) Ω

Klasicka´ implementace ba´zovyćh prvku˚, ktera´ byla uka´zańa v podkapitole Galerkinova metoda a kapitole Stokesova u´loha vede na nelinea´rnı´ syste´m algebraickyćh rovnic. Linearizace teˇchto rovnic pomocı´ Newtnovy metody na´m da´va´ konecˇneˇ dimenziona´lnı´ analogii (4.31)

55 hleda´me odchylky δXh0 ∈ XhE a δp ∈ M h takove´, zˇe Rk (vh ) = c(uh , δuh , vh ) + c(δuh , uh , vh ) + ν Z rh (qh ) = qh (∇ · ∇uh ).

Z

∇δuh : ∇vh − Ω

Z

δph (∇ · vh ), Ω

(4.37)

Ω

pro vsˇechna vh ∈ Xh0 a qh ∈ M h . Kde Rk (vh ) a rh (qh ) tvorˇ´ı nelinea´rnı´ rezidua diskre´tnı´ formulace (4.35) a (4.36) (v prˇedchozı´ rovnici jsme pro prˇehlednost vynechali cˇlen k). Odstraneˇnı´m vy´razu c(δuh , uh , vh ) dosta´va´me diskre´tnı´ analogii Picardovi iterace. Uzˇitı´m vhodnyćh ba´zovyćh funkcı´ (u´lohu rˇesˇ´ıme na Q2 − Q1 prvcıćh) lze kazˇdou funkci rychlosti vyja´drˇit jako uh =

nu X

nu +n∂ X

~j + uj φ

j=1

~j , uj φ

δuh =

j=nu +1

nu X

~j , ∆uj φ

(4.38)

j=1

Fixujeme koeficienty uj : j = nu + 1, . . . , nu + n∂ tak, zˇe druhy´ vy´raz v (3.24) prˇedepisuje ~j hodnoty nehomogennı´ Dirichletovy podmıńky na ∂ΩD . Poznamenejme, zˇe mnozˇina φ reprezentuje ba´zove´ funkce prostoru, v neˇmzˇ hleda´me rˇesˇenı´ pro rychlost a oproti tomu mnozˇina ψj reprezentuje ba´zove´ funkce prostoru, v neˇmzˇ hleda´me rˇesˇenı´ pro tlak. Odtud ph =

np X

pk ψk ,

δph =

k=1

np X

∆pk ψk .

(4.39)

k=1

Dosazenı´m do (4.37) dosta´va´me soustavu linea´rnıćh rovnic ∆u f νA + N + W B T = . B 0 ∆p g

(4.40)

Matice A a B jsou definovańy jako (3.27) a (3.28). Matice N prˇedstavuje vektor konvekce a matice W prˇedstavuje derivace rychlosti ve smeˇru. Obeˇ za´visı´ na aktua´lnı´m odhadu diskre´tnı´ rychlosti uh . Prˇicˇemzˇ Z ~j ) · φ ~i, N = [nij ], nij = (uh · ∇φ (4.41) Ω

W = [wij ],

wij =

Z

~ j · ∇uh ) · φ ~i, (φ

(4.42)

Ω

pro vsˇechna i, j = 1, . . . , nu . Poznamenejme, zˇe matice W je symetricka´. Prava´ strana (4.40) prˇedstavuje vektory, cozˇ jsou nelinea´rnı´ rezidua spojena´ s aktua´lnı´m odhadem rychlosti uh a ph , definovana´ na´sledujıćı´m zpu˚sobem Z Z ~ ~i f = [fi ], fi = f · φi − uh · ∇uh · φ Ω Ω Z Z ~ i + ph (∇ · φ ~ i ), −ν ∇uh : ∇φ (4.43) Ω Ω Z g = [gk ], gk = ψk (∇ · uh ). (4.44) Ω

56

Syste´m (4.40) nazveme diskre´tnı´ Newtonu˚v proble´m. Pro Pickardovu iteracˇnı´ metodu vypustı´me matici W, cˇ´ımzˇ dosta´va´me νA + N B T ∆u f = . (4.45) B 0 ∆p g Syste´m (4.40) nazveme diskre´tnı´ Ossenu˚v proble´m. Protozˇe rˇesˇ´ıme syste´m (4.40) jako dvourozmeˇrnou u´lohu, lze prˇepsat (4.40) do tvaru      νA + N + Wxx Wxy BxT ∆ux fx  Wyx νA + N + Wyy ByT   ∆uy  =  fy  , (4.46) ∆p g Bx By 0

kde matice N je n × n skala´rnı´ matice konvekce Z ~j ) · φ ~i, N = [nij ], nij = (uh · ∇φ

(4.47)

Ω

ktera´ je studovańa v podkapitole konvekce- difuze. Matice Wxx , Wxy , Wyx a Wyy jsou rˇa´du n × n a prˇedstavu˚jı´ slabou derivaci aktua´lnı´ rychlosti ve smeˇru x a y naprˇ. definovańı´m x-ove´ slozˇky rychlosti uh jako ux dosta´va´me Z ∂ux Wxy,ij = [wxy,ij ], wxy,ij = φi φj . (4.48) Ω ∂y

4.7

´ lohy k rˇesˇenı´ U

´ lohu Budeme rˇesˇit totozˇnou u´lohu, jako u Stokesovy rovnice a to proudeˇnı´ v kaviteˇ. U rˇesˇ´ıme na sı´ti 256 × 256, pomocı´ dvou Pickardovyćh iteracˇnıćh kroku˚ a cˇtyrˇech Newtonovyćh iteracˇnıćh kroku˚. Pote´ je splneˇna podmıńka konvergence kterou jsme si zvolili, tak aby norma rezidujı´ (4.37) byla mensˇ´ı nezˇ na´mi pozˇadovana´ prˇesnost rˇesˇenı´ (1e − 8). Na obra´zku 18 lze videˇt rˇesˇenı´ u´lohy.

57

tlakove pole

Proudnice 1 −93 0.5 −93.5 0

−94

−0.5

−94.5 1 1

−1 −1

0 −0.5

0

0.5

1

0 −1 −1

Obra´zek 18: Rˇesˇenı´ Navierovy-Stokesovy u´lohy na 256 × 256 sı´ti

58

5

Za´veˇr

V te´to praći jsme si prˇedstavili za´klady numericke´ho rˇesˇenı´ Navierovyćh-Stokesovyćh rovnic pomocı´ smı´sˇene´ metody konecˇnyćh prvku˚ a uka´zali jsme si pouzˇitı´ te´to metody pro rˇesˇenı´ dvourozmeˇrne´ modelove´ u´lohy. Hlavnı´m cı´lem bylo porozumeˇnı´ modelovańı´ proudeˇnı´ tekutin a to od obecnyćh fyzika´lnıćh za´konitostı´ azˇ po obecne´ matematicke´ vyja´drˇenı´ vlastnostı´ tekutin, kdy v prvnı´ kapitole byly prˇedstaveny rovnice popisujıćı´ idea´lnı´ tekutiny v klidu azˇ po Navierovu-Stokesovu rovnici, popisujıćı´ turbulentnı´ (vy´rˇive´) proudeˇnı´. Pote´ jsme ve druhe´ kapitole formulovali za´kladnı´ vlastnosti variacˇnıćh metod pro jednorozmeˇrnou a dvourozmeˇrnou u´lohu. Z variacˇnı´ formulace jsme prˇesˇli na metodu konecˇnyćh prvku˚, specia´lneˇ na Galerkinovu metodu, odkud uzˇ byl jen kru˚cˇek ke smı´sˇene´ formulaci metody konecˇnyćh prvku˚, u nichzˇ byly popsańy jejich vlastnosti tak, aby bylo mozˇne´ formulovat vhodne´ modelove´ u´lohy pro numerickou implementaci. Trˇetı´ kapitola je veˇnovańa Stokesoveˇ u´loze. V prvnı´ pa´sti trˇetı´ kapitoly byly formulovańy vlastnosti Stokesovy u´lohy, za kteryćh je u´loha rˇesˇitelna´. Da´le byly uka´zańy vhodne´ prvky pro diskretizaci u´lohy, vcˇetneˇ oveˇrˇenı´ podmıńek stability danyćh prvku˚. Nynı´ uzˇ nic nebrańilo samotne´ numericke´ implementaci Stokesovy u´lohy, jejı´zˇ vy´sledek lze najı´t na konci kapitoly. Poslednı´m a steˇzˇejnı´m te´matem diplomove´ praće byla Navierova-Stokesova u´loha, lisˇ´ıcı´ se od Stokesovy u´lohy prˇedevsˇ´ım konvekcˇnı´m cˇlenem, ktery´ je popsań v cˇa´sti konvekce-difuze. Da´le jsme se v te´to kapitole snazˇili uka´zat alesponˇ za´kladnı´ za´konitosti teorie stability NS rovnic. Za´veˇr kapitoly je zameˇrˇen na Pickardovu a Newtonovu iteracˇnı´ metodu pro rˇesˇenı´ nelinea´rnıćh soustav, ktere´ jsou potrˇeba k numericke´mu rˇesˇenı´ NS rovnic. V za´veˇru praće bylo uka´zańo numericke´ rˇesˇenı´ modelove´ u´lohy. Numericka´ implementace Navierovyćh-Stokesovyćh rovnic na´m da´va´ silny´ apara´t pro modelovańı´ proudeˇnı´. Doka´zˇe na´m poskytnout prˇesne´ informace o tekutineˇ v danyćh mı´stech, cozˇ je obrovska´ prˇednost prˇed laboratornı´m meˇrˇenı´m. Dokonce i prˇed metodou konecˇnyćh objemu˚, ktera´ ma´ sve´ limity, pokud pozˇadujeme velkou prˇesnost, cenu, jizˇ musı´me zaplatit, je vy´pocˇetnı´ na´rocˇnost. Dalsˇ´ı nespornou vy´hodou vlastnı´ implementace oproti pouzˇitı´ softwaru pro proudeˇnı´ (naprˇ. ANSYS) je analy´za chyb a modelovanı´ specifickyćh u´loh, ktere´ nejsou implementovańy v dane´m softwaru. Diplomova´ praće tvorˇ´ı za´klad pro dalsˇ´ı rozvoj v mnoha oblastech. Numerickou cˇa´st lze rozsˇ´ırˇit o implementaci dalsˇ´ıho rozmeˇru (myslı´me trojrozmeˇrnou u´lohu). Da´le mu˚zˇeme uvazˇovat o nestaciona´rnı´ u´loze nebo o slozˇiteˇjsˇ´ıch u´lohaćh, ktere´ by vyzˇadovaly urcˇitou formu rozlozˇenı´ oblasti nebo paralelizace. V teorii stability lze implementovat i dalsˇ´ı prvky. Z teˇch stabilnıćh jmenujme naprˇ. Q2 −P0 , P2 −P1 atd., z nestabilnıćh naprˇ. Q1 −Q1 . Dalsˇ´ım smeˇrem mu˚zˇe by´t implementace vhodnyćh metod a prˇedpodmıńeˇnı´ pro rychle´ rˇesˇenı´ soustav, jezˇ vznikajı´ prˇi diskretizaci u´lohy. Nesmı´me take´ opomenout matematickou analy´zu, z nı´zˇ byly uka´zańy jen za´klady. Pozornost by si jisteˇ zaslouzˇily i odhady chyb. Martin Hasal

59

6

Literatura

[1] AXSELSSON, O. Iterative solution methods. Cambridge University Press, Cambridge, 1994. [2] BLAHETA, R. Matematicke´ modelovańı´ a metoda konecˇnyćh prvku˚. Vysoka´ sˇkola bańˇska´ - Technicka´ univerzita Ostrava a Za´padocˇeska´ univerzita v Plzni, 2011. [3] BOUCHALA, J. Funkciona´lnı´ analy´za. Pozna´mky z prˇedna´sˇek. [4] BOUCHALA, J. Nelinea´rnı´ funkciona´lnı´ analy´za. Pozna´mky z prˇedna´sˇek. [5] BOUCHALA, J. Variacˇnı´ metody. Vysoka´ sˇkola bańˇska´ - Technicka´ univerzita Ostrava a Za´padocˇeska´ univerzita v Plzni, 2010. [6] BRAMBLE, J.H. - PASCIAK, J.E. A preconditioning technique for indefinite systems resulting from mixed approximations of eliptic problems. Math. Comp.,50 (1988). Dostupne´ z . [7] BRDICˇKA, M. - SAMEK, L. - SOPKO, B. Mechanika kontinua. Academia, Praha, 2005. ISBN 80-200-1344-X (3. revidovane´ vydańı´). [8] CHORIN, J. A. - MARSDEN, J. E. A Mathematical Introduction to Fluid Mechanics. Marsden, Springer, 1993. ISBN 0-387-97918-2 (third editon). [9] Conservation of momentum[online].[cit. 2011-05-01]. Dostupne´ .

z

[10] DEVLIN, K. Proble´my pro trˇetı´ tisıćiletı´, Sedm nejveˇtsˇ´ıch nevyrˇesˇenyćh ota´zek matematiky. Nakladatelstvı´ Dokorˇań a Argo, 2005. [11] DOERING, CH. R. - GIBBON, J. D. Applied Analysis of the Navier-Stokes equations. Cambridge texts in applied mathematics, Cambridge university press, 1995. ISBN 0-521-44557-4. ´ BKOVA ´ , S. a kol. Mechanika tekutin. VSˇB-TU Ostrava, 2007. [12] DRA [13] ELMAN, H. C. - SILVESTER, D. J. - WATHEN, A. J. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Numerical Mathematics and Scientific Computation, Oxford University Press, New York, 1994. [14] GLOWINSKI, R. Numerical Methods for Fluids (Part 3). Elsevier, 2003. HANDBOOK of NUMERICAL ANALYSIS. ISBN 0-444-51224-1. [15] GREENOUGH, C. - ROBINSON, K. The Finite Element Library[online].c2000,[cit. 2011-05-01]. Dostupne´ z .

60

[16] KUZMIN, D. A Guide to Numerical Methods for Transport Equations[online].c2010,[cit. 2011-05-01]. Dostupne´ z . [17] LIESEN, J. - PARLETT, B. N. On nonsymmetric saddle point matrices that allow conjugate gradient iterations. Numerische Mathematik manuskript No. [18] MACUR, J. Dynamicke´ syste´my[online].[cit. 2011-05-01]. Dostupne´ z . [19] MEEK, J. L. Isoparametric Quadrilateral and Hexahedron Elements[online].[cit. 2011-05-01]. Dostupne´ z . [20] NOSKIEVICˇ, J. a kol. Mechanika tekutin. SNTL - Nakladatelstvı´ technicke´ literatury, Praha, 1987. ´, [21] POKORNY M. Navier-Stokesovy rovnice. 2008. Dostupne´ .

z

[22] RANNACHER, R. Finite Element Methods for the Incompressible Navier-Stokes Equations.Institute of Applied Mathematics University of Heidelberg 1999.INF 293/294, D-69120. Dostupne´ z . [23] REKTORYS, K.Variacˇnı´ metody v inzˇeny´rskyćh proble´mech a v proble´mech matematicke´ fyziky. Academia Praha 1994. ISBN 80-200-0714-8 (2. cˇeske´ vydańı´). [24] SPURK, J. H. - AKSEN, N. Fluid Mechanics. Springer, 2008. ISBN 978-3-540-73536-6 (second editon). [25] SˇI´STEK, J. Stabilization of finite element method for solving incompressible viscous flows. Praha, 2004. Diplomova´ praće na Strojnı´ fakulteˇ Cˇeske´ho vysoke´ho ucˇenı´ technicke´ho. Vedoucı´ diplomove´ praće Prof. RNDr. Pavel Burda CSc. [26] TEMAM, R. Navier-Stokes equations theory and numerical analysis.NorthHolland Publishing Company, 1997.ISBN: 07204 2840 8. Dostupne´ z <math.univ-lille1.fr/˜calgaro/ENS/biblio-m2ma/T.pdf>. [27] THOMASSET, F. Implementation of Finite Element Methods for Navier-Stokes equations. Springer series in computational physics, Springer-Verlag New York Inc., 1981.

61

A

Apendix

A.1 Opera´tory Uvedeme si neˇktere´ opera´tory pro dvourozmeˇrny´ prostor, ktere´ jsou vyuzˇ´ıvańy v te´to praći a nebyli uvedeny v prˇedchozıćh cˇa´stech u · v = u x v x + uy v y

∇u =

"

∇·u=

∂ux ∂x ∂uy ∂x

∂ux ∂y ∂uy ∂y

#

∂ux ∂uy + ∂x ∂y

∂uy ∂vy ∂ux ∂vx ∂ux ∂vx ∂uy ∂vy + + + ∂x ∂x ∂y ∂y ∂x ∂x ∂y ∂y # " ∂2u 2 x + ∂∂yu2x T ∂x2 ∆u = ∇ · (∇u) = ∂ 2 uy ∂2u + ∂y2y ∂x2 " # x x ux ∂u + uy ∂u ∂x ∂y (u · ∇)u = ∂u ∂u ux ∂xy + uy ∂yy

∇u : ∇u =

(A.1)

(A.2) (A.3) (A.4)

(A.5)

(A.6)

A.2 Vlastnosti opera´toru Zde budou uvedeny jen za´kladnı´ vlastnosti ktere´ vyuzˇ´ıva´me v te´to praći. Vıće v [23], kde lze nale´zt vsˇechny pojmy. Definice A.1 Opera´tor se nazy´va´ omezeny´ (ohranicˇeny´), ve sve´m definicˇnı´m oboru DA , lze-li najı´t takove´ cˇ´ıslo K ≥ 0, zˇe pro vsˇechna u ∈ DA platı´ kAuk ≤ Kkuk. Nejmensˇ´ı30 z cˇ´ısel K, jezˇ splnˇuje vztah (A.7) nazy´va´me normou opera´toru A. Zrˇejmeˇ platı´ kAuk ≤ kAkkuk.

(A.7)

(A.8)

Definice A.2 Necht’ DA je linea´l husty´ v H (Hilbertoveˇ prostoru). Opera´tor A linea´rnı´ v DA se nazy´va´ symetricky´ v DA , jestlizˇe pro kazˇdou dvojici prvku˚ u, v z DA platı´ (Au, v) = (u, Av). 30

Lze doka´zat, zˇe toto nejmensˇ´ı cˇ´ıslo existuje.

(A.9)

62

Definice A.3 Opera´tor A se nazy´va´ pozitivnı´ ve sve´m definicˇnı´m oboru DA , je-li symetricky´31 a platı´-li pro vsˇechna u ∈ DA (Au, u) ≥ 0, (A.10) prˇicˇemzˇ (Au, u) = 0 ⇒ u = 0 v DA .

(A.11)

Jestlizˇe mimoto existuje konstanta C > 0 takova´, zˇe pro vsˇechna u ∈ DA je dokonce (Au, u) ≥ C 2 kuk2 ,

(A.12)

nazy´va´ se opera´tor A pozitivneˇ definitnı´ v DA . Veˇta A.1 Greenova veˇta v obecne´ formeˇ, prˇes oblast G. Z Z Z ∂c ∂b b dx = b c vl dS − c dx, ∂x ∂x l l G Γ G

(A.13)

kde b ∈ W21 , c ∈ W21 a vl je l = 1, . . . , n prˇedstavuje sourˇadnice jednotkove´ho vektoru vneˇjsˇ´ı norma´ly.

A.3 Veˇty Veˇta A.2 Poincare–Friedrichsova nerovnost. Prˇedpokla´dejme, zˇe Ω ⊂ R2 je obsazˇena ve cˇtverci R de´lky L (a v prˇ´ıpadeˇ, zˇe ∂ΩN dS 6= 0, tj. zajisˇteˇnı´ dostatecˇneˇ hladke´ hranice). Vzhledem k tomu, R zˇe ∂ΩD dS 6= 0, dosta´va´me, zˇe kvk ≤ Lk∇vk (A.14) pro vsˇechna v ∈ HE1 0 . L nazveme Poinrare´ho konstantou. Du˚kaz, lze zı´skat z Friedrichsovy nerovnosti viz [23].

A.4 Prostory Definice A.4 Metricky´m prostorem L2 na oblasti Ω mnozˇinu funkcı´, jejichzˇ integra´l(ve smyslu Lebesgueova integra´lu) z druhe´ mocniny funkce v oblasti A je konecˇny´, tj. Z L2 (Ω) := u : Ω → R| u2 < ∞ . (A.15) Ω

A prˇ´ıslusˇna´

L2

norma kuk :=

31

Z

u Ω

2

1/2

(A.16)

Lze uka´zat, zˇe v prˇ´ıpadeˇ komplexnı´ho Hilbertova prostoru vyply´va´ symetricˇnost opera´toru A prˇ´ımo z pozˇadavku (A.10). U rea´lne´ho Hilbertova prostoru je trˇeba pozˇadavek symetricˇnosti zvla´sˇteˇ vyslovit.

63

Definice A.5 Sobolevovy´m prostorem H 1 na oblasti Ω nazveme vsˇechny funkce pro neˇzˇ platı´ ∂u ∂u 1 2 H (Ω) := u : Ω → R|u, , ∈ L (Ω) . (A.17) ∂x ∂y A prˇ´ıslusˇna´ H 1 norma kuk :=

Z

2

u + Ω

!1/2 Z X N ∂u 2 Ω

1

∂xi

(A.18)

Navierových-Stokesových rovnic. Numerical Solution of Navier-Stokes equations using the finite element method

Recommend Documents