Počítačová mechanika tekutin - cvut.czmarian.fsik.cvut.cz/~furst/PMT/pmt.pdf · 2014-03-03 ·...

Výchova studentů pro aplikace řešené na vý-konných počítačích

České vysoké učení technické v PrazeFakulta strojní

Počítačová mechanika tekutin

Jiří Fürst

TENTO PROJEKT JE SPOLUFINANCOVÁN EVROPSKÝMSOCIÁLNÍM FONDEM A STÁTNÍM ROZPOČTEM ČESKÉREPUBLIKY A ROZPOČTEM HLAVNÍHO MĚSTA PRAHY

Preliminary version – 3. března 2014

Obsah

1 Eulerovy rovnice 21.1 Linearizace Eulerových rovnic . . . . . . . . . . . . . . . . . . 21.2 Okrajové podmínky pro Eulerovy rovnice . . . . . . . . . . . 3

1.2.1 Podmínky pro linearizované rovnice . . . . . . . . . . 3

2 Numerické řešení Eulerových rovnic 52.1 Rusanovovo schéma . . . . . . . . . . . . . . . . . . . . . . . . 52.2 VanLeerovo schéma . . . . . . . . . . . . . . . . . . . . . . . . 6

3 Potenciální proudění 73.1 Stacionární potenciální proudění ve 2D . . . . . . . . . . . . . 73.2 Linerizovaná rovnice pro poruchový potenciál . . . . . . . . . 8

3.2.1 Příklad: subsonické proudění GAMM kanálem . . . . 9

4 Tvorba sítí 134.1 Strukturované sítě . . . . . . . . . . . . . . . . . . . . . . . . 13

4.1.1 Algebraické sítě . . . . . . . . . . . . . . . . . . . . . . 134.1.2 Eliptické sítě . . . . . . . . . . . . . . . . . . . . . . . 164.1.3 Hyperbolické sítě . . . . . . . . . . . . . . . . . . . . . 17

4.2 Nestrukturované sítě . . . . . . . . . . . . . . . . . . . . . . . 194.2.1 Delaunayovská triangulace oblasti . . . . . . . . . . . 19

4.3 Hybridní sítě . . . . . . . . . . . . . . . . . . . . . . . . . . . 234.4 Adaptivní strukturované sítě . . . . . . . . . . . . . . . . . . 244.5 Tvorba sítí ve 3D . . . . . . . . . . . . . . . . . . . . . . . . . 24

Poznámka na úvodTento text slouží především jako poznámky pro výuku předmětu Počítačovámechanika tekutin na FS ČVUT. Jedná se pouze o určité doplnění existu-jící literatury (např. série skript Numerická simulace proudění vydávaná naÚTM FS ČVUT). Některé kapitoly jsou zatím spíše osnovou přednášky a narozšíření se stále pracuje. Soubor pmt.pdf obsahuje tento dokument v PDFformátu.

1


Kapitola 1

Eulerovy rovnice

Eulerovy rovnice popisují prodění nevazké stlačitelné tekutiny. Vzniknou zúplného systému Navierových-Stokesových rovnic zanedbáním členů popisu-jících vazké efekty a vedení tepla. Z tohoto důvodu je lze použít pouze tam,kde tyto efekty nehrají rozhodující roli. Celý systém Eulerových rovnic lzezapsat v tzv. konzervativním tvaru jako

Wt + F (W )x +G(W )y +H(W )z = 0, (1.1)

kde

W = [ρ, ρu, ρv, ρw, e]T , (1.2)

F (W ) =[ρu, ρu2 + p, ρuv, ρuw, (e+ p)u

]T, (1.3)

G(W ) =[ρv, ρuv, ρv2 + p, ρvw, (e+ p)v

]T, (1.4)

H(W ) =[ρw, ρuw, ρvw, ρw2 + p, (e+ p)w

]T. (1.5)

Tlak p je dán stavovou rovnicí, která pro případ ideálního plynu vede navztah

p = (γ − 1)[e− 1

2ρ(u2 + v2 + w2)]. (1.6)

V dalším textu se budeme zabývat pouze jedním problémem a to vhod-nou volbou okrajových podmínek pro Eulerovy rovnice. Ostatní problémy(teorie Eulerových rovnic a numerické metody jejich řešení) jsou popsány vevýše uvedené sérii skript.

1.1 Linearizace Eulerových rovnicPřeveďme nejprve Eulerovy rovnice pro jednorozměrný případ do primitiv-ních proměnných U = [ρ, u, p]T . Je tedy

Ut +AUx = 0, (1.7)

2


kde

A =

u ρ 00 u 1/ρ0 γp u

. (1.8)

Provedeme-li linearizaci kolem konstantního stavu ρ, u, p a dostáváme

Ut + AUx = 0, (1.9)

kde

A =

u ρ 00 u 1/ρ0 γp u

. (1.10)

Vlastní čísla matice A jsou u − a, u a u + a, kde a2 = γp/ρ je kvadrátrychlosti zvuku.

Vlastní vektory zapsané po sloupcích do matice R jsou

R =

ρ/a 1 ρ/a−1 0 1ρ/a 0 ρ/a

(1.11)

a její inverze je

R−1 =

0 −1/2 1/(2ρa)1 0 −1/a2

0 1/2 1/(2ρa)

. (1.12)

1.2 Okrajové podmínky pro Eulerovy rovniceV následující sekci si ukážeme jeden z přístupů pro konstrukci okrajovýchpodmínek založený na linearizovaných Eulerových rovnicích. K pochopenítěchto podmínek je třeba znát principy řešení lineárních hyperbolickýchsystémů PDR prvního řádu. Tato problematika je opět popsána ve výšezmíněných skriptech.

1.2.1 Podmínky pro linearizované rovnice

Při popisu okrajových podmínek se omezíme nejprve na jednorozměrný pří-pad s rychlostí u ≥ 0 a vyjdeme z linearizovaných Eulerových rovnic pro pri-mitivní proměnné. Z nich můžeme snadno určit charakteristické proměnnéjako

V = R−1U = [p/(ρa)− u, ρ− p/a2, p/(ρa) + u]T . (1.13)

Podmínka na vstupu, subsonický případ

U subsonického vstupu pouze jedna charakteristika vystupuje z proudovéhopole, tedy v1 = const.. Je tedy třeba zadat dvě veličiny, které nám spolu sextrapolovanou hodnotou v1 umožní určit všechny složky toku hranicí.

3


Zadaný celkový tlak a teplota V tomto případě jsou zadány hodnotycelkového tlaku p0 a celkové T0. Z izentropických vztahů

T0 = (1 + γ − 12 M2)Tb, (1.14)

p0 = (1 + γ − 12 M2)

γγ−1 pb (1.15)

určíme hodnoty tlaku a teploty na vstupu. Velikost rychlosti pak určíme zcharakteristické podmínky

ρa(ub − u1)− (pb − p1) = 0. (1.16)

Podmínka na pevné stěně

Na pevné stěně platí, že u = 0. To znamená, že pouze jedna charakteris-tika vstupuje do oblasti a měli bychom tedy předepsat jednu veličinu. Provýpočet toku hraniční stěnou je zapotřebí znát pouze hodnotu tlaku. Před-pokládejme, že stěna je umístěná na pravém okraji oblasti. Pak hodnotatlaku na stěně se určí z podmínky v3 = const. a ub = 0, tedy

pb = pN + ρauN , (1.17)

kde uN je rychlost směrem ke stěně v poslední buňce.

Podmínka na výstupu, subsonický případ

Pro případ subsonického výstupu se zadavá jediná veličina. Jedna z možnostíje zadat hodnotu tlaku pb. Zbylé dvě veličiny dopočítáme z charakteristic-kých proměnných. Tedy

(pb − pN )− a2(ρb − ρN ) = 0, (1.18)(pb − pN ) + ρa(ub − uN ) = 0, (1.19)

neboli

ρb = ρN + (pb − pN )/a2, (1.20)ub = uN − (pb − pN )/(ρa). (1.21)

4


Kapitola 2

Numerické řešení Eulerovýchrovnic

V této části se budeme zabývat problémem numerického řešení Eulerovýchrovnic. Budeme vycházet z metody konečných objemů (viz např. [3], [2], [4]).

Budeme uvažovat nejprve jednorozměrný systém Eulerových rovnic vetvaru

Wt + F (W )x = 0, (2.1)

kde W = [ρ, ρu, e]T a F (W ) = [ρu, ρu2 + p, (e+ p)u].V tomto případě budeme uvažovat numerické schéma ve tvaru

Wn+1i = Wn

i −∆t∆xi

[H(Wn

i ,Wni+1)−H(Wn

i−1,Wni )], (2.2)

kde H(WR,WL) je tzv. numerický tok.Poté přejdeme ke dvourozměrnému systému ve tvaru

Wt + F (W )x +G(W )y = 0, (2.3)

kde W = [ρ, ρu, ρv, e]T , F (W ) = [ρu, ρu2 + p, ρuv, (e + p)u] a G(W ) =[ρv, ρuv, ρv2 + p, (e+ p)v].

Numerické schéma bude v tomto případě

Wn+1i = Wn

i −∆t|Ωi|

∑j

H(Wni ,W

nj , ~Sij). (2.4)

2.1 Rusanovovo schémaJedno z nejjednodušších schémat. Numerický tok pro jednorozměrný případje

H(WR,WL) = 12 [F (WL) +W (WR)− (|u|+ a) (WR −WL)] . (2.5)

5


2.2 VanLeerovo schémaVanLeerovo schéma patří mezi tzv. flux-splitting metody, kdy je tok F (W )rozdělen na dvě části F+(W ) a F−(W ) s kladnými resp. zápornými vlast-ními čísly příslušných Jakobiho matic a numerický tok je poté počítán jakoH(WL,WR) = F+(WL) + F−(WR).

Pro případ Eulerových rovnic je toto rozdělení navíc provedeno zvlášťpro tzv. konvektivní část toku F c(W ) = u[ρ, ρu, e+ p]T a pro tlakovou částF p(W ) = [0, p, 0]T .

Označme Ψ = [ρ, ρu, e+ p]T . Pak konvektivní část toku pro VanLeerovoschéma je

Hc(WL,WR) =MLaLΨL +MRaRΨR, (2.6)

a tlaková část toku je

Hp(WL,WR) = [0,PL, 0]T + [0,PR, 0]T , (2.7)

a celkový tok je tedy

H(WL,WR) = Hc(WL,WR) +Hp(WL,WR). (2.8)

VeličinyM a P jsou definovány vztahy

ML =

0, pro ML ≤ −114(1 +ML)2, pro − 1 < ML < 1ML, pro 1 ≥ML

(2.9)

PL =

0, pro ML ≤ −1pL4 (1 +ML)2(2−ML), pro − 1 < ML < 1pL, pro 1 ≥ML

(2.10)

a

MR =

MR, pro MR ≤ −1−1

4(1−MR)2, pro − 1 < MR < 10, pro 1 ≥MR

(2.11)

PR =

pR, pro MR ≤ −1pR4 (1−MR)2(2 +MR), pro − 1 < MR < 1

0, pro 1 ≥MR

(2.12)

kde ML/R = uL/RaL/R

.

6


Kapitola 3

Potenciální proudění

Systémy Eulerových resp. Navierových-Stokesových rovnic popisují poměrněširokou třídu tekutin, jejich řešení však může být komplikované a tím pádemi časově velmi náročné. Proto lze v některých případech použít zjednodu-šené modely. Jedním z takovýchto zjednodušených modelů je model tzv.potenciálního proudění. Předpokládejme, že:

• proudění je nevazké,

• nevířivé,

• a izentropické.

Pro jednoduchost dále předpokládejme, že řešíme stacionární problémve 2D (tyto předpoklady nejsou nutné, nám však usnadní zápis rovnic, pro3D nestacionární případ odkazujeme čtenáře např. na [1]).

3.1 Stacionární potenciální proudění ve 2DJe-li proudění stacionární a nevazké, je popsáno soustavou Eulerových rovnic

(ρu)x + (ρv)y = 0, (3.1)(ρu2 + p)x + (ρuv)y = 0, (3.2)(ρuv)x + (ρv2 + p)y = 0. (3.3)

Rozderivujeme rovnici (3.1) a dostaneme

ρxu+ ρux + ρyv + ρvy = 0. (3.4)

Je-li proudění izentropické, je z definice rychlosti zvuku c2 = (∂p∂ρ)s gradienthustoty ρx = 1

c2 px a ρy = 1c2 py. Tedy

1c2 (pxu+ pyv) + ρ(ux + vy) = 0. (3.5)

7


Z rovnice (3.2) vyjádříme px jako

px = −(ρu2)x − (ρuv)y = −(ρu)xu− ρuux − (ρv)yu− ρvuy == −ρuux − ρvuy. (3.6)

Podobně z (3.3) je

py = −(ρuv)x − (ρv2)y = −(ρu)xv − ρuvx − (ρv)yv − ρvvy == −ρuvx − ρvvy. (3.7)

Dosadíme derivace tlaku do rovnice (3.5), výslednou rovnici vydělímehustotou a vynásobíme c2 a dostáváme

− u2ux − uvuy − uvvx − v2vy + c2(ux + vy) = 0. (3.8)Je-li proudění nevířivé, existuje potenciál rychlosti, tj. u = Φx a v = Φy.

Dosadíme tyto vztahy do poslední rovnice a dostaneme(c2 − (Φx)2

)Φxx − 2ΦxΦyΦxy +

(c2 − (Φy)2

)Φyy = 0. (3.9)

Tato rovnice se nazývá rovnice potenciálního proudění nebo též potenciálnírovnice. Jedná se o nelineární skalární rovnici druhého řádu. Je-li

• u2 + v2 < c2, jedná se o rovnici eliptického typu,

• u2 + v2 > c2, jedná se o rovnici hyperbolického typu,

• u2 + v2 = c2, jedná se o rovnici parabolického typu.

3.2 Linerizovaná rovnice pro poruchový potenciálDalším zjednodušením předchozí rovnice je linearizace potenciální rovnicepomocí metody malých poruch. Předpokládejme, že do volného proudu okonstantní rychlosti (u, v) = (U∞, 0) vložíme těleso o velmi malé tloušťce.To způsobí v rovnoměrném proudu určitou poruchu. Rychlost tedy bude

u = U∞ + u′, (3.10)v = v′, (3.11)

kde u′, v′ jsou poruchy rychlosti. Nechť φ je potenciál těchto poruch (tj.u′ = φx a v′ = φy). Potom

Φ = U∞x+ φ. (3.12)

Dosadíme tento potenciál do rovnice (3.9) a zanedbáme všechny členy ob-sahující vyšší mocniny φ. Dále z izentropických vztahů pro rychlost zvukuc2

0 = c2 + γ−12 (u2 + v2) dostaneme

c2 = c2∞ +O(φ), (3.13)

8


a dostáváme lineární rovnici pro potenciál malých poruch

(1−M2∞)φxx + φyy = 0. (3.14)

3.2.1 Příklad: subsonické proudění GAMM kanálem

Jako jednoduchý příklad řešení potenciálního proudění si uveďme subsonicképroudění tzv. GAMM kanálem. Jedná se o kanál délky 3 a výšky 1. Na jehodolní stěně je mezi body [1, 0] a [2, 0] překážka ve tvaru kruhového obloukuprocházejícího bodem [1.5, 0.1].

Na vstupu budeme uvažovat proudění o parametrech M∞ = 0.5, u =U∞ = 1 a v = 0. Na horní a dolní stěně uvažujeme podmínku neprostupnosti,tj. (u, v) ·~n = 0. Tyto podmínky transformujeme na okrajové podmínky proporuchový potenciál φ:Vstup je-li u = U∞, je φx = 0,

Stěny podmínka neprostupnosti je 0 = (u, v) · (n1, n2) = (U∞ + φx, φy) ·(n1, n2) a tedy ∂φ/∂n = −U∞n1,

Výstup Ačkoliv jsme při popisu úlohy nezmiňovali žádnou podmínku navýstupu, pro korektní formulaci je třeba i zde zadat vhodnou okrajovoupodmínku. Předepíšeme zde např. podmínku v = 0, co6 lze přepsatjako φ = 0.

Celý výpočet ještě zjednodušíme tak, že využijeme předpokladu, že pře-kážka je velmi tenka a aproximujeme okrajové podmínky na dolní stěněpouze s prvním řádem přesnosti (vzhledem k tloušťce překážky). Tak budemožno provést výpočet na jednoduché kartézské síti a tvar překážky budememodelovat pouze zvolenou okrajovou podmínkou.

Program Potential/gamm-simple.cpp implemetnuje řešení tohoto pro-blému pomocí Gaussovy-Seidelovy iterační metody. Výsledek (tj. hodnotyporuchového potenciálu) jsou uloženy v souboru Potential/phi.plt. Tentosoubor lze zobrazit pomocí komečního programu Tecplot nebo pomocí volnědostupného programu Visit .

Obrázek 3.1 ukazuje uživatelské rozhraní programu Visit. Po jeho spuš-tění najdeme v levém okně soubor phi.plt a stiskneme open. Pomocí nabídkyPlots/Mesh zobrazíme výpočetní síť (viz obr. 3.2). Dále můžeme zobrazitrozložení vypočtené hodnoty poruchového potenciálu pomocí volby Plots/-Pseudocolor (obr. 3.3). Pokud bychom chtěli zobrazit například rychlosti,muzeli bychom je nejprve z poruchového potenciálu dopočítat. To můžemeudělat pomocí okna Controls/Expression. Obrázek 3.4 ukazuje postup přivýpočtu rychlosti jako součtu konstantního vektoru (U∞, 0) = (1, 0) a gra-dientu potenciálu a dále výpočet velikosti rychlosti. Pomocí nabídky Plots/-Vector můžeme zobrazit vektory rychlosti (obr. 3.5). Vybereme-li v řídícímokně graf zobrazující potenciál, můžeme pomocí možnosti Variables změnitzobrazení potenciálu na zobrazení velikosti rychlosti (obr. 3.6).

9


Potential/gamm-simple.cpp

Potential/phi.plt

Obrázek 3.1: Uživatelské rozhraní programu Visit

Obrázek 3.2: Zobrazení sítě

10


Obrázek 3.3: Zobrazení poruchového potenciálu

Obrázek 3.4: Výpočet rychlosti a její velikosti

11


Obrázek 3.5: Zobrazení vektorů rychlosti

Obrázek 3.6: Zobrazení velikosti rychlosti

12


Kapitola 4

Tvorba sítí

V této kapitole se budeme věnovat některým metodám tvorby sítí. Pokusímese vysvětlit principy těchto metod na relativně jednoduchých dvourozměr-ných případech. Problematika tvorby sítí ve trojrozměrném prostoru přisložitější geometrii se již rozsahem značně vymyká z možností tohoto kurzua případnému čtenáři doporučujeme například http://www.hpc.msstate.edu/publications/gridbook/ nebo [5].

Algoritmy popisované v této publikaci jsou doplněny demonstračnímiprogramy odladěnými v jazyce Python s rozšířením pro numerické výpočtyNumeric a pro zobrazovaní sítí používaji program Gnuplot. Tento je pronázorné animace ovladán z jazyka python pomocí modulu Gnuplot.py.

Pro správné fungování je zapotřebí modul Grids/grid-base.py.Algoritmy jsou testovány na následujících oblastech:

• jednotkový čtverec Grids/ctverec.py,

• kanál s půlkruhohou překážkou Grids/channel.py,

• kanál s půlkruhohou překážkou se zjemněním Grids/channel-ns.py,

• GAMM kanál Grids/gamm.py,

• GAMM kanál se zjemněním Grids/gamm-ns.py.

4.1 Strukturované sítě

4.1.1 Algebraické sítě

Jednou z nejjednodušších metod tvorby strukturovaných sítí je tzv. alge-braická metoda. Ta vychází z přímé geometrické konstrukce sítě ze zadanýchbodů na hranici oblasti. V následujícím textu se budeme pokoušet vytvo-řit dvourozměrnou strukturovanou síť o N timesM buňkách. Vrcholy sítěbudeme označovat xi,j , kde i = 0, .., N a j = 0, ..,M .

13


http://www.hpc.msstate.edu/publications/gridbook/


www.python.org

http://numeric.scipy.org/

www.gnuplot.info

gnuplot-py.sf.net

Grids/grid-base.py

Grids/ctverec.py

Grids/channel.py

Grids/channel-ns.py

Grids/gamm.py

Grids/gamm-ns.py

Nejprvme si představme, že máme zadány pouze body na dolní a horníhranici oblasti, tj. jsou zadány body xi,0 a xi,M pro i = 0, ..., N . Body uvnitřoblasti pak můžeme získat vhodnou interpolací těchto hraničních bodů. Vnejjednodušším případě můžeme použít lineární interpolaci a máme tedy

xi,j = (1− ηj)xi,0 + ηjxi,M , (4.1)

kde1 ηj = j/M a indexy jsou i = 0, .., N a j = 1, ..M − 1. Tato jednodu-chá metoda vytváří síť bez jakéhokoliv zjemnění ve směru indexu j. Navíclevý a pravý okraj sítě oblasti musí být tvořen úsečkou. Program Grid-s/alg_simple.py ukazuje ukazuje implementaci tohoto algoritmu. Pro jehopoužití nejprve vytvořte jedním z výše uvedených programů popis hraniceoblasti a potom spusťte tento program.

Určitým vylepšením této metody je zavedení vhodnějších interpolačníchmetod. Lineární interpolaci v (4.1) je možné nahradit interpolací danoufunkcemi βk(η), tj. např.

xi,j = β0(ηj)xi,0 + β1(ηj)xi,M , (4.2)

kde β0(η) +β1(η) = 1 a β0 je klesající funkce, pro níž β0(0) = 1 a β0(1) = 0.Velikost derivace funkce β0 určuje výšku buněk v dané vrstvě. Čím je velikostderivace větší, tím větší je i výška buněk.

Dalším rozšířením je interpolace pomocí Lagrangeova polynomu. Mějměv tomto případě zadány nejen body ve vrstvách 0 a M , ale ještě v několikadalších vrstvách uvnitř sítě. Označme tyto vrstvy včetně krajních jako jk,kde k = 1, ..,K. Potom Lagrangeovskou síť dostaneme jako

xi,j =∑k

βk(ηj)xi,jk , (4.3)

kde βk(ηjl) = δkl.Velmi užitečná může též být interpolace Hermitovská. Zde máme v ně-

kolika vrstvách zadány nejen souřadnice uzlů, ale také první (nebo vyšší)derivace podle η (neboli směry síťových čar). V nejjednodušším, přesto všakvelmi užitečném, případě jsou zadány body x a směry síťových čar spolu svelikostí buněk xξ na hranici. Potom

xi,j = β00(ηj)xi,0 + β1

0(ηj)(xη)i,0 + β01(ηj)xi,M + β1

1(ηj)(xη)i,M . (4.4)

Funkce βlk jsou přitom zvoleny tak, že

β00(0) = 1, (β0

0)′(0) = 0, β00(1) = 0, (β0

0)′(1) = 0,β1

0(0) = 0, (β10)′(0) = 1, β1

0(1) = 0, (β10)′(1) = 0,

β01(0) = 0, (β0

1)′(0) = 0, β01(1) = 1, (β0

1)′(1) = 0,β1

1(0) = 0, (β11)′(0) = 0, β1

1(1) = 0, (β11)′(1) = 1.

(4.5)

1Pozor na celočíselné dělení!

14


Grids/alg_simple.py

Grids/alg_simple.py

Všechny výše uvedené interpolace však nebraly v úvahu tvar levé a pravéstrany oblasti. Předpokládejme nyní, že jsou zadány body na všeh čtyřechstranách oblasti. Pak lze body uvnitř sítě získat tzv. metodou transfinitníinterpolace. Ztotožněme zadané body xi,j na hranicích oblasti s hodnotamifunkce x(ξi, ηj). Nechť

U(ξ, η) = α0(ξ)x(0, η) + α1(ξ)x(1, η), (4.6)

je jednorozměrná interpolace ve směru indexu i, funkce α0,1 mohou být např.α0(ξ) = 1− ξ a α1(ξ) = ξ. Podobně

V (ξ, η) = β0(η)x(ξ, 0) + β1(ξ)x(ξ, 1), (4.7)

je interpolace ve směru indexu j. K tomu, aby součet U + V byl vhodnouinterpolací, je třeba odečíst korekci

UV (ξ, η) = α0(ξ)β0(η)x(0, 0) + α1(ξ)β0(η)x(1, 0)++ α1(ξ)β1(η)x(1, 1) + α0(ξ)β1(η)x(0, 1), (4.8)

tedyx(ξ, η) = U(ξ, η) + V (ξ, η)− UV (ξ, η). (4.9)

Obrázek 4.1 ukazuje síť v GAMM kanále získané pomocí výše uvedenémetody pomocí programu Grids/alg_tfi_lin.py.

(a) Celá oblast (b) Detail

Obrázek 4.1: Algebraická síť v GAMM kanále se zjemněním u stěn vytvořenálineární TFI interpolací.

Tato jednoduchá metoda však selhává pro případ, kdy je požadovánozjemnění sítě v blízkosti hranice. V tomto případě lze postupovat následu-jícím způsobem: pro body na hranici nejprve vypočteme délky křivek a tonásledujícím způsobem:

• pro levou hranici je u0,j = 0 a v0,0 = 0 a v0,j = v0,j−1 + ||x0,j −x0,j−1||/

∑Ml=1 ||x0,l − x0,l−1||,

• pro pravou hranici je uN,j = 1 a vN,0 = 0 a vN,j = vN,j−1 + ||xN,j −xN,j−1||/

∑Ml=1 ||xN,l − xN,l−1||,

15


Grids/alg_tfi_lin.py

• pro dolní hranici je vi,0 = 0 a u0,0 = 0 a ui,0 = ui−1,0 + ||xi,0 −xi−1,0||/

∑Nl=1 ||xl,0 − xl−1,0||,

• pro horní hranici je vi,M = 1 a u0,M = 0 a ui,M = ui−1,M + ||xi,M −xi−1,M ||/

∑Nl=1 ||xl,M − xl−1,M ||.

Potom se použije výše uvedená transfinitní interpolace (4.9) pro výpočethodnot ui,j a vi,j z hraničních hodnot a hodnoty parametrů u a v se použijína místě ξ a η při interpolaci uzlů sítě, viz program Grids/alg_tfi_bsg.py.

Pomocí tohoto algoritmu byla získána například síť na obrázku 4.2.


Obrázek 4.2: Algebraická síť v GAMM kanále se zjemněním u stěn vytvořenálineární TFI interpolací s uvažováním rozdělení na hranicích.

4.1.2 Eliptické sítě

Eliptická síť se získá řešením Poissonovy rovnice

∆ξ = P (x, y), (4.10)∆η = Q(x, y), (4.11)

se okrajovými podmínkami určenými daným rozdělením bodů na hranicíchoblasti. Síťové čáry zde odpovídají izočarám veličin ξ a η.

Pro účely numerického řešení je vhodnější tuto rovnici transformovat dokřivočarých souřadnic ξ, η jako (viz [5])

g22 (xξξ + Pxξ)− 2g12xξη + g11 (xηη +Qxη) = 0, (4.12)

kde

g11 = xξ · xξ = x2ξ + y2

ξ , (4.13)g12 = xξ · xη = xξxη + yξyη, (4.14)g22 = xη · xη = x2

η + y2η. (4.15)

Význam členů P aQ je podrobněji diskutován v http://www.hpc.msstate.edu/publications/gridbook/. Obecně řečeno kladné hodnoty P posou-vají síťové čáry ve směru růstu ξ (doprava), záporné naopak. Podobně Q

16


Grids/alg_tfi_bsg.py



posouvá síťové čáry ve směru souřadnice η (nahoru nebo dolů). Vhodnouvolbou těchto zdrojových členů lze dosáhnout například zjemnění sítě nebokolmost síťových čar v blízkosti hranice.

Obrázek 4.3 ukazuje síť získanou řešením rovnice (4.12) bez uvažovánízdrojových členů P a Q pomocí programu Grids/eli_simple.py. Obrázek4.4 ukazuje síť získanou řešením problému se zdrojovými členy zaručujícímiortogonalitu a zachování velikosti buněk u hranice.


Obrázek 4.3: Eliptická síť v GAMM kanále se zjemněním u stěn, při tvorběsítě nebyly uvažovány zdrojové členy. Všimněte si, že síť vůbec nerespektujepožadované zjemnění u stěn.


Obrázek 4.4: Eliptická síť v GAMM kanále se zjemněním u stěn, zdrojovéčleny dle kapitoly 6.3.2 z [5]

4.1.3 Hyperbolické sítě

Další možností pro tvorbu sítí je formulovat požadavek ortogonality síťovýchčar jako

xξ · xη = xξxη + yξyη = 0. (4.16)Tato rovnice nám při zadaném směru xξ umožnuje určit směr vektoru xη. Jevšak třeba doplnit rovnici, pomocí které učíme délku tohoto vektoru (nebolivelikost kroku sítě ve směru kolmém na stěnu). Touto rovnicí může býtnapříklad

xξ × xη = xξyη − xηyξ = V (ξ, η), (4.17)

17


Grids/eli_simple.py

kde V je plocha buňky.Linearizujeme-li rovnice (4.16) a (4.17) kolem stavu x, tj. x = x + x′,

dostaneme

xξxη + yξ yη + xξx′η + yξy

′η + x′ξxη + y′ξ yη = 0,

xξ yη − xηyξ + xξy′η − xηy′ξ + x′ξ yη − x′ηyξ = V (ξ, η).

Vzhledem k tomu, že pro stav s pruhem platí xξxη + yξ yη = 0, můžeme tutorovnici přičíst k první rovnici posledně uvedeného systému. Dále využijemetoho, že x + x′ = x a dostáváme

xξxη + yξyη + xηxξ + yηyξ = 0. (4.18)

Podobně pro druhou rovnici soustavy využijeme vztah xξ yη−xηyξ = V (ξ, η)a dostáváme

xξyη − xηyξ + xξ yη − xηyξ = V (ξ, η) + V (ξ, η). (4.19)

Tyto dvě rovnice můžeme přepsat do vektorového tvaru jako

Bxη + Axξ = f , (4.20)

kde

A =[xη yηyη −xη

], B =

[xξ yξ−yξ xξ

], f =

[0

V + V

]. (4.21)

Tento systém můžeme upravit na tvar

xη + Cxξ = g, (4.22)

kde C = B−1A a g = B−1f . Přímým výpočtem lze zjistit, že

B−1 = 1x2ξ + y2

ξ

[xξ −yξyξ xξ

],

C = 1x2ξ + y2

ξ

[xξxη − yξ yη xξ yη + yξxηyξxη + xξ yη yξ yη − xξxη

],

g = V + V

x2ξ + y2

ξ

[−yξxξ

].

Matice C je symetrická a má tedy reálná vlastní čísla a systém (4.22) jetedy hyperbolický. To nám dává návod na sestavení numerické metody provýpočet hodnot x(ξ, η) pro η > 0 při zadané počáteční podmínce x(ξ, 0).

Vlastní čísla matice C jsou

λC = ±√

2x2ξ + y2

ξ

√x2ξ x

2η + y2

ξ y2η. (4.23)

18


Výpočet j + 1 vrstvy sítě ze zadané j-té vrstvy tedy odpovídá jednomukroku metody pro řešení hyperbolického systému (4.22). Pro jeho řešení lzevyužít buď explicitní nebo implicitní metodu. Obě metody vyžadují stabili-zaci pomocí nějakého druhu umělé vazkosti. Explicitní variantu centrálníhoschématu lze zapsat jako

xi,j+1 = xi,j −12Ci,j (xi+1,j − xi−1,j) + gi,j + εi,j (xi+1,j − 2xi,j + xi−1,j) .

(4.24)Tato metoda však vyžaduje splnění podmínky stability ∆η ≤ ∆ξ/ρC. Vnašem případě je ∆η = ∆ξ = 1, takže nejsme schopni volbou vhodného ∆ηdosáhnout stability metody. Naproti tomu implicitní metoda ve tvaru

... (4.25)

4.2 Nestrukturované sítěDalším typem sítí jsou tzv. nestrukturované sítě. Na rozdíl od strukturova-ných sítí je u těchto sítí pro určení sousedních vrcholů nebo buněk mít kdispozici informaci o topologii sítě.

Ve dvourozměrném případě se nejčastěji využívají sítě tvořené trojúhel-níkovými nebo čtyřúhelníkovými buňkami.

4.2.1 Delaunayovská triangulace oblasti

Jedna z metod tvorby nestrukturované trojúhelníkové sítě je tzv. Delauna-yovská triangulace. Ta je vytvořena tak, že uvnitř kruhu opsaného libovol-nému trojúhelníku sítě se nenachází žádný vrchol triangulace.

Algoritmus pro tvorbu triangulace byl popsán na přednášce. Zde si uká-žeme příklad použití volně dostupného programu triangle. Tento programpracuje v několika režimech. Pro naše účely bude nejužitečnější režim, kdyprogram automaticky vygeneruje síť uvnitř oblasti se zadanou hranicí. Projednoduchost předpokládejme, že oblast je jednoduše souvislá, tj. neobsa-huje žádné díry. V tomto případě připravíme vstupní formát typu *.poly.Ten obsahuje body a úseček na hranici oblasti a má následující strukturu:

• První řádek obsahuje: počet vrcholů ve vstupním souboru, dimenzeoblasti (musí být 2), počet atributů, počet hraničních značek (0 nebo1).

• Následují řádky obsahující: číslo vrcholu, jeho souřadnice x a y, atri-buty vrcholu (pokud byl nastaven počet atributů větší než 0), hraničníznačka (pokud byl nastaven počet začek 1).

• Další sekce souboru začíná počtem hraničních úseček spolu s počtemhraničních značek.

19


http://www.cs.cmu.edu/~quake/triangle.html

• Dále následují řádky obsahující číslo hraniční úsečky, indexy jejíchkrajních vrcholů a případne hraniční značku.

• Třetí sekce popisuje díry v oblasti a začína řádkem obsahujícím početděr (pro jednoduše souvislé oblasti je zde 0).

• Následující rádky obsahují číslo díry a souřadnice libovolného boduuvnitř díry.

• Další sekce jsou nepovinné a pro naše účely zatím zbytečné.

Vstupní soubor obsahuje mimo informací o pozici bodů také atributy a hra-niční značky. Atributy jsou libovolná čísla, která generátor sítě při síťováníinterpoluje na nově vzniklé body. Lze tak například nechat přepočítávatřešení při adaptaci sítě. My budeme zatím používat soubory bez atributů.

Naproti tomu hraniční značky pro nás budou velmi užitečné. Očíslujemesi ve vstupním souboru hranice oblasti a tato čísla použijeme jako hraničnízvonky. To nám později umožní rozpoznat, která hrana či bod náleží jakéhranici.

Síť ve čtverci Jako první příklad si ukážeme tvorbu sítě v oblasti tvaručtverce. Tu můžeme popsat vstupním souborem Grids/ctverec.poly s násle-dujícím obsahem:

# V souboru jsou 4 body bez atributu a hranicnich znacek4 2 0 01 0.0 0.02 1.0 0.03 1.0 1.04 0.0 1.0## Hranice oblasti je tvorena 4 useckami bez atributu# 1 hranicni znacku nastavime takto:# leva hranice - znacka 10# prava hranice - znacka 20# dolni hranice - znacka 30# horni hranice - znacka 404 11 1 2 302 2 3 203 3 4 404 1 4 10## V oblasti nejsou diry0

20


Grids/ctverec.poly

Síť vytvoříme pomocí příkazu

triangle -p ctverec.poly

Vzniklou síť lze zobrazit pomocí programu showme, který je součástí balíkutriangle příkazem

showme ctverec.1.ele

Nově vytvořená síť je uložena v souborech ctverec.1.*, kde

• soubor ctverec.1.ele obsahuje v prvním řádku počet trojúhelníků,počet uzlů troujúhelníku (většinou 3, ale pro tzv. kvadratické prvky jeto 6) a počet atributů. Další řádky obsahují číslo trojúhelníku a indexyjeho uzlů v orientaci proti směru hodinových ručiček.

2 3 01 4 1 22 2 3 4

• Soubor ctverec.1.poly obsahuje popis hranice ve stejném formátu,jako měl vstupní soubor. Na rozdíl od něj však neobsahuje souradnicevrcholů trojúhelníků a jejich počet je tedy 0.

0 2 0 14 1

1 2 1 302 3 2 203 4 3 404 1 4 10

0

• Soubor ctverec.1.node obsahuje souřadnice vrcholů s atributy a hra-ničními značkami ve stejném formátu, jako měl vstupní soubor.

4 2 0 11 0 0 302 1 0 303 1 1 204 0 1 40

Vzniklá síť je zobrazena na obrázku 4.5. Vlevo je zobrazena základní síť bezzjemnění. Pro síť uprostřed byla pomocí parametru -a0.002 předepsánamaximální plocha trojúhelníka a nakonec vpravo byla pomocí parametru -qpožadována podmínka minimálního úhlu 20. Tvorbu sítě lze ovlivnit ještěcelou řadou dalších voleb. Jejich popis lze snadno získat spuštěním

triangle -h

21


Obrázek 4.5: Síť v oblasti tvaru čtverce. Vlevo vytvořená pomocítriangle -p ctverec.poly, uprostřed pomocí triangle -a0.002 -pctverec.poly a vpravo pomocí triangle -a0.002 -q -p ctverec.poly.

Síť v GAMM kanále Nyní si ukážeme tvorby sítě v oblasti tvaru kanálus překážkou na dolní stěně (tzv. GAMM kanál). Programem Grids/gamm-tri.py nejprve vytvoříme vstupní soubor Grids/gamm.poly. V něm je po-psáno rozložení bodů na hranicích kanálu odpovídající kroku h = 0.1. Po-kud by síť byla tvořena rovnostrannými trojúhelníky, pak plocha každého znich by byla S = h2√3/4 ≈ 0.0044. Přidáme tedy parametr -a0.0044. Po-žadujme dále, aby všechny úhly byly větší než 25, tedy -q25 a aby programnepřidával žádné další body na hranici -Y. Síť tedy vytvoříme pomocí

triangle -a0.0044 -q25 -Y -p gamm.poly

Výsledná síť je znázorněna na obrázku 4.6.

Obrázek 4.6: Nestrukturovaná síť v GAMM kanále vytvořená pomocí pro-gramu triangle.

Poznamenejme, že program triangle nabízí celou řadu dalších funkcí,které lze využít například pro adaptaci sítě. Na druhou stranu je třeba siuvědomit, že kvalita výsledné sítě je silně ovlivněna algoritmem pro určenípolohy nově vkládaného bodu. Například programem delaunay2 byla prostejně zadané hranice vytvořena síť na obrázku 4.7. Ta obsahuje pouze 357

2Tento program je dostupný pouze na požádání u autora tohoto textu.

22


Grids/gamm-tri.py

Grids/gamm-tri.py

Grids/gamm.poly

uzlů a 632 téměř rovnostranných trojůhelníků oproti 572 vrcholům a 1060trojůhelníkům z programu triangle.

Obrázek 4.7: Nestrukturovaná síť v GAMM kanále vytvořená pomocí pro-gramu delaunay.

4.3 Hybridní sítěNa správně vytvořené strukturované síti je většinou snadné získat velmikvalitní řešení a to i pro případ mezních vrstev, kde jsou používány velmiploché buňky. Na druhou stranu je jejich tvorba pro složité oblasti velmikomplikovaná. Na druhou stranu nestrukturované sítě lze vytvořit bez vět-ších problémů pro téměř libovolný tvar oblasti, avšak použití trojůhelníků vmezních vrstvách je stále předmětem diskusí.

Proto se jako vhodný kompromis jeví tzv. hybridní sítě, které používajístrukturovanou (často hyperbolickou) síť v blízkosti hranice oblasti (tj. tam,kde lze očekávat mezní vrstvy) a nestrukturovanou síť ve zbytku oblasti.Obrázek 4.8 ukazuje detail hybridní sítě v okolí náběžné a odtokové hranyturbínové lopatky.

Obrázek 4.8: Hybridní síť v okolí náběžné a odtokové hrany turbínové lo-patky.

23


4.4 Adaptivní strukturované sítěJednou z možností, jak zlepšit kvalitu numerického řešení, je provést dosta-tečné zjemnění sítě. Pokud bychom však síť zjemňovali rovnoměrně v celéoblasti, vedlo by to k neúměrnému nárůstu počtu buněk sítě. Proto je vý-hodnější síť zjemňovat pouze v oblastech, kde lze očekávat velkou chybu.Určitým problémem však je, tyto oblasti identifikovat. Obrázek 4.9 ukazujenestrukturovanou síť použitou pro výpočet obtékání turbínové mříže SE1050. Vypočtené isočáry Machova čísla jsou na obrázku 4.10. Tento výsle-dek byl analyzován a výpočetní síť byla upravena tak, že buňky ve kterýchdocházelo k velké změně entropie ve smeru proudu byly zjemněny. Na tétosíti byl znovu proveden výpočet a celý cyklus upravy sítě a výpočtu bylaplikován ještě jednou. Obrázky 4.11 a 4.12 ukazují síť po dvou adaptacícha získané izočáry Machova čísla. Z obrázků je vidět mnohem lepší rozlišenístruktury rázových vln na zjemněné síti. Lokálně lokálně zjemněná síť ob-sahovala 11440 buněk. Pokud bychom chtěli získat odpovídající rozlišení vcelé oblasti pomocí globálního zjemnění, potřebovali bychom zhruba dvoj-násobné množství buňek.

4.5 Tvorba sítí ve 3DV předcházející části jsme se věnovali především problematice tvorby sítípro dvourozměrné úlohy. Mnohé z uvedených metod lze rozšířit i na troj-rozměrný případ. Při řešení reálných problémů se ovšem často setkávámes komplikovanou geometrií, jejíž popis je dán výstupem z CAD softwaru.Je-li geometrie složitější, nevystačíme si většinou s jednoblokovou sítí a jetřeba síť tvořit po částech. Výsledný software pro tvorbu sítí je tedy velmikomplikovaný a v důsledku také velmi drahý.

Na tomto místě si ukážeme pouze několik obrázků, které prezentují ně-které kroky procesu tvorby sítě pomocí komerčního software ICEM CFD.Síť byla vytvořena ing. P. Furmánkem následujícím způsobem:

• nejprve byla v programu ICEM CFD vytvořen popis tvaru křídla (vizobr. 4.13),

• dále byla oblast rozdělena na několik bloků (viz obr. 4.14),

• poté byla v každém bloku zvlášť vytvořena síť pomoci TFI (viz obr.4.15),

• nakonec byla celá síť vyhlazena pomocí eliptického generátoru sítě auložena do souboru pro výpočet.

24


Obrázek 4.9: Původní nestrukturovaná síť v okolí lopatky SE 1050.

25


Obrázek 4.10: Izočáry Machova čísla získané na původní nestrukturovanésíťi.

26


Obrázek 4.11: Dvakrát adaptivně zjemněná síť v okolí lopatky SE 1050.

27


Obrázek 4.12: Izočáry Machova čísla získané na dvakrát adaptivně zjemněnénestrukturované síťi.

28


Obrázek 4.13: Popis geometrie křídla v softwaru ICEM CFD.

Obrázek 4.14: Rozdělení oblasti na několik bloků v softwaru ICEM CFD.

29


Obrázek 4.15: Vícebloková síť kolem křídla vytvořená v ICEM CFD.

30


Literatura

[1] Rudolf Dvořák. Transsonické proudění. Academia Praha, 1986.

[2] J. Fořt, J. Fürst, K. Kozel, and P. Louda. Numerické metody řešeníproudových polí ii, 2003.

[3] Karel Kozel and Jiří Fürst. Numerické metody řešení problémů prouděníI. Skriptum ČVUT, srpen 2001. ISBN 80-01-02384-2.

[4] Karel Kozel, Jiří Fürst, and Petr Louda. Numerické metody řešení pro-blémů proudění III. Skriptum ČVUT, srpen 2003.

[5] Joe F. Thompson, Bharat K. Soni, and Nigel P. Weatherill, editors.Handbook of Grid Generation. CRC Press, 1999. ISBN 0-8493-2687-7.

31


Date post:	29-Feb-2020
Category:	Documents
Upload:	others
View:	18 times
Download:	0 times

Počítačová mechanika tekutin - cvut.czmarian.fsik.cvut.cz/~furst/PMT/pmt.pdf · 2014-03-03 ·...

Documents