ČESKÁ ZEMĚDĚLSKÁ UNIVERZITA V PRAZE
FAKULTA ŽIVOTNÍHO PROSTŘEDÍ
DIPLOMOVÁ PRÁCE
Modelování proudění vody na měrném přelivu
Vedoucí práce: Ing. Jiří Pavlásek, Ph.D. Diplomant: Roman Kožín 2009
Prohlášení Prohlašuji, že jsem tuto diplomovou práci vypracoval samostatně pod vedením Ing. Jiřího Pavláska, Ph.D., a že jsem uvedl všechny literární prameny ze kterých jsem čerpal. V Praze 30.7.2009
…………………………
Poděkování Rád bych poděkoval vedoucímu diplomové práce, Ing. Jiřímu Pavláskovi, Ph.D., za rady a připomínky. Dále pak Ing. Martinu Kantorovi a Ing. Aleši Macálkovi, kteří mi pomáhali s praktickou částí řešení ve FLUENTu. Ing. Petru Baštovi za pomoc při interpolaci gridu a Šárce Brdlíkové za pomoc při měření s totální stanicí.
Abstrakt Diplomová práce se zabývá 3D simulací neustáleného proudění v korytě Ptačího potoka na Šumavě, kde byl vybudován měrný profil s Thomsonovým přelivem. Měrný profil i koryto byly zaměřeny za účelem vytvoření výpočetní sítě sloužící k simulaci proudění v potoce a na měrném přelivu. Prostřednictvím variantních výpočtů byly určeny body kozumpční křivky přelivu a její zpřesnění při extrémním průtoku. K tomu bylo využito matematického modelování CFD a softwarového prostředí FLUENT. Přínos práce je v ověření přesnosti měření průtoků a zpřesnění konzumpční křivky v oblasti extrémních průtoků. Klíčová slova: CFD modelování, otevřená koryta, konzumpční křivka
Abstract This thesis deals with 3D simulation of unsteady flow in open channel of Ptaci Potok in Sumava. A V-notch scharp-crested weir was built up there with a hydrometric profile. The weir and the channel were located in order to set up a numerical grid for the flow simulation. Computational fluid dynamics (CFD) and FLUENT software was used for the simulation. Using the variant calculus the points of discharge curve and its extrapolation for extreme discharge can be determined. This type of simulation can verify the actual values of discharge measurements and improves the accuracy of the discharge curve for extreme discharges. Key words: CFD modeling, open channel, discharge curve
Obsah 1. Úvod .......................................................................................................... 8 2. Cíle práce................................................................................................... 8 3. Charakteristika území ................................................................................ 9 3.1 Charakteristiky povodí ........................................................................ 10 3.2 Měrný přeliv........................................................................................ 10 4. Mechanika tekutin .................................................................................... 11 4.1 Newtonské tekutiny ............................................................................ 11 4.2 Nenewtonské tekutiny ........................................................................ 12 4.3 Fyzikální vlastnosti tekutin.................................................................. 12 4.4 Hydrostatika ....................................................................................... 14 4.5 Hydrodynamika .................................................................................. 15 4.5.1 Základní pojmy z hydrodynamiky................................................. 15 4.5.2 Režim proudění ........................................................................... 18 4.5.2.1 Ustálené rovnoměrné proudění ................................................ 18 4.5.2.2 Ustálené nerovnoměrné proudění ............................................ 19 4.5.2.3 Neustálené proudění................................................................. 20 4.5.2.4 Proudění bystřinné, kritické a říční............................................ 21 4.5.3 Rovnice proudění tekutin ............................................................. 24 4.5.3.1 Rovnice kontinuity..................................................................... 24 4.5.3.2 Eulerova rovnice ideální tekutiny .............................................. 26 4.5.3.3 Navier-Stokesova rovnice pro nestlačitelnou tekutinu............... 28 4.5.3.4 Bernoulliho rovnice ................................................................... 29 4.5.4 Proudění v korytech a jeho řešení ............................................... 31 4.5.4.1 Ustálený rovnoměrný průtok ..................................................... 32 4.5.4.2 Ustálený nerovnoměrný průtok ................................................. 33 4.5.4.3 Neustálený průtok ..................................................................... 36 5. Přepad ..................................................................................................... 37 5.1 Přepad přes ostrou hranu................................................................... 37 5.1.1 Obecný tvar rovnice přepadu....................................................... 39 6. Thomsonův přeliv..................................................................................... 40 6.1 Omezení použitelnosti ........................................................................ 43
6.2 Konzumpční křivka při extrémních průtocích...................................... 44 7 Matematické CFD modelování.................................................................. 45 7.1 Výpočetní síť ...................................................................................... 45 7.2 Gambit................................................................................................ 46 7.2.1 Kvalita sítě ................................................................................... 47 7.3 Fluent ................................................................................................. 48 7.3.1 Metoda konečných objemů .......................................................... 49 7.3.2 Numerické řešení turbulence ....................................................... 52 7.3.3 Vícefázové proudění .................................................................... 56 7.3.4 Solver........................................................................................... 57 7.3.5 Konvergence................................................................................ 59 8 Interpolační metoda IDW .......................................................................... 60 9 Metodika ................................................................................................... 62 9.1 Sběr a zpracování dat ........................................................................ 62 9.2 Tvorba geometrie a sítě...................................................................... 63 9.3 Vlastní výpočet ................................................................................... 69 9.3.1 Nastavení výpočtu ....................................................................... 69 10. Výsledky a diskuze ................................................................................ 74 11 Závěr....................................................................................................... 78 12 Použitá literatura ..................................................................................... 79 13 Příloha 1.................................................................................................. 80
Použité značky p – tlak [Pa] T – teplota [°C, °K] Cp – měrné teplo [J/Kg/K]
ρ − hustota [Kg/m3] m – hmotnost [Kg] V – objem [m3]
β – objemová roztažnost χ − objemová stlačitelnost [m2.N-1] η − dynamická viskozita [N.s.m-2] ν – kinematická viskozita [m2.s-1] σ – povrchové napětí [N.m-1] τ – tangenciální napětí [N.m-2] Re – Reynoldsovo číslo [-] R – hydraulický poloměr [-] Fr – Froudovo číslo [-] S – obsah, plocha [m2] v – rychlost [m.s-1] Q – průtok [m3.s-1] Qm – hmotnostní průtok [kg.s-1] F – síla [N] a – zrychlení [m.s-2] l – délka [m] g – gravitační zrychlení [m.s-2]
ζ − ztrátový součinitel [-] ez –ztrátová energie [-] h, z, H – výška, hloubka [m] Ce – přepadový součinitel [-] Θ – úhel v koruně přelivu [°]
α – rychlostní koef. [-] E – specifická energie [-] i0, ie, if – sklon dna, čáry energie, čáry ztrát [-] m – sklon břehů [-] b, B – šířka [m] Kh – koef. vlastnosti kapalin [-] C – Chezyho rychlostní součinitel [m0.5.s-1] n – Manningův součinitel drsnosti [-]
1. Úvod Řešení odtoku z malých lesnických a zemědělských povodí a jeho měření je problém, kterým se zabývá současná hydrologie. Právě zpřesňování
měření
průtoků,
zvláště
pak
měření
popř.
simulace
povodňových vln jsou faktory, které nám pomohou poskytnout informace o chování daného povodí. Diplomová práce se zabývá simulací proudění v korytě Ptačího potoka, přesněji jeho části v pramenné oblasti, kde byl vybudován měrný profil s Thomsonovým přelivem. Ptačí potok odvodňuje experimentální povodí Modrava 2, které se nachází na severním svahu Malé Mokrůvky na Šumavě. Modravská povodí byla vybudována Katedrou vodního hospodářství a Katedrou biotechnických úprav krajiny FLE ČZU v roce 1998 v rámci výzkumných aktivit grantového projektu VaV 620/6/97 „Obnova biodiverzity a stability lesních ekosystémů v pásmu přirozeného rozšíření smrku na území NP Šumava“. V současné době jsou spravována Katedrou vodního hospodářství a environmentálního modelování. (www.kvhem.cz)
2. Cíle práce Práce si klade za cíl zaměřit měrný profil a koryto Ptačího potoka. Vytvořit výpočetní síť a nasimulovat proudění v potoce a na měrném přelivu. Prostřednictvím variantních výpočtů určit měrnou křivku (kozumpční křivku) přelivu a její zpřesnění při extrémních průtocích. K tomu bude využito matematického modelování s numerickými metodami CFD (Computational Fluid Dynamics) a softwarového prostředí FLUENT, který s těmito metodami pracuje. Následně vyhodnocení výsledků simulace a jejich porovnání s výsledky měření a ověření tak správné funkce měrného přelivu.
8
3. Charakteristika území Experimentální povodí Modrava 2 se nachází na severním svahu Malé Mokrůvky v pramenné oblasti Ptačího potoka (hydrologické pořadí povodí 108-01-002), 5 km jižně od Filipovy Huti, na hranici s Bavorskem. Po kůrovcové kalamitě byla v této lokalitě povolena těžba napadeného smrkového porostu. Původní smrkový porost byl starý přibližně 160 let a na části plochy se vyskytoval porost starý 26 let. Porost rovnoměrně pokrýval celou plochu povodí. Po těžbě byla paseka zalesněna smrkem a částečně jeřábem a klenem. V současné době tvoří povrch terénu vysazené a náletové dřeviny, travní porost, tlející větve a pařezy, které zde zbyly po těžbě. Na povodí se jako půdní typy vyskytují především podzoly nebo kryptopodzoly s velkým zastoupením skeletu ve všech půdních horizontech. Hloubka půdního profilu je 0,6–0,8 m (www.kvhem.cz). Z hydrologického hlediska jde o oblast srážkově nadprůměrnou. Spadne zde v průměru 1224 mm za rok. V oblasti je půda obvykle velmi saturována, takže po vydatnějších deštích dochází k povrchovému odtoku. Ten má pak za následek extrémní průtoky. Zatím největší průtok byl změřen 8.8.2008 o hodnotě 277 l.s-1. Průměrný průtok činí 2,68 l.s-1. Co se týče teploty je oblast chladná a její průměrná teplota je 5,52 °C. Minimum je -17,4 a maximum 31,5 °C.
9
Obr. 3.1 Povodí Modrava 2 (www.kvhem.cz, 2009)
3.1 Charakteristiky povodí
Plocha povodí: 0,16 Km2
Min. nadmořská výška: 1197 m.n.m.
Max.nadmořská výška: 1330 m.n.m.
Délka údolnice: 0,745 Km
Sklon svahů: 0,21
(www.kvhem.cz 2009)
3.2 Měrný přeliv Průtok v profilu je vypočítáván z přepadové výšky na trojúhelníkovém Thomsonově přelivu, která je automaticky snímána tlakovým čidlem v časovém kroku 1 hod. Přeliv je s bočními kontrakcemi jak je vidět na (Obr. 3.2).
10
Obr. 3.2 Měrný přeliv (www.kvhem.cz 2009)
4. Mechanika tekutin 4.1 Newtonské tekutiny Jsou tekutiny, které mají lineární závislost mezi tangenciálním napětím a rychlostním gradientem ve směru kolmém k proudu (Kolář et al. 1966). Tzn., že u nich platí Newtonův zákon dynamické viskozity (Obr. 4.2).
τ =η
dv [ N.m − 2 ] dy
(4.1)
11
Obr. 4.1. Rozdělení tekutin (Kolář et al. 1966)
Mezi tyto tekutiny patří většina plynů i kapalin o nízké molekulární tíze, např. voda.
Obr. 4.2 Newtonův zákon dynamické viskozity (Drábková et al. 2007)
4.2 Nenewtonské tekutiny Nemají vztah mezi tečným napětím a gradientem rychlosti lineární. Mezi takovéto tekutiny patří tekutiny dilatantní (silně koncentrované suspenze), pseudoplastické (roztoky polymerů jako jsou polyethylen a polystyren) a Binghamovy plastické hmoty (řídké kaše, bahno, kaly a pasty) (Kolář et al. 1966).
4.3 Fyzikální vlastnosti tekutin
Skupenství vody může být pevné, kapalné a plynné. Bod tání a výparu se obecně mění s tlakem a samozřejmě s teplotou. U vody existuje tzv. trojný bod vody, ve kterém může voda existovat ve všech třech skupenstvích současně (Obr. 4.3). V tomto bodě je rovnováha
12
mezi pevnou, kapalnou a plynnou fází. Souřadnice pro TP jsou p=611,73 Pa a T=273,16 °K (Kolář et al. 1966).
Obr. 4.3 Trojný bod vody (www.wikipedia.org 2009)
Na (Obr. 4.3) je ještě důležitý kritický bod. Ve kterém může existovat voda ve skupenství plynném i kapalném současně. V tomto bodě je tedy rovnováha mezi kapalnou a plynnou fází. Hodnoty pro CP jsou T=374 °C a p=22,064 MPa (navajo.cz 2009).
Měrné teplo vody je teplo Cp, které 1 Kg vody potřebuje k ohřátí o 1°C. Se pohybuje v rozmezí 4217,8 – 4216 J.Kg-1.K-1 pro teploty 0°C – 100°C při tlaku 1013,25 hPa (Kolář et al. 1966).
Hustota, neboli měrná hmotnost je definována jako
ρ=
[
m Kg.m −3 V
]
(4.2)
Nejvyšší hustota vody za normálního tlaku nastává při teplotě 3,98 °C (Kolář et al. 1966).
Objemová roztažnost tekutin je definována jako
β=
1 ∂V V0 ∂T
(4.3)
(Kolář et al. 1966).
13
Objemová stlačitelnost tekutin je definována jako
χ =−
1 ∂V [m 2 .N −1 ] V0 ∂p
(4.4)
(Kolář et al. 1966).
Tepelná vodivost tekutin vyjadřuje schopnost látky vést teplo, nemění-li částice svou polohu vzhledem ke zdroji tepla (Kolář et al. 1966).
Viskozita (vazkost) tekutin vzniká v důsledku tečných napětí a tedy tření, které je způsobeno pohybem sousedních vrstev tekutiny s různými rychlostmi. V kapalině je způsobená kohezí částic a v plynech výměnou hybnosti mezi vrstvami s různou rychlostí. Dynamická viskozita je tedy definována jako
η =τ
dy [ N.s.m − 2 ] dv
(4.5)
Kinematickou viskozitu vyjadřuje poměr dynamické viskozity a hustoty. Tedy:
ν=
η 2 −1 [m .s ] ρ
(4.6)
Se stoupající teplotou viskozita kapalin klesá, kdežto u plynů roste (Kolář et al. 1966).
Povrchové
napětí
volného
povrchu
kapaliny
je
způsobeno
molekulárními silami, které se jej snaží zmenšit. Napětí je dáno vztahem
σ=
dF dl
[N.m −1 ]
(4.7)
vyjadřující účinek kohezních sil mezi molekulami kapaliny vztažený na jednotku délky uzavřené hranice (Kolář et al. 1966).
4.4 Hydrostatika Zabývá se zákony tlaku a jeho rozdělení v kapalinách, které jsou v klidu vzhledem ke stěnám nádoby jež je obsahuje (Kolář et al. 1966). Znamená to tedy, že tvar objemu kapaliny se nemění. Touto částí mechaniky se však práce
nezabývat
nebude.
Práce
je
zaměřena
na
část
druhou
–
hydrodynamiku. 14
4.5 Hydrodynamika Zabývá se prouděním kapalin. Proudění reálných kapalin je složitý problém proto se zavádí zjednodušení ve formě ideální nevazké kapaliny. Proudění se může vyšetřovat v prostoru, rovině nebo po křivce, známé také jako 3D, 2D a 1D proudění (Drábková et al. 2007).
4.5.1 Základní pojmy z hydrodynamiky
Laminární proudění: částice tekutiny se pohybují v tenkých vrstvách, aniž se přemísťují po průřezu viz. (Obr. 4.4). U laminárního proudění v potrubí je rychlostní profil rotační paraboloid viz. (Obr. 4.5).
Obr. 4.4 Laminární proudění (Drábková et al. 2007)
Obr. 4.5 Rychlostní profil (Drábková et al. 2007)
Turbulentní proudění: částice tekutiny mají kromě podélné rychlosti také turbulentní (fluktuační) rychlost, jíž se přemísťují po průřezu viz. (Obr. 4.6). Částice tekutiny neustále přecházejí z jedné vrstvy do druhé, přičemž dochází k výměně kinetické energie a jejich rychlosti po průřezu se značně vyrovnávají. Protože při přemístění částic dochází též ke změně hybnosti, což se projevuje brzdícím účinkem, bude výsledný odpor proti pohybu větší než odpovídá smykovému napětí od vazkosti při laminárním proudění. Rychlostní profil turbulentního proudu v potrubí se proto více podobá obdélníku, a to
15
tím více, čím větší je turbulence (Drábková et al. 2007), tj. čím větší je Reynoldsovo číslo Re viz. (Obr. 4.7).
Obr. 4.6 Turbulentní proudění (Drábková et al. 2007)
Obr. 4.7 Rychlostní profil (Drábková et al. 2007)
Trajektorie je pomyslná čára po které probíhá částice tekutiny. Za ustáleného proudění se trajektorie s časem nemění naopak u neustáleného mohou být trajektorie v každém časovém okamžiku jiné (Drábková et al. 2007).
Proudnice jsou obálkou vektorů rychlostí a jejich tečny udávají směr vektoru rychlosti. U neustáleného proudění vytvářejí proudnice různé částice a nejsou totožné s drahami částic. U ustáleného proudění se nemění rychlosti s časem, a proto mají proudnice stále stejný tvar a jsou totožné s drahami částic (Drábková et al. 2007).
16
Obr. 4.8 Proudnice rychlosti v čerpadle (Soukal et Sedlář 2009)
Proudová trubice je tvořena svazkem proudnic, které procházejí zvolenou uzavřenou křivkou k. Plášť proudové trubice má stejné vlastnosti jako proudnice viz. (Obr. 4.9).
Obr. 4.9 Proudová trubice (Drábková et al. 2007)
Protože směr rychlosti je dán tečnami k proudnicím, je v každém bodě pláště proudové trubice normálová složka rychlosti nulová vn = 0 . Nemůže tedy žádná částice projít stěnou proudové trubice. Proudová trubice rozděluje prostorové proudové pole na dvě části. Částice tekutiny nemohou přetékat z jedné části proudového pole do druhého, a proto platí, že všechny částice protékající průřezem S proudové trubice, musí protékat libovolnými průřezy S1, S2 téže proudové trubice. Jestliže průřez proudové trubice S → 0 , dostane se proudové vlákno. Proudová trubice představuje pomyslné potrubí (Drábková et al. 2007).
17
Reynoldsovo číslo charakterizuje daný proud a režim proudění. Vypočítá se dle vztahu
Re =
v.l
ν
(4.8)
kde v je rychlost, l je charakteristická délka (u otevřených koryt se za l dosazuje hydraulický poloměr
R=
S O
(4.9)
O je omočený obvod), ν je kinematická viskozita. Pro proudění v korytech není hodnota Re rozdělující laminární (Obr. 4.4) a turbulentní (Obr. 4.6) režim. Jedná se spíše o interval <530 - 3450> charakterizující zónu přechodu. Kde může být proudění jak laminární tak turbulentní. Re < 530 zaručuje proudění laminární a Re > 3450 proudění turbulentní (Boor et al. 1968).
Froudovo číslo je dáno výrazem Fr =
v2
(4.10)
g. y
Dle Froudova čísla se rozlišuje proudění bystřinné, kritické a říční, viz. níže.
4.5.2 Režim proudění Kromě proudění laminárního a turbulentního, jejichž definice je uvedena výše, se rozděluje proudění ještě do následujících kategorií.
4.5.2.1 Ustálené rovnoměrné proudění Je neproměnné časově i místně, tedy:
∂v ∂v ∂Q ∂Q =0, =0, = 0, =0 ∂t ∂x ∂t ∂t Může vzniknout jen v pravidelných prizmatických korytech stálého sklonu, jehož všechny příčné řezy jsou stejné a je stálý průtok. Hladina je rovnoběžná se dnem (při zanedbání místních ztrát), takže sklony hladiny a dna se rovnají. Jelikož jsou střední rychlosti ve všech průřezech stejné, bude i čára energie rovnoběžná se dnem (Kunštátský et Patočka 1971). Čili jak
18
uvádí Kolář (1966), jsou-li v rovnováze síly působící pohyb kapaliny a síly tento pohyb brzdící. Pohyb kapaliny způsobuje gradient tlaku nebo složka gravitačního zrychlení působící ve směru proudění.
4.5.2.2 Ustálené nerovnoměrné proudění Tento pohyb můžeme ještě dále rozdělit do dvou kategorií.
Zvolna se měnící Charakterizuje při stálém průtoku zvolna se měnící střední rychlost a tedy v prizmatickém korytě změnu hloubky proudu. Tento pohyb vzniká v prizmatickém kanálu kde je např. překážka proudění jako je jez nebo změna spádu. Vzdutí a snížení vznikající při tomto druhu pohybu závisí především na tření kapaliny o stěny koryta. Pohyb definují tyto rovnice
∂Q = 0, ∂t
∂Q = 0, ∂x
∂v = 0, ∂t
∂v ≠ 0. ∂x
Vlastnosti proudu a koryta se nemění po uvažovanou dobu. Proudnice se
považují
za
rovnoběžné
–
což
umožňuje
předpoklad
hydrostatického rozdělení tlaků v celém průtočném průřezu a dále možnost zanedbat svislou a příčnou složku vektoru rychlosti a vyšetřovat pohyb jako rovinný. Pro střední rychlost platí vztah
v=
Q S
(4.11)
Hloubka po svislici nebo po kolmici ke dnu je prakticky stejná. Součinitel drsnosti nezávisí na hloubce (Kolář et al. 1966).
Náhle se měnící Tento pohyb se vyznačuje především velkou křivostí proudnic, příkladem je třeba vodní skok. V místě náhlé změny křivosti se vytvoří oblast silné turbulence. Díky zakřivení proudnic dochází při proudění k šikmému rozdělení tlaků, které již pak nelze považovat za hydrostatické (Sturm 2001). Rychlá změna je v krátkém úseku, takže tření je zanedbatelné. Náhlá změna pohybu závisí na geometrii překážek. Rozdělení rychlostí v proudu není pravidelné takže neplatí rovnice (4.11). Tvoří se víry a válce takže účinná plocha proudu není
19
dána pevnými stěnami, ale plochou mezi vírovými oblastmi (Kolář et al. 1966).
4.5.2.3 Neustálené proudění Nastává jestliže průtok, rychlost, průtočná plocha a hloubka proudu jsou proměnné, závislé na poloze a čase. Tedy:
∂Q ≠ 0, ∂t
∂Q ≠ 0, ∂x
∂v ≠ 0, ∂t
∂v ≠0 ∂x
Základní typy tohoto proudění jsou oscilační pohyb a translační pohyb.
Oscilační pohyb Je charakterizován kmitáním částic vody kolem rovnovážné polohy bez přenosu průtoku od místa vzniku vlnového pohybu (vlny na hladině, vlny vyvolané nárazem apod.).
Translační (vlnový) pohyb Je kromě vychýlení hladiny z původní polohy charakterizován přenosem průtoku od místa vzniku translačního pohybu (povodňové vlny). Jinými slovy translační vlna je neustálený pohyb v podélném směru, který vyvolává změny průtoku, rychlosti a hloubky v čase (Kolář et al. 1966). Rozlišují se čtyři charakteristiky translačního pohybu. Čelo vlny: přechod mezi původním prouděním a prouděním vyvolaným změnou polohy hladiny či průtoku. V případě pomalu proměnného proudění je čelo vlny prakticky nepostřehnutelné, kdežto u rázových vln zasahuje konečnou oblast. V případě vln poměrně vysokých (ymax : y0 > 1,8 m) vzhledem k hloubce původního proudění je čelo tvořeno překlápějícím se provzdušněným vodním válcem. Pro 1,8 < ymax : y0 < 1,4 je čelo tvořeno menším provzdušněným válcem a zvlněnou hladinou v kratším úseku. Pro ymax : y0 < 1,4 je čelo tvořeno řadou vln, obdobným oscilačním vlnám, jejichž amplituda se postupně směrem od původního proudění zmenšuje, přičemž hladina kolísá kolem určité střední polohy, v kterou postupně přechází (Kolář et al. 1966). Tělo vlny: je oblast za čelem vlny, kde dochází mezi krajním profilem čela vlny a obecným profilem těla vlny ke změně průtoku ΔQ. 20
Ten je definován jako vlnový průtok, tj. průtok přenášený vlnou v této oblasti (Kolář et al. 1966). Absolutní postupivost čela vlny je rychlost, kterou postupuje čelo vlny po proudu nebo proti proudu původniho proudění vzhledem k pozorovateli stojícímu na břehu (Kolář et al. 1966). Relativní postupivost čela vlny je rychlost, kterou postupuje čelo vlny vzhledem k pozorovateli pohybujícímu se rychlostí původního proudění. Pohybuje-li se translační vlna po hladině v klidu, potom je absolutní postupivost rovna relativní postupivosti. (Kolář et al. 1966).
4.5.2.4 Proudění bystřinné, kritické a říční Jak již bylo uvedeno výše, bystřinné, kritické a říční proudění rozděluje Froudovo číslo. Toto číslo je velmi spjato s energií proudění a tak jeho odvození začne definicí specifické energie čili energetické výšky. Specifická energie (energetická výška) Tento termín poprvé zavedl v roce 1912 Bakhmeteff a definoval ho jako
P +z=h ρg
(4.12)
kde h je hloubka. Takže výška specifické energie se definuje jako v2 H0 = h +α 2g
(4.13)
kde α je koeficient rozdělení rychlosti. Je vidět, že specifická energie se rovná součtu hloubky v korytě a rychlostní výšky. Ovšem za předpokladu, že proudnice jsou nezakřivené a rovnoběžné. Platí-li rovnice (4.11) potom
H0 = h +α
Q2 2gS 2
(4.14)
kde S je plocha průřezu a Q je průtok tímto průřezem. Z této rovnice je vidět, že za konstantního průtoku v dané části koryta je specifická energie funkcí pouze hloubky vody (Bos et al. 1976). Grafické vynesení závislosti hloubky vody h a specifické energie dává křivku specifické energie viz. (Obr. 4.10).
21
Obr. 4.10 Vztah hloubky na specifické energii (Bos et al. 1976)
Pro daný průtok a příslušnou specifickou energii jsou dvě „alternativní“ hloubky. V bodě C je specifická energie na minimu pro daný průtok a dvě „alternativní“ hloubky se rovnají. Tato hloubka se nazývá kritickou a značí se hc. Vztah mezi touto minimální specifickou energií a kritickou hloubkou udává diferenciální rovnice, kde průtok Q je konstantní.
dH 0 Q 2 dS v 2 dS = 1−α = 1 − α dh gS dh gS 3 dh
(4.15)
dosadí se dS = B.dh , potom je tedy dH 0 v2B = 1−α dh gS
Je-li specifická energie minimální, platí
(4.16)
dH 0 = 0 a může se psát dh
vc2 S c = α g Bc Tato
rovnice
platí
pouze
za
(4.17)
předpokladu
ustáleného
proudění
s rovnoběžnými proudnicemi a koryta s malým sklonem dna. Je-li rychlostní koeficient α roven jedné, kriterium pro kritické proudění je následující:
vc2 S c = g Bc
(4.18)
čili
22
vc = g
Sc = ghc Bc
(4.19)
kde vc je kritická rychlost, Sc –průřez, Bc – šířka. Výraz pro vc tedy udává Froudovo číslo, které je v tomto případě rovno jedné.
Fr =
vc
(4.20)
ghc
(Bos et al. 1976) Je-li hloubka větší než hloubka kritická proudění se nazývá podkritické (říční) a Fr < 1, je-li nižší než kritická hloubka, proudění je nadkritické (bystřinné) a Fr > 1. Bystřinné proudění tedy charakterizuje malá hloubka a velká rychlost a říční proudění naopak větší hloubka a menší rychlost (Bos et al. 1976). Při říčním proudění je rychlost vody menší než kritická, tedy menší než rychlost šíření vln, které proto mohou postupovat po hladině směrem po proudu i proti němu. Naopak při bystřinném proudění nemůže vlna postupovat proti proudu (Kunštátský et Patocka 1971). Jestliže nastane rychlá změna v hloubce proudu z vyšší hladiny na hladinu nižší, nastává tzv. hydraulický propad. Na druhou stranu, stoupne-li rychle hladina z nižší úrovně na vyšší nastává tzv. hydraulický (vodní) skok, který se projevuje turbulencemi (Bos et al. 1976). A je vždy provázen značnou ztrátou energie. Zpravidla nastává při přechodu z bystřinného do říčního proudění. (Kunštátský et Patočka 1971). Vodní skok je ilustrován na (Obr. 4.11).
23
Obr. 4.11 Vodní skok (Sturm 2001)
Turbulentní víry disipují energii hlavního proudu, mimo to se disipuje také kinetická energie turbulence. Proto je kinetická energie velmi malá na konci vodního skoku. Pro výpočet vodního skoku je vhodné používat rovnice hybnosti, protože přesný matematický popis tohoto proudění je prakticky nemožný. Výpočet a detailnější popis uvádí (Sturm 2001).
4.5.3 Rovnice proudění tekutin 4.5.3.1 Rovnice kontinuity Rovnice kontinuity, často nazývaná také rovnice spojitosti, vyjadřuje obecný fyzikální zákon o zachování hmotnosti. Pro elementární objem, kterým proudí tekutina, musí být hmotnost tekutiny konstantní m = konst. , a tedy celková změna hmotnosti nulová dm = 0 . Celkovou změnu hmotnosti lze dělit na lokální a konvektivní, kde lokální (časová) změna probíhá v elementárním objemu samém (tekutina se stlačuje nebo rozpíná) a konvektivní změna je způsobena rozdílem hmotnosti přitékající a vytékající tekutiny z elementárního objemu. Součet konvektivní a časové změny průtoku je roven nule. Rovnici kontinuity je možné definovat také tak, že rozdíl vstupující hmotnosti do kontrolního objemu a vystupující hmotnosti z kontrolního objemu je roven hmotnosti, která se v tomto kontrolním objemu
24
akumuluje
(Drábková
et
al.
2007)
Tímto
kontrolním
objemem
dV = dx.dy.dz tedy protéká tekutina o rychlosti v = (v x , v y , v z )
Změny způsobené konvekcí - hmotnostní průtok elementem plochy dS je dán vztahem →→
dQm = ρ v n .dS
(4.21)
kde vektor rychlosti se skalárně násobí normálovým vektorem vzhledem k ploše dS , protože průtok je definován v kolmém směru na průtočnou plochu dS . Celkový průtok plochou S je tedy určen plošným integrálem →→
Qm = ∫ ∫ ρ v n.dS
(4.22)
S
Ten se pomocí Gaussova – Ostrogradského věty o divergenci vektoru převede na objemový viz. rovnice (4.23). →→ →→ ⎛ ∂(ρvx ) ∂(ρvy ) ∂(ρvz ) ⎞ ⎟⎟.dx d y dz Qm = ∫∫ρ v n.dS = ∫∫∫div(ρ v n)dx d y dz = ∫∫∫⎜⎜ + + x y z ∂ ∂ ∂ ⎠ V V ⎝ S
(4.23)
(Drábková et al. 2007).
Změny časové - hmotnost je také definována vztahem m = ρV . Jelikož hustota nemusí být v celém objemu konstantní, definuje se na elementárním objemu dm = ρ ⋅ dV , potom je tedy celková hmotnost objemu
m = ∫∫∫ ρ .dV
(4.24)
V
Průtok za čas t je dán vztahem. Qm = ∫∫∫ V
∂ρ ∂ρ dV = ∫∫∫ d x d y d z ∂t ∂t V
(4.25)
Podle zákona zachování hmotnosti (hmotnostního průtoku) platí, že součet konvektivní a časové změny je roven nule. ⎛ ∂ ( ρv x ) ∂ ( ρv y ) ∂ ( ρv z ) ⎞ ∂ρ ⎜⎜ ⎟.d x d y d z = 0 d d d + + + x y ∫∫∫ ∫∫∫ ∂t ∂x ∂y ∂z ⎟⎠ V V ⎝
(4.26)
Jelikož tato rovnice platí pro libovolný objem V, může se rovnice kontinuity zapsat v diferenciálním tvaru 25
∂ρ ∂ ( ρv x ) ∂ ( ρv y ) ∂ ( ρv z ) + + + =0 ∂t ∂x ∂y ∂z
(4.27)
Při proudění kapalin se předpokládá vzhledem k minimálním změnám hustoty že
ρ = konst a
∂ρ =0 ∂t
Potom se dá rovnice kontinuity zapsat ve zjednodušeném tvaru ∂v x ∂v y ∂v z + + =0 ∂x ∂y ∂z
(4.28)
nebo ve vektorovém tvaru →
div ( v ) = 0
(4.29)
(Drábková et al. 2007).
4.5.3.2 Eulerova rovnice ideální tekutiny Rovnice
vyjadřuje →
rovnováhu →
sil
hmotnostních
(objemových),
→
tlakových a setrvačných FS = FO + FP (Obr. 4.14).
Obr. 4.12 Rozdělení sil na elementární objem (Drábková et al. 2007)
26
Diferenciál síly hmotnostní je dán vztahem →
→
→
d FO = a .dm = ρ a .dV
(4.30)
a diferenciál síly tlakové udává vztah →
→
d F p = p n .dS
(4.31)
Celková síla je pak dána objemovým integrálem pro sílu hmotnostní →
→
FO = ∫∫∫ ρ a .dV
(4.32)
V
a plošným integrálem pro sílu tlakovou →
→
F p = ∫∫ − p n.dS
(4.33)
S
Diferenciál setrvačné síly je dán zrychlením →
Dv Dt
(4.34)
→
Dv je tzv. substanciální derivace a rozpis pro jednu složku např. x vypadá Dt
následovně:
Dv x ∂v x ∂v x ∂x ∂v x ∂y ∂v x ∂z ∂v x ∂v x ∂v x ∂v = + + + = + vx v y x vz ∂t ∂x ∂t ∂y ∂t ∂z ∂t ∂t Dt ∂x ∂y ∂z
(4.35)
Pro všechny tři složky je tedy zrychlení →
→
→ Dv ∂v → = + ( v .grad ) v Dt dt →
(4.36)
→
kde ( v .grad ) v představuje zrychlení konvektivní. Převede – li se plošný integrál na integrál objemový dle Gaussovy – Ostrogradského věty, →
→
F p = ∫∫ − p n.dS = − ∫∫∫ gradp.dV S
(4.37)
V
může se psát →
→ Dv ρ .dV = ∫∫∫ ρ a .dV − ∫∫∫ gradp.dV ∫∫∫ Dt V V V
(4.38)
Jelikož vztah platí pro libovolný objem tekutiny, bude platit i pro výraz stojící u integrálu.
27
Tedy: →
Dv → 1 = a − gradp Dt ρ
(4.39)
Nebo ještě lépe →
→ → ∂v → 1 + ( v .grad ) v = a − gradp ρ dt
(4.40)
Což je Eulerova rovnice ve vektorovém tvaru pro neustálené proudění ideální tekutiny (Drábková et al. 2007). Podobné odvození uvádí i (Kolář et al. 1966) nebo (Boor et al. 1968).
4.5.3.3 Navier-Stokesova rovnice pro nestlačitelnou tekutinu N-S rovnice vyjadřuje rovnováhu sil při proudění kapaliny, kdy setrvačná →
→
síla →
→
je
rovna
součtu
síly
hmotnostní,
tlakové
a
třecí
→
FS = FO + FP + Ft . Třecí síla Ft je způsobena viskozitou tekutiny. Základ pro
N-S rovnici tvoří Eulerova rovnice →
→ → ∂v → 1 + ( v .grad ) v = a − gradp dt ρ
(4.41)
→
N-S rovnice se získá přičtením členu νΔ v , který představuje sílu potřebnou k překonání viskózního tření tekutiny. Jeho rozpis pro složku x je následující ⎛ ∂ 2v
∂ 2v
∂ 2v ⎞
νΔv x = ν ⎜⎜ 2x + 2x + 2x ⎟⎟ ∂y ∂z ⎠ ⎝ ∂x
(4.42)
Výsledný vztah N-S rovnice pro nestlačitelnou tekutinu je pak následující →
→ → → ∂v → 1 + ( v .grad ) v = a − gradp + νΔ v dt ρ
(4.43)
Při řešení proudění se zpravidla určuje rozložení rychlostí a tlaků. V systému N-S rovnice a rovnice kontinuity jsou neznámé čtyři veličiny. Složky rychlosti vx, vy, vz a tlak p. Pro řešení těchto rovnic tedy musíme znát vnější zrychlení →
a , hustotu tekutiny ρ a okrajové podmínky. Rovnice se řeší numericky
metodou konečných objemů nebo konečných prvků (Drábková et al. 2007).
28
4.5.3.4 Bernoulliho rovnice Bernoulliho rovnice se dá odvodit z rovnice Navier – Stokesovy, která vyjadřuje rovnováhu sil při proudění reálných tekutin. →
→ → → ∂v → 1 + ( v .grad ) v = a − gradp + νΔ v dt ρ
(4.44)
Předpokladem je jednorozměrné proudění v trubici, tedy předchozí rovnice se zjednoduší tak, že se uvažuje pouze jeden souřadný směr. Vektor rychlosti má jen jednu souřadnici a rozměr je označen l. ∂v 1 ∂v 2 1 ∂p ∂ 2v dl + dl = a. cos ϕ .dl − dl + ν 2 dl ρ ∂l ∂t 2 ∂l ∂l
(4.45)
Jednotlivé členy rovnice odpovídají jednotlivým energiím. Zleva to je energie zrychlení (při neustáleném proudění), kinetická energie, potenciální energie, tlaková energie a ztrátová energie. Integrací výše uvedené rovnice a zavedením potenciálu silového pole dU=a.cosϕ.dl se dostane následující rovnice ∂v 1 ∂v 2 1 ∂p ∂ 2v dl + dl + dl − ν ∫ ∂t ∫ 2 ∂l ∫ ρ ∂l ∫ ∂l 2 dl − ∫ dU .dl = 0
(4.46)
Obr. 4.13 Schéma proudění v trubici
Pro nestlačitelné proudění se vyčíslí integrály pro průřez 1 a 2 proudové trubice.
∫
2
1
⎛ v 2 − v12 ∂v dl + ⎜⎜ 2 ∂t ⎝ 2
2 2 ∂ v ⎞ 1 ⎟⎟ + ( p 2 − p1 ) − ∫ ν 2 dl − (U 2 − U 1 ) = 0 1 ∂l ⎠ ρ
(4.47)
Vyčíslení integrálu vyjadřujícího třecí síly je obtížné, proto se prakticky určuje 29
poloempirickými vztahy a označuje se ez. Představuje práci třecích sil na jednotku hmotnosti proudící tekutiny, což je rozptýlená (disipovaná) energie, nebo též ztrátová energie spotřebovaná na překonání hydraulických odporů na úseku 1 – 2 proudové trubice. Tato ztrátová energie zmenšuje mechanickou energii (tlakovou + kinetickou + polohovou) tekutiny a mění se v teplo (Drábková et al. 2007). Bernoulliho rovnice pro proudění skutečné tekutiny za neustáleného režimu a za předpokladu netlakového proudění, kde působí pouze tíhové zrychlení a = -g a tedy U = -gh je následující 2 ∂v v12 p1 v2 p + + gh1 = 2 + 2 + gh2 + ∫ dl + e z 1 ∂t 2 ρ 2 ρ
(4.48)
Ztrátová energie ez se může vyjádřit jako násobek kinetické energie
v2 ez = ζ 2
(4.49)
nebo tlakové ztrátové energie
ez =
pz
(4.50)
ρ
popřípadě ztrátové výšky ez = hzg
(4.51)
Srovnáním uvedených vztahů se dostane
p z = hz ρg = ζ
v2 ρ 2
(4.52)
ζ je ztrátový součinitel a závisí na druhu hydraulického odporu či ztráty. Bernoulliho rovnice pro proudění skutečné tekutiny za ustáleného režimu, kde
∂v = 0 má tvar ∂t p1
ρ
+
v12 p v2 + gh1 = 2 + 2 + gh2 + e z 2 ρ 2
(4.53)
Bernoulliho rovnice vyjádřena ve výškách, tj. polohové, tlakové, kinetické a ztrátové je na (Obr. 4.14). Rozdíl mezi čarou celkové energie H a čarou mechanické energie představuje rozptýlenou (ztrátovou) energii.
30
Obr. 4.14 Beroulliho rovnice vyjádřená ve výškách (Drábková et al. 2007)
4.5.4 Proudění v korytech a jeho řešení Proudění v říčních korytech je proudění o volné hladině, které vzniká, když omočený obvod netvoří uzavřenou křivku (Kolář et al. 1966). Volná hladina je místem, kde se stýká proud kapaliny s ovzduším (atmosférickým tlakem). Zpravidla se jedná o turbulentní proudění. U otevřených profilů není tečné napětí rozděleno stejnoměrně podél omočeného obvodu a rozdíl tlaků způsobuje v jednotlivých profilech druhotná proudění, která zkreslují rozdělení rychlostí a zvětšují ztráty při proudění. U vodní hladiny jsou tato proudění navíc ovlivněna prouděním vzduchu (což se však často zanedbává) (Drábková et al. 2007). Rychlost proudění se mění jak s hloubkou tak po šířce koryta. Maximální rychlost ovšem není uprostřed koryta na hladině, ale jak je zřejmé z izočar rychlosti (Obr. 4.15), je oblast maximální rychlosti posunuta pod hladinu, což je způsobeno brzděním hladiny o okolní prostředí, tedy o vzduch. Ztráta rychlosti po obvodu je způsobena třením vody o dno a stěny koryta (Drábková et al. 2007). Tření vzniká v důsledku drsnosti, která se zpravidla dost liší v jednotlivých úsecích, kynetě a bermách. Drsnost, která je způsobena jak elementy ve dně (písek, štěrk, kameny), tak vegetací při březích či jinými překážkami jako větve, kmeny apod. je těžko určitelná. Základním a důležitým parametrem po popis proudění v korytě je Froudovo číslo (Sturm 2001).
31
Obr. 4.15 Rychlostní profil v korytě (Drábková et al. 2007)
Komplexnost proudění v korytech často bývá řešena kombinací teorie a experimentů, tak jako v ostatních odvětvích mechaniky tekutin. Musí být splněny základní principy kontinuity, zachování energie a hybnosti. Výsledný vztah pro řešení je často komplikovaný, zejména jsou-li průřezy značně variabilní. Dříve se řešilo neustálené proudění pomocí numerických nomogramů a grafických závislostí jednotlivých proměnných kvůli nelinearitě řídících rovnic a složité geometrii. Dnes nám pro tato řešení pomáhají výkonné počítače spolu s numerickou matematikou (Sturm 2001).
4.5.4.1 Ustálený rovnoměrný průtok Rozdělení hydrostatického tlaku ukazuje (Obr. 4.16). Na hladině (1) je hydrostatický tlak p1 = 0, kdežto v bodě (2) je tlaková výška
p2 = h0 − z tudíž ρg
tlak p 2 = ρg (h0 − z ) . To platí za předpokladu nezakřivených proudnic (Bos et al. 1976).
Obr. 4.16 Rozdělení hydrostatického tlaku (Bos et al. 1976)
32
Obr. 4.17 Rovnoměrný průtok v korytě
Hladina vody je v tomto případě rovnoběžná se dnem koryta viz. (Obr. 4.17). Bernoulliho rovnice pro body na hladině 1 a 2 je stejná jako rovnice (4.53). Pro ztráty třením platí vzorec: hz = z = λ
Δl v 2 d 2g
(4.54)
kde d = 4R a R je hydraulický poloměr. λ.je součinitel tření. Poměrný spád koryta je dán vzorcem i=
λ v2 λ v2 z = = Δl d 2 g 8 g R
(4.55)
Střední rychlost rovnoměrného průtoku v korytě je potom v=
8g
λ
i.R = C i.R
(4.56)
což je Chezyho rovnice a C je tedy rychlostní součinitel. Standardně se ale C vyjadřuje empirickými vztahy na základě drsnosti koryta. Uvedu vztah dle Manninga 1
1 C = R6 n
(4.57)
kde n je stupeň drsnosti (Drábková et al. 2007). C je tedy funkcí Reynoldsova čísla, tvaru profilu, charakteru stěn a Froudova čísla (Kolář et al. 1966).
4.5.4.2 Ustálený nerovnoměrný průtok V místech, kde se spád koryta mění, takže z ≠ hz , vzniká pohyb nerovnoměrný. Při proměnném spádu se průtočná rychlost v a tím i hloubka h mění po délce koryta viz. (Obr. 4.18), nikoliv však v závislosti na čase (Drábková et al. 2007).
33
Obr. 4.18 Nerovnoměrný průtok v korytě (Drábková et al. 2007).
Energie kapaliny v libovolném průřezu je opěr dána Bernoulliho rovnicí, ovšem s tím rozdílem, že musíme uvážit rozdíl výšek hladin h1, h2 . Obecné koryto se rozdělí na úseky délky ΔL, v nichž lze přepokládat, že se průtočné průřezy a tedy i rychlosti spojitě mění z horního úseku k úseku dolnímu. Vychází se opět z rovnice (4.53) a platí tedy vtah Δz +
αv12 2g
=Z+
αv 22 2g
(4.58)
kde Δz je Δz = Z +
α (v d2 − v h2 ) 2g
(4.59)
Δz je rozdíl hladin v obou průřezech, v1 je rychlost v horním průřezu a v2 je rychlost v průřezu dolním. Z je celková ztrátová výška, která se skládá ze ztrát třením Zt a ztrát místních Zm, které vznikají zúžením či rozšířením koryta. Ztráty třením Zt se dají vyjádřit pomocí Chezyho jako
Z t = iΔL =
v s2 Q2 Δ L = ΔL C s2 Rs C s2 Rs S s2
(4.60)
Ztráty místní se vyjádří
Z m = ±ζ
vs 2g
(4.61)
kde ζ je součinitel místní ztráty a znaménka ± ukazují, jedná-li se o zúžení či rozšíření ve směru proudu. Více podrobností v (Kunštátský et Patocka 1971) a (Roub et Pech 2003). (Sturm 2001) uvádí výpočet, kde změny hloubky jsou počítány na základě specifický změn vzdáleností. Opět z Bernoulliho rovnice (4.53) vyplívá 34
H = z+h+
v2 2g
(4.62)
kde z je výška dna od srovnávací roviny a h je hloubka a v je střední rychlost. Změna této výšky ve směru x je
dH dE = −ie = −i0 + dx dx
(4.63)
kde ie je sklon čáry energie, i0 je sklon dna a E je specifická energie.
dE = i0 − i e dx
(4.64)
Jak je vidět z rovnic, specifická energie roste nebo klesá v závislosti na sklonu dna a čáry energie. Platí, že
dE dE dh = dx dh dx
(4.65)
dE 2 = 1 − Fr dh
(4.66)
a z rovnice (4.16) vyplívá, že
následně je možno již určit změnu hladiny podél osy x, tedy
dh i0 − ie = dx 1 − Fr 2
(4.67)
kde Fr je froudovo číslo. Výpočet hloubky h se provede z rovnice
hi +1 − hi =
xi + 1
i0 − ie
∫ 1− F
xi
r
2
dx =
xi + 1
∫ f (h)dx
(4.68)
xi
Z rovnice je vidět, že integrand obsahuje neznámou hloubku h. Je tedy nutné použít alternativy – rozvoj Taylorovy řady pro hi+1, který je následující
hi +1 = hi + f (hi )Δx , kde f (h) =
dh dx
(4.69)
Tato metoda se nazývá Eulerova metoda nebo také metoda prvního řádu, která v obecnosti potřebuje malý časový krok k dosažení dobré přesnosti. Pro dosažení vyšší přesnosti se používá metoda Runge-Kutta, která má přesnost čtvrtého řádu. Více ve (Sturm 2001).
35
4.5.4.3 Neustálený průtok Pro popis neustáleného proudění potřebujeme dvě diferenciální rovnice
reprezentující
principy
kontinuity
a
hybnosti.
Tato
forma
diferenciálních rovnic se nazývá Saint-Venantova rovnice. Většina případů použití těchto rovnic nemá analytické řešení a tak se řeší numericky. Teorie je založená na pohybu tzv. translačních vln. Při říčním proudění se mohou vlny šířit i proti směru proudu. Účelem získání řešení řídících rovnic je popis rychlosti proudění a hloubky jako funkci prostoru a času. Výsledný tvar Saint-Venantovy rovnice je následující ∂Q ∂ ⎛ Q 2 ⎞ ∂ ⎟ + ( ghc S ) = gS (i0 − i f ) + q L v L cos φ + ⎜β ∂t ∂x ⎜⎝ S ⎟⎠ ∂x
(4.70)
kde Q je průtok, S je plocha, hc je kritická hloubka, i0 je sklon dna, if je sklon ztrátové výšky, qL je boční přítok na jednotku délky, vL rychlost bočního přítoku, φ úhel, který svírá vektor rychlosti bočního přítoku s osou x (Sturm 2001). Celé odvození Saint-Venantovy rovnice lze nalézt ve (Sturm 2001) nebo (Havlík et al. 1992). Toto platí pouze pro 1D proudění. Pro 2D proudění, uvádí Havlík et al. (1992) soustavu tří diferenciálních rovnic. 2 2 2 ∂ (v x h ) ∂ (v x h ) ∂ (h + z ) ∂ (hv x v y ) gv x v x + v y 1 + + gh + + − τ wx − fhv y = 0 ρ ∂t dx ∂x ∂y Cx 2 2 2 ∂ (h + z ) ∂ (hv x v y ) gv x v x + v y 1 + + gh + + − τ wy − fhv x = 0 ρ ∂t dy ∂y ∂x Cy 2 ∂h ∂ (v x h) ∂ (v y h) + + =0 ∂t ∂x ∂y
∂ (v y h )
∂ (v y h ) 2
kde Cx je rychlostní součinitel ve směru osy x, Cy ve směru osy y. τwx resp.
τwy je tečné napětí na hladině vlivem větru v ose x resp. v ose y. f je geostropický parametr f=2ωrsinφ, kde ωr je úhlová rychlost otáčení země a
φ je zeměpisná šířka. Pro 3D proudění se obvykle užívá rovnice Navier-Stokesova doplněná o další rovnice resp. software který tyto rovnice řeší numericky. Více v kapitole 7.
36
5. Přepad Přepad je výtok otvorem ve stěně, jehož hrany mohou mít různý geometrický tvar, je-li přepadová výška menší než výška otvoru ve stěně nebo není-li obrys otvoru uzavřen. Přepadová výška h je výška hladiny nad nejnižším místem otvoru – korunou přelivu. Většinou se měří ve vzdálenosti od přelivu 3 – 10 h. Přepad je jev, který vzniká na hydrotechnické konstrukci – přelivu. Stěna přelivu může být tenká ostrohranná nebo masivní široká jak je to vidět u různých hydrotechnických staveb viz. (Kolář.et al. 1966) Základní úlohou je určení přepadové výšky a na jejím základě průtok v daném korytě.
5.1 Přepad přes ostrou hranu Tvar přepadového paprsku přes ostrou hranu závisí na poměru tlaku vzduchu v prostoru pod paprskem k atmosférickému tlaku vzduchu. Při malé přepadové výšce přilne paprsek ke stěně, čímž vznikne paprsek lpící (Obr. 5.1a).
Obr 5.1 Tvary přepadového paprsku (Kolář et al. 1966)
Ten se vytvoří, měníme-li průtok pomalu od malého k většímu. Jeho výkonnost může být až o 30% větší než přepad dobře zavzdušněný jako následek podtlaku v prostoru A. Přibývá-li přepadové výšky, paprsek se od stěny odlepí, a není-li umožněn volný přístup vzduchu, vznikne pod paprskem nestabilní podtlak. Jeli h < 0,4s je podtlakem paprsek stlačován a hladina v prostoru pod paprskem je výše než dolní hladina (Obr. 5.1b). Takový paprsek se nazývá snížený nebo stlačený. Zvětší-li se dále přepadová výška a h > 0,4s , bude všechen vzduch pod paprskem postupně vysáván proudící vodou a dolní voda stoupne, až zaplní celý prostor pod paprskem. Vznikne paprsek spodem ponořený (Obr. 5.1c). Nebrání-li nic 37
volnému
přístupu vzduchu pod paprsek, vytvoří se paprsek volný (Obr.
5.1d). Varianty na (Obr. 12.1a – c), vznikají při nedokonale zavzdušněném přepadu. Tyto jevy jsou při měření průtoků nežádoucí. Na základě dat, která změřil (Howe 1955 in Bos et al. 1976), byli autoři schopni najít vztah, který udává množství vzduchu potřebné pro dokonalé zavzdušnění přepadajícího paprsku.
q air = 0.1
qw ⎛ yp ⎜⎜ ⎝ h1
⎞ ⎟⎟ ⎠
1.5
(5.1)
kde qw je průtok přes přeliv, h1 je přepadová výška a yp je výška vodní hladiny v prostoru pod přepadajícím paprskem, jak ukazuje (Obr. 5.2).
Obr. 5.2 Potřeba zavzdušnění přepadu (Bos et al. 1976)
Za předpokladu dobrého zavzdušnění je přepadový průtok při volném paprsku konstantní a dá se přesně určit. Jedná se o tzv. dokonalý přepad. S tím samozřejmě souvisí i podmínka, že hladina spodní vody nesmí ovlivňovat hladinu vody horní a tím i velikost přepadajícího množství.Tzn. je-li hladina dolní vody níže než koruna přelivu Proto jedině tento paprsek je vhodný k měření průtoků (Bos et al. 1976). V opačném případě vzniká přepad nedokonalý, zatopený. V tomto případě uvedené rovnice neplatí a je třeba zavést opravný součinitel, více o problematice uvádí např. (Kolář et al. 1966).
38
Obr. 5.3 Rozdělení tlaku a rychlosti při dokonalém přepadu (Bos et al. 1976)
Při dokonalém přepadu s plně zavzdušněným prostorem pod přepadovým paprskem nacházíme atmosférický tlak v obou bodech 1 i 2, zatímco v toku před přelivem existuje určité rozdělení hydrostatického tlaku. Tato odchylka 2
je hlavně způsobena hodnotou integrálu
v2 ∫1 gr dn viz. (Obr.5.3).
Klesající piezometrická výška následkem dostředivého zrychlení gr nutně vyvolává vzrůst výšky rychlostní, tedy
v 22 v12 2 v 2 − = dn 2 g 2 g ∫1 gr
(5.2)
(Bos et al. 1976).
5.1.1 Obecný tvar rovnice přepadu Základní rovnice pro přepad bez uvažování přítokové rychlosti se stanoví jako výtok otvorem ve svislé stěně (Obr. 5.4) předepsaný pro integrační meze 0,h tedy jako h
Q = C e 2 g ∫ z y.dz
(5.3)
0
Kde Ce je přepadový součinitel jehož hodnota se mění s výškou h a podle různých vlastností přepadu. Uvážíme-li vliv přítokové rychlosti, dostaneme rovnici h ⎛ v2 ⎞ Q = C e 2 g ∫ ⎜⎜ z + 0 ⎟⎟ y.dz 2g ⎠ 0 ⎝
(5.4)
kde v0 je přítoková rychlost.
39
Obr. 5.4 (Kolář et al. 1966)
Práce se však věnuje přepadu vznikajícímu přes trojúhelníkový přeliv. Tento typ přelivu je zvlášť vhodný pro měření malých průtoků. Při výpočtu průtoků vycházíme ze základní rovnice (5.2), kterou integrujeme pro
y = 2(h − z )tg
Θ 2
(5.5)
je-li trojúhelník rovnoramenný, platí
Θ 8 Q = C e 2 g ⋅ tg ⋅ h 2 15 2 5
(5.6)
kde Θ je úhel v koruně přelivu. Speciální případ trojúhelníkového přelivu je jeho pravoúhlá varianta Θ = 90° zvaná Thomsonův přeliv.
6. Thomsonův přeliv Je jedním z nejpřesnějších měrných přelivů užívaných k měření průtoků v přírodě, zejména na malých vodních tocích. Přeliv a jeho detail zachycují (Obr. 6.1) a (Obr. 6.2).
40
Obr. 6.1 Thomsonův přeliv (Bos et al. 1976)
Obr 6.2 Detail koruny Thomsonova přelivu (Bos et al. 1976)
Přepad má velikost součinitele Ce = 0,5926. Rovnice pro výpočet průtoku je tedy
Q = 1,4h
5 2
(6.1)
Tento vzorec odvodil Thomson již v roce 1861. Zpřesnění výpočtu uvádí (Gourley et Crimp 1915 in Kolář et al. 1966), kde průtok se vypočte dle vztahu
Q = 1,32tg
Θ 2, 47 h 2
(6.2)
A nejpřesnější výsledky jsou v rozmezí výšky přepadového paprsku 0,06m < h < 0,65m, kde však h + s > 3h. 41
Další zpřesnění přinášejí (Kindsvater et Carter in Bos et al. 1976), kteří rovnici (5.6) modifikovali na
Q = Ce
Θ 8 2 g tan he2.5 15 2
(6.3)
kde he je efektivní výška, která se rovná h1 + Kh. Kh reprezentuje vlastnosti kapaliny a je funkcí úhlu Θ. Jak ukazuje (Obr 6.3).
Obr. 6.3 Závislost Kh na úhlu v koruně přelivu (Bos et al. 1976)
(Bos et al. 1976) uvádí, že součinitel přepadu Ce platí pro rozmezí teplot 5 – ⎡h p ⎤ 30 °C a je funkcí tří proměnných C e = f ⎢ 1 , , Θ⎥ . Jestliže jsou poměry ⎣ p B1 ⎦
h1 p ≤ 0.4 a ≤ 0.2 , pak závisí Ce pouze na úhlu Θ jak je vidět na (Obr. 6.4). p B
Obr. 6.4 Závislost Ce na úhlu v koruně přelivu (Bos et al. 1976) 42
Budou-li poměry
h1 p > 0.4 a > 0.2 pak pro odvození Ce pro Thomsonův B p
přeliv můžeme použít grafu na (Obr.6.5).
Obr. 6.5 Koeficienty Ce v závislosti na různých h1/p (Bos et al. 1976)
Koeficienty Ce jsou v případě (Obr. 6.4) s přesností na 1% a v případě (Obr. 6.5) na 1-2%. Tolerance Kh se pohybuje v rozmezí ±0.3mm (Bos et al. 1976).
6.1 Omezení použitelnosti Použitelnost vztahu dle Kindsvater a Carter je možné pouze za předpokladů, které uvádí (Bos et al. 1976):
poměr
h1 by měl být roven nebo menší než 1.2 p
poměr
h1 by měl být roven nebo menší než 0.4 B
výška přepadající hladiny by měla být v rozmezí 0.05 až 0.6 m
šíře koryta by měla být aspoň 60 cm
úhel přelivu mezi 25 a 100 °
spodní voda by měla zůstat pod korunou přelivu
43
6.2 Konzumpční křivka při extrémních průtocích Stanovení konzumpční křivky přelivu pro extrémní průtoku může být v některých případech komplikované. Při překročení návrhového průtoku, dochází ke zkreslení měření a výsledek je tak zatížen určitou chybou. Aby se minimalizovala velikost chyby, musí být použit vhodný koncepční přístup. Ke
stanovení
a
ověření
platnosti
konzumpční
křivky
v oblasti
extrémních průtoků se nabízí několik metod.
matematická extrapolace
hydraulický výpočet
matematické CFD modelování
fyzikální modelování
Obr. 6.7 Konzumpční křivka (Kantor 2007)
Matematická extrapolace spočívá v proložení stávající konzumpční křivky polynomem n-tého stupně a v následné extrapolaci hodnot do oblasti průtoků extrémních. Výhodou této metody je její nenáročnost a je vhodná v případě, kdy proudění není ovlivněno vnějšími vlivy (různé překážky, velká nerovnost terénu apod.) Hydraulický výpočet je založen na výpočtu rovnice přepadu viz. výše, zohledňující jen základní charakteristiky proudu. Tato metoda je použitelná pouze u jednoduchých geometrií.
44
Matematické CFD modelování je založené na metodě konečných prvků resp. konečných objemů a stále více se uplatňuje v nejrůznějších odvětvích. Je možné ji aplikovat na hydrotechnické výpočty, avšak některé úlohy mohou být značně náročné – záleží na charakteru modelu a použitých omezení. Nutností využití této metody je dobré výpočetní zázemí, čas a znalost z oboru CFD. Fyzikální modelování slouží k ověření proudění na přelivu za pomoci modelu ve zmenšeném měřítku. Metoda vychází z mechanické podobnosti dvou fyzikálně stejnorodých dějů. Následně umožňuje extrapolaci výsledků získaných výzkumem do skutečnosti. Metoda je náročná na materiální a prostorové zázemí laboratoří. Více v (Čábelka et Gabriel 1987). (Kantor 2007)
7 Matematické CFD modelování Matematické modelování je jedna z možností jak popsat reálný jev matematickými vztahy. Jednou z oblastí matematického modelování je CFD (Computational Fluid Dynamics), což je modelování proudící tekutiny. CFD je numerické modelování metodou konečných prvků (2D) nebo konečných objemů (3D). Proto je nezbytné před samotným výpočtem vytvořit výpočetní síť pro požadovanou úlohu.
7.1 Výpočetní síť Síť představuje systém rozdělení výpočtové oblasti na dílčí na sebe navazující 2D (v rovině) nebo 3D (v prostoru) elementy. Tyto elementy mohou být trojúhelníky či čtverce ve 2D nebo čtyřstěny či šestistěny ve 3D. Takto vysíťovaná oblast je základem matematického modelování. Neboť samostatný matematický model (systém matematických vztahů) je pouze „pasivním“ nástrojem, který nabývá smyslu až ve chvíli, kdy je aplikován na konkrétní problém - výpočtovou oblast pokrytou sítí (Kozubková 2008). Platí zde několik zásad: výpočet
je
o
to
náročnější,
čím
více
rovnic
je
v rámci
matematického modelu do výpočtu zahrnuto 45
výpočet je o to náročnější, čím více má výpočetní oblast buněk výpočet je o to náročnější, čím méně kvalitní je výpočtová síť Počet buněk patří k hlavním limitujícím faktorům matematického modelování. Jsou případy, kdy se počty výpočetních buněk pohybují v řádu milionů. Proto je cílem s ohledem na budoucí čas tento počet redukovat na nutné minimum. Obrovský nárůst představuje hlavně vytváření mezních vrstev. Redukování počtu elementů by však nemělo být na úkor kvality sítě. Elementy by měly mít rovněž přiměřenou velikost, aby bylo možné jimi v dostatečné míře zachytit modelovaný fyzikální děj (například turbulentní vírové struktury). Proto se používá zhušťování sítě v místech, která jsou z hlediska proudění tekutin pro řešitele zajímavá nebo pro výpočet stěžejní a naopak použití řidší sítě v místech jiných. Zvláštním případem zhuštění buněk je vytvoření tzv. mezní vrstvy v blízkosti stěn, která má za úkol zachytit velké změny fyzikálních veličin u stěny. Zhušťování buněk by mělo být plynulé. Pokud by byla změna ve velikosti buněk provedena příliš velikou skokovou změnou, projevilo by se to znatelně na průběhu výpočtu (problémy s konvergencí úlohy) i konečném výsledku výpočtu (chybný výsledek v daném místě výpočtové oblasti) (Kozubková 2008).
7.2 Gambit Je program na tvorbu geometrií a výpočetních sítí ve 2D i 3D. Umožňuje také importovat geometrie z nejrůznějších CAD/CAE programů. Strukturovaná síť – starší typ, který je složen ze čtyřstěnů nebo šestistěnů ve 3D či z obdélníků nebo čtyřúhelníků ve 2D. Kde základním pravidlem je, že hranice prvků musí sousedit pouze s jedinou hranicí sousedního elementu, nelze tedy libovolně zhušťovat síť – je analogií k metodě konečných diferencí. Výsledná oblast je pak obdélník nebo kvádr. Nestrukturovaná síť – novější typ, kde konečným objemem je kvádr, čtyřstěn, prizmatický a pyramidový prvek. Také umožňuje libovolné zahuštění v bodě, při hraně a v prostoru.
46
Obr. 7.1 Typy elementů (Kantor, 2007)
Jednotlivé typy objemů či elementů je možno kombinovat k dosažení nejlepšího vyplnění prostoru. Pokud je to možné, je výhodné používat obdélníky (2D) a šestistěny (3D), protože jsou při stejné velikosti hrany prostorově úspornější než ostatní elementy. Je také s nimi dosažena vyšší stabilita výpočtu a rychlejší konvergence. Při síťování složitých geometrií se používají trojúhelníkové a čtyřstěnové elementy. Tzv. polyhedron tvoří přechod mezi šesti a čtyřstěny.
7.2.1 Kvalita sítě Kvalita sítě se posuzuje dle velikosti buněk s ohledem na přesnost výpočtu, vhodnosti prostorového uspořádání buněk, kvality buněk z hlediska nesouměrnosti – Skewness a poměru hran prvků – Aspect ratio. Nejvýznamnějším
kriteriem
pro
posouzení
kvality
buňky
je
nesouměrnost, kdy se posuzuje, jak hodně se buňka svým tvarem blíží ideálnímu pravidelnému geometrickému tvaru. Pokud je buňka jakkoliv deformována, je její kvalita horší. Obecně se kvalita každé buňky vyjadřuje bezrozměrným číslem v rozsahu 0 – 1, kde 0 znamená výsledek nejlepší a naopak 1 výsledek nejhorší, tedy problematickou buňku pro výpočty. Tato hodnota se nazývá „míra skosení buňky“ (skewness).
47
Obr. 7.2 Princip posuzování kvality buňky (Kozubková, 2008)
Pro určení kvality 2D buňky, resp. míry její deformace slouží následující vztah
Skewness =
S optimal − S real S optimal
(7.1)
kde Soptimal představuje „optimální plochu buňky“, Sreal „reálnou plochu buňky“ a Skewness je „míra deformace buňky“ vztahující se ke 2D trojúhelníkovému schématu sítě. U jiných schémat se vychází z obdobné logiky. Kozubková (2008) uvádí, že výsledná hodnota by neměla přesáhnout 0.85. Pokud by se tak stalo, je třeba buňku nebo schéma sítě upravit, aby nebyla ohrožena realizovatelnost a přesnost výpočtu. Pro určení kvality 3D buňky čtyřstěnu platí vztah
Skewness =
Voptimal − Vreal Voptimal
(7.2)
kde Voptimal představuje „optimální objem buňky“, Vreal „reálný objem buňky“. Kozubková (2008) uvádí, že výsledná hodnota by neměla přesáhnout 0.9., ale Gambit toleruje deformaci až do hodnoty 0.97. Pokud by však byla hodnota vyšší, je třeba buňku nebo schéma sítě upravit v zájmu dobré konvergence výpočtu. Kvalita sítě se dá testovat v preprocesoru Gambit pomocí příkazu Examine Mesh. Okno poskytuje několik možností, jak kontrolovat lokálně či globálně kvalitu sítě v celé výpočtové oblasti.
7.3 Fluent Je program obsahující matematické modely postihující široké možnosti potřebné k modelování proudění, turbulence, přenosu tepla a reakcí pro průmyslové aplikace. Ty sahají od proudění vzduchu kolem leteckých profilů
48
po spalování v pecích, od modelování probublávání po ropné plošiny, od toku krve po výrobu polovodičů a od návrhu ventilace místností po úpravu a čištění vody. Speciální modely, které dávají softwaru možnosti modelovat multifyzikální úlohy, umožňuje rozšíření působnosti tohoto programu (www.techsoft-eng.cz). Fluent řeší parciální diferenciální rovnice metodou konečných objemů, která spočívá ve třech základních bodech
dělení oblasti na diskrétní objemy užitím obecné křivočaré sítě
bilancování neznámých veličin v individuálních konečných objemech a diskretizace
numerické řešení diskretizovaných rovnic
Fluent definuje diskrétní konečné objemy užitím non-staggered schématu, kdy všechny proměnné jsou uchovávány ve středech konečných objemů.
7.3.1 Metoda konečných objemů Integrace diferenciálních rovnic je vysvětlena na rovnicích o jedné proměnné v ustáleném režimu. Rovnice kontinuity:
∂v =0 ∂x
(7.3)
1 ∂p ∂ ⎡ ∂v ⎤ ∂v 2 =− + ν +S ∂x ρ ∂x ∂x ⎢⎣ ∂x ⎥⎦
(7.4)
Rovnice zachování hybnosti
Rovnice pro přenos skalární veličiny ξ
∂vξ ∂ ⎡ ∂ξ ⎤ = + Sξ αξ ∂x ∂x ⎢⎣ ∂x ⎥⎦
(7.5)
49
Obr. 7.3 Souřadnicové schéma se speciálním značením buněk pro 1D a 3D model (Kozubková, 2008)
Integrací těchto rovnic přes konečné objemy se převedou výchozí diferenciální rovnice na objemový integrál (dV=dx.dy.dz, dA=dy.dz) a užitím Gaussovy – Ostrogradského věty na plošný. Velká písmena označují středy konečných objemů a malá písmena hranice, tj. stěny mezi konečnými objemy, viz (obr. 7.3). Kozubková (2008) uvádí tuto diskretizaci na výsledný algebraický tvar. ∂v
∂v
∫ ∂x dV = ∫ ∂x dxdydz = ∫ vdA = (vA)
V
V
− (vA) w
e
(7.6)
V
Integrace rovnice kontinuity (7.3) vede na tvar
(vA) e − (vA) w = 0
(7.7)
Fyzikálně, výrazy na levé straně označují rozdíl objemových průtoků
Qe − Q w = 0
(7.8)
Integrací rovnice zachování hybnosti (7.4) se získá Qe v e − Q w v w = −
1
ρ
⎛
⎞
⎛
⎠
⎝
⎞
( pe − p w )A + ⎜⎜ ve v E − v P ⎟⎟ A − ⎜⎜ v w v P − vW ⎟⎟ A + SΔV ⎝
Δxe
Δx w
⎠
(7.9) kde A je plocha a S je zdrojový člen.
50
Rovnice pro skalární veličinu (7.5) se upraví shodným postupem na tvar ⎛ ξ − ξP ξ − ξW Qeξ e − Qwξ w = ⎜⎜ α e E −αw P Δx w Δx e ⎝
⎞ ⎟⎟ A + S ξ ΔV ⎠
(7.10)
V předchozích rovnicích se používají veličiny a koeficienty jednak definované ve středech konečných objemů a jednak na stěnách těchto objemů (např. rychlost v rovnici (7.9). To je určitá nevýhoda a je nutné sjednotit ukládání veličin pouze ve středech konečných objemů. Pokud tato veličina bude určena na stěně, použije se interpolační schéma pro interpolaci této veličiny do středu buňky. Pro ilustraci je použito nejjednodušší schéma - aritmetický průměr a diferenční schéma se zjednoduší. Např. rovnici (7.9) lze upravit následovně:
Qe
⎛ v − vE v −v v −v vE + vP 1 − Qw P W = − ( pe − p w )A + ⎜⎜ ve P − vw P W Δxe Δxw ρ 2 2 ⎝
⎞ ⎟⎟ A + SΔV (7.11) ⎠
Rovnice pro skalární veličinu (obecnou proměnou) se upraví takto: Qe
ξP + ξE 2
− Qw
ξ P − ξW 2
⎛ ξ − ξE ξ − ξW = ⎜⎜ α e P −αw P Δxe Δx w ⎝
⎞ ⎟⎟ A + S ξ ΔV ⎠
(7.12)
Pak lze pro tuto obecnou rovnici v jednorozměrném případě vyjádřit ξP pomocí hodnot v sousedních konečných objemech následujícími úpravami
⎛ Qe Qw ⎛ ⎛ A A ⎞ A Qe ⎞ A Qw ⎞ ⎜⎜ − ⎟⎟ξ P = ⎜⎜ − α e + αe −αw − ⎟⎟ξ E + ⎜⎜α w + ⎟⎟ξW + Sξ ΔV (7.13) Δxe Δxw ⎠ Δxe 2 ⎠ 2 ⎝ 2 ⎝ ⎝ Δxw 2 ⎠ Platí že
AP ξ P = AE ξ E + AW ξ W + S C kde AP = − AE − AW − S C Potom rovnici (7.13) lze ještě upravit na
ξ P ∑ (− Ai − S P ) = ∑ Ai ξ i + S C i
(7.14)
i
kde součet se provede přes sousední buňky (v jednorozměrném případě je i= E, W; v trojrozměrném případě i=N, S, E, W, F, B,). Ai jsou koeficienty, které obsahují příspěvky od konvektivních, difúzních a zdrojových členů a SC a SP jsou složky linearizovaných zdrojových členů a Sξ = SC + SP.ξP. Použité označení je patrné z (obr. 7.3). Rovnice řešené ve Fluentu jsou rozšířením předchozích na třídimenzionální křivočarý souřadný systém. Každá iterace sestává z kroků, které jsou zobrazeny diagramem na (obr. 7.4). 51
Obr. 7.4 Algoritmus řešení ve Fluentu (Kozubková, 2008)
pohybové rovnice pro neznámé složky rychlosti jsou řešeny s užitím hodnot tlaků tak, aby se aktualizovalo rychlostní pole
rychlosti určené v předchozím bodě nemohou splňovat rovnici kontinuity, proto se určují tzv. tlakové korekce a následně i korekce rychlostního pole
pomocí nových hodnot rychlostí se řeší rovnice pro turbulentní energii k a disipaci ε
řeší se další rovnice pro určení teploty a dalších skalárních veličin
aktualizují se fyzikální vlastnosti kapalin (např. viskozita)
kontrola konvergence
7.3.2 Numerické řešení turbulence Numerické řešení turbulence se rozděluje na: přímou numerickou simulaci (DNS – direct numerical simulation) simulaci velkých vírů (LES – large eddy simulation) Navier – Stokesova rovnice s Reynoldsovým zprůměrováním (RANS – Reynolds Average Navier – Stokes)
52
Jak ukazuje (obr. 7.5) jednotlivé metody řeší (počítají) vírovou kaskádu až do určité mezní velikosti Δ, víry menší jsou pak již odhadnuty (domodelovány). Samozřejmě, že čím menší víry jsou počítány, tím hustší musí být výpočetní síť. Tím se však stává výpočet náročnějším.
Obr 7.5 Metody řešení turbulence a jejich přesnost (Kantor 2007, Kozubková, 2008)
Jak je vidět z (obr. 7.5) turbulentní proudění je charakterizováno fluktuací rychlostního pole. Jednotlivé výpočetní přístupy mají větší či menší přesnost. FLUENT
používá
spoustu
turbulentních
modelů.
Z hlediska
časové
náročnosti je asi nejlepší použít metodu RANS spolu se standardním k-ε modelem jejichž rovnice a předpoklady jsou uvedeny níže. Reynoldsovy podmínky vi = v i + vi′ vi′ = 0 vi = vi + vi′ = vi
kde v i a vi′ jsou střední a fluktuační rychlost pro i=1,2,3. Pro tlak a další skalární veličiny platí
ξ = ξ +ξ′
53
RANS rovnice je následující
∂p ∂ρvi =0 + ∂xi ∂t ∂ρvi ∂ρvi v j ∂p ∂ + =− + ∂t ∂x j ∂xi ∂x j
⎡ ⎛ ∂vi ∂v j 2 ∂vi ⎞⎤ ∂ ⎟⎥ + + − ρ vi′v ′j − δ ij ⎢ μ ⎜⎜ ⎢⎣ ⎝ ∂x j ∂xi 3 ∂xi ⎟⎠⎥⎦ ∂x j
(
)
(7.15) kde Rij = − ρ vi′v′j je Reynoldsovo napětí. Samotná RANS rovnice je neřešitelná a je třeba nahradit Reynoldsova napětí Boussinesqovou hypotézou. Tato hypotéza předpokládá, že podobně jako při laminárním proudění, kdy platí v zjednodušeném dvourozměrném proudění pro smykové napětí Newtonův vztah, jsou turbulentní napětí a turbulentní toky úměrné gradientu střední rychlosti, teploty, koncentrace apod. Její vztah je následující
⎛ ∂v ∂v j Rij = μ t ⎜ i + ⎜ ∂x ⎝ j ∂xi
⎞ 2⎛ ⎞ ⎟ − ⎜ ρ i + μ t ∂vi ⎟δ ij ⎟ 3⎜ ∂xi ⎟⎠ ⎝ ⎠
(7.16)
kde μt je turbulentní viskozita a δ je Kroneckerovo delta (δij = 0; i≠j). 7.3.2.1 Standardní k-ε model Je nejjednodušší kompletní model turbulence skládající se ze dvou rovnic. Rovnici pro k – turbulentní kinematickou energii a rovnici pro ε – disipaci turbulentní kinematické energie. Tento model je ekonomický, co se týče náročnosti na čas CPU a poskytuje rozumnou přesnost pro široký rozsah turbulentního proudění. Rovnice modelu jsou následující
∂ρk ∂ρkvi ∂ + = ∂t ∂xi ∂x j
⎡⎛ μt ⎢⎜⎜ μ + σk ⎣⎢⎝
⎞ ∂k ⎤ ⎟⎟ ⎥ + Gk + Gb − ρε − YM + S k ⎠ ∂x j ⎦⎥
∂ρε ∂ρεvi ∂ + = ∂t ∂xi ∂x j
⎡⎛ μt ⎢⎜⎜ μ + σε ⎢⎣⎝
⎞ ∂ε ⎤ ε ε2 ⎟⎟ + C ( G + C G ) − C + Se ρ ⎥ k 1e 3e b 2e k k ⎠ ∂x j ⎥⎦
(7.17)
kde Gk reprezentuje vytváření turbulentní kinetické energie následkem gradientu střední rychlosti, Gb reprezentuje vytváření turbulentní kinetické energie následkem vztlaku, Cx jsou konstanty, σx jsou Prandtlova čísla pro k a ε, Sx jsou uživatelem definované členy. Více ve FLUENT User’s guide.
54
7.3.2.2 Řešení v blízkosti pevné stěny V rámci řešení turbulence v blízkosti stěny jsou dvě možnosti. Buď použít
metodu
standardní
stěnové
funkce
nebo
metodu
přímého
numerického řešení. Metoda stěnové funkce je vhodná v případě hrubé sítě v blízkosti stěny. Tato metoda je založena na modelování proudu v oblasti viskózní a přechodové vrstvy. Je doporučeno, aby střed buňky přiléhající ke stěně náležel do intervalu 50≤y+≤500. Metoda přímého numerického řešení řeší přímo proudění ve viskózní a přechodové vrstvě, nutností je ovšem zahuštění sítě v blízkosti stěny tak, aby ve viskózní vrstvě byly alespoň 4 buňky (Kantor, 2007).
Obr. 7.6 Možnosti řešení proudění v blízkosti stěny
Obr. 7.7 Závislost rychlosti u na vzdálenosti od stěny y
55
Hodnota y+ je definována vztahem
y+ =
ρuτ Δy μ
(7.18)
7.3.3 Vícefázové proudění FLUENT umožňuje řešit vícefázové proudění kapalin, plynů a pevných látek a jejich vzájemnou interakci. Pro vícefázové proudění nemísitelných látek slouží model VOF (volume of fluid), který se tedy hodí pro sledování rozhraní dvou látek. Tzn., že je možné jej využít pro prouděné o volné hladině – fáze voda x vzduch. Model VOF počítá proudění pro všechny fáze obdobně, mění se pouze objem fáze ve výpočetní buňce αq
αq = 0 – výpočetní buňka neobsahuje danou fázi αq = 1 – výpočetní buňka je zaplněna danou fází 0 < αq < 1 – výpočetní buňka obsahuje obě fáze (rozhraní fází) Rozhraní fází je řešeno rovnicí kontinuity pro jednotlivé fáze n ⋅ 1 ⎡ ∂α q ρ q ⎛ ⋅ ⎞⎤ + ∇(α q ρ q v q ) = S α q + ∑ ⎜ m pq − mqp ⎟⎥ ⎢ ρ q ⎣ ∂t ⎠⎦ p =1 ⎝ ⋅
(7.19)
⋅
kde mqp je hmotnostní tok z fáze q do fáze p a m pq je hmotnostní tok z fáze p do fáze q. Defaultně je zdrojový člen Sα q roven nule, ale uživatel jej může definovat pro každou fázi zvlášť. Rovnice pro podíl jednotlivých fází není řešena pro primární fázi, která je počítána na základě předpisu: n
∑α q =1
q
=1
(7.20)
Časová diskretizace při užití implicitního schématu je následující:
α qn +1 ρ qn +1 − α qn ρ qn Δt
n ⋅ ⎡ ⎛ ⋅ ⎞⎤ V + ∑ ρ qn +1U nf +1α qn,+f1 = ⎢ S α q + ∑ ⎜ m pq − mqp ⎟⎥V ⎠⎦ f p =1 ⎝ ⎣
(
)
(7.21)
Implicitní schéma dovoluje užití 4 interpolačních metod pro interpolaci hodnoty na stěně buňky do středu buňky a naopak. Jsou to First Order Upwind, Second Order Upwind, QUICK a Modified HRIC.
56
7.3.4 Solver 7.3.4.1 Interpolační metody
First Order Upwind – hodnota ve středu buňky je stejná jako na stěně buňky.
Second Order Upwind (SOU) – interpoluje veličinu φ podle výrazu: r (7.22) φ f , SOU = φ + ∇φ ⋅ r
r kde ∇φ je gradient veličiny φ a r je vektor posunutí ze středu na stěnu buňky.
QUICK – je založen na váženém průměru SOU a středové interpolace proměnné φ. ⎡
⎤ ⎡ S + 2S c ⎤ Sd Sc Sc φP + φ E ⎥ + (1 − θ )⎢ u φP − φW ⎥ Sc + Sd ⎦ Su + Sc ⎣ Sc + Sd ⎣ Su + Sc ⎦
φe = θ ⎢
(7.23)
θ = 1/8
Obr. 7.8 Princip interpolace
Modified HRIC – je přesnější pro VOF výpočty v porovnání s QUICK a SOU. Je však náročnější na výpočetní čas. Více o interpolačních schématech ve Fluent User’s guide. 7.3.4.2 Výpočet gradientu veličin Získání gradientu veličin je nutné nejen pro dopočtení hodnot na stěnách buňky, ale také pro získání rychlostí. FLUENT poskytuje tři metody pro jeho výpočet.
Green – Gauss Cell – Based: stěnová hodnota veličiny φ f je získána
aritmetickým
průměrem
ze
středových hodnot
sousedních buněk.
57
φf =
φ c 0 + φ c1
(7.24)
2
Green – Gauss Node – Based: φ f je počítáno z nodových hodnot dané stěny 1 φf = Nf
Nf
∑φ
n
(7.25)
n
kde Nf je počet nodů na stěně.
Least Squares Cell – Based: změna hodnot mezi středy buněk c0 a ci podél vektoru ri je vyjádřena jako
(∇φ )c 0 ⋅ Δri = (φci − φc 0 )
(7.26)
Obr. 7.9 Princip výpočtu gradientu
7.3.4.3 Výpočetní algoritmy pro tlak – rychlost FLUENT nabízí řadu algoritmů pro výpočet tlakového a rychlostního pole. Základní rozdělení je na segregovaný (Segregated) a sdružený (Coupled) algoritmus. Segregovaný algoritmus obsahuje přístupy SIMPLE, SIPLEC a PISO. Sdruženým algoritmem je Coupled. Jejich detailní popis lze nalézt ve FLUENT User’s Guide. Hlavní rozdíl tkví v tom, že Segregated algoritmus řeší rovnice odděleně. Tato semi – implicitní metoda ústí v pomalou konvergenci. Naproti tomu Coupled algoritmus řeší rovnice dohromady – je tedy rychlejší a lépe konverguje. Dává robustnější a efektivnější řešení v případě ustáleného proudění a je nutný pro neustálené proudění, je-li kvalita sítě malá. Je však náročnější na operační paměť počítače, a to zhruba 1,5x.
58
7.3.5 Konvergence 7.3.5.1 Reziduály Při simulaci proudění pomocí programu FLUENT je velmi důležité získat konvergentní řešení. Mírou konvergence jsou reziduály, které představují maximum rozdílu dvou odpovídajících si veličin ve stejném bodě sítě ve dvou po sobě následujících iteracích. Reziduály jsou vyhodnocovány pro všechny počítané veličiny v každém kroku iterace a zobrazovány pro vybrané veličiny (Kozubková 2008).
Obr. 7.10 Iterace při numerickém výpočtu (Kozubková 2008)
Takže např. pro rovnici
AP ξ P = AE ξ E + AW ξ W + S C
(7.27)
je reziduál definován součtem přes všechna P následovně
R = ∑ AE ξ E + AW ξ W + S C − AP ξ P
(7.28)
P
R je nenormalizovaný reziduál, který má fyzikální rozměr odpovídající rozměru každého členu rovnice a číselně se reziduály, např. pro tlak a rychlost, mohou lišit o řády. Proto se obvykle používá normalizovaný reziduál definovaný následovně:
R=
∑Aξ E
P
E
+ AW ξ W + S C − AP ξ P
∑Aξ P
(7.29) P
P
Reziduály lze vyhodnocovat graficky v každém kroku iterace a snižující se hodnota reziduálu svědčí o dobře konvergující úloze. Obecně řešení velmi
59
dobře konverguje, když se normalizované reziduály snižují řádově k hodnotě 1.10-3 (Kozubková 2008). 7.3.5.2 Relaxační faktory Z důvodu nelinearity diferenciálních rovnic není obecně možné získat hodnoty všech proměnných řešením původně odvozených aproximačních diferenčních schémat. Konvergence lze však dosáhnout užitím relaxace, která redukuje změny každé proměnné v každé iteraci. Jednoduše řečeno, nová hodnota ξP,i+1 v konečném objemu obsahujícím bod P závisí na staré hodnotě z předešlé iterace ξP,i, nové hodnotě z aktuální iterace ξP,i+1,vyp (resp. vypočtené změně ΔξP= ξP,i+1,vyp - ξP,i) a relaxačním parametru α є <0,1> následovně
ξ P ,i +1 = αξ P ,i +1,vyp + (1 − α )ξ P ,i
(7.30)
Tyto relaxační parametry se mohou nastavit pro všechny počítané proměnné. Zvláště pro rychlosti se nastavují velmi malé, řádově desetiny až setiny. Přitom je vhodné během výpočtu tyto hodnoty měnit a tím urychlovat konvergenci, tzn. jestliže změny reziduálů jsou velké při přechodu od jedné iterace k druhé, nastaví se malý relaxační faktor a tím se tlumí nelinearity, pokud se změny reziduálů stávají konstantní, je vhodné relaxační faktory zvětšit (Kozubková 2008).
8 Interpolační metoda IDW Metoda IDW (inverse distance weighted) je založena odhadu hodnoty interpolovaného bodu metodou vážených průměrů z naměřených hodnot okolních bodů, kde váhy jsou reprezentovány inverzními vzdálenostmi od okolních bodů. Základní rovnice pro metodu IDW je ∧
n
z ( p0 ) = ∑ λi z ( pi ) i =1
(8.1)
∧
kde
z ( p0 ) představuje odhad výšky interpolovaného bodu p0 a z( pi )
reprezentuje výšku okolního bodu pi a λi je váha, která přiřazuje hodnotě
60
z( pi ) poměrnou důležitost při procesu interpolace, tzn. jakou měrou ovlivní ∧
výška z( pi ) výsledný odhad výšky z ( p0 ) . Výpočet váhy se řídí rovnicí
λi =
wi
(7.2)
n
∑w i =1
i
Jak bylo uvedeno, hodnocení vlivu okolních bodů na výsledný odhad v bodě interpolovaném, je založen na inverzních vzdálenostech. Tedy čím blíže k interpolovanému bodu se vzorový bod nachází, tím větší váhu získá. Hodnota wi je dána rovnicí
wi = d 0−,iβ
(8.3)
kde d0,i-β je vzdálenost mezi body p0 a pi, β je exponent definovaný uživatelem. Sloučením rovnic (8.1), (8.2) a (8.3) obdržíme tvar n
∧
z ( p0 ) =
∑ z( i =1
pi )
d 0−,βi
n
∑ d 0−,βi
(8.4)
i =1
(Bartier et al., 1996 in Bašta 2008). Zápornost exponentu β způsobuje inverzi délek. Dále platí n
∑λ i =1
i
=1
(8.5)
Výsledný odhad je tedy závislý na počtu n měřených (okolních) bodů, které jsou použity v procesu interpolace, a na hodnotě exponentu β. Počet bodů je obvykle volen z intervalu od jednoho do třiceti, exponent pak z intervalu od jedné do pěti. Konkrétní hodnoty obou parametrů je třeba zvolit tak, aby interpolační technika poskytovala co nejlepší odhady. Volba jejich hodnot se provádí na základě výsledků procesu kalibrace parametrů IDW, který je tedy nutno provést před vlastním procesem interpolace. Proces kalibrace parametrů interpolační techniky se provádí tak, že si pro každý bod ze souboru naměřených dat vypočteme odhady jejich výšek, a to na základě všech kombinací hodnot parametrů n a β. Je-li n = (1; 2; …; 30) a β= (1; 2; 61
…; 5), získáme celkem 150 hodnot odhadů výšek pro každý naměřený bod. Vybrány jsou potom takové kombinace hodnot parametrů n a β, jejichž odhady poskytují celkově nejlepší odhady. Nejlepší kombinaci parametrů lze odvodit buď jednu, která bude platit plošně pro celé zájmové území, nebo (např. při členitější morfologii terénu) na základě prostorové analýzy odchylek mezi odhady a naměřenými hodnotami vytvořit distribuovaný model těchto parametrů měnících se v závislosti na poloze (Bašta 2008).
9 Metodika 9.1 Sběr a zpracování dat Zaměření měrného profilu bylo provedeno totální stanicí Topcon. Koryto a okolí měrného přelivu se měřilo ve čtvercové síti s hustotou přibližně 0.5 x 0.5 m. Změřené body byly přeneseny z totální stanice do počítače programem Geoman. Následně se musely změřené body interpolovat do pravidelné sítě, vhodné pro import do programu Gambit. Před interpolací byly ze souboru dat vyjmuty body zaměřeného přelivu, aby nedošlo ke zkreslení. Samotná interpolace byla provedena pomocí naprogramovaných
skriptů
v programu R, které pocházejí z diplomové práce Petra Bašty (2008). Interpolovalo se metodou IDW kde každý bod byl vypočítán z okolí 7 bodů.
Obr. 9.1 Kritérium pro výběr interpolačních parametrů
62
Pomocí předprogramovaných skriptů byla odvozena jedna kombinace parametru n a β pro celé zájmové území a to n=7 a β=1 viz. (Obr. 9.1). Výsledkem byl vyinterpolovaný grid 0.2 x 0.2 m viz. (Obr. 7.2), který byl pak použit pro vlastní vytvoření výpočetní sítě.
Obr. 9.2 Povrch vytvořeného gridu
9.2 Tvorba geometrie a sítě Výpočetní síť byla vytvořena v software GAMBIT a to tak, že výsledný vyinterpolovaný grid ve formě xyz byl načten do Gambitu přes Open – import – vertex.
Obr. 9.3 Grid zájmového území 20 x 20 cm
63
Následně bylo nutné grid upravit na stejný počet řádků i sloupců aby se dala vytvořit plocha (face) pomocí funkce Create Face From Vertex Rows. Dále byly importovány body přelivu. Tyto body se pospojovaly úsečkami – Create Edge, kde spodní část přelivu byla protažena až pod plochu tvořící dno, tak aby se obě plochy protínaly. Následně byla vytvořena plocha přelivu funkcí Create Face From Edges. Už zbývá jen tyto plochy (oříznout) spojit do geometrie, aby na sebe přesně navazovaly. K tomu slouží funkce Subtract Faces.
Obr. 9.4 Tvorba dna a přelivu
Aby se s geometrií dalo dále pracovat je nutné plochu umístit do souřadnic tak, aby středem koryta procházela osa x a vstup ležel v souřadnicích 0,0,0. K tomu je určená funkce Align Face, více v Gambit User’s Guide. Jelikož toto území je příliš velké na výpočet je třeba ho oříznout tak, aby zůstalo koryto s přelivem – oblast hlavního zájmu. To se provede opět funkcí Subtract Faces podle (Obr. 9.5). Vytvořenou pomocnou plochou ořízneme jak dno tak přeliv. 64
Obr. 9.5 Postup ořezu
Všechny zbylé plochy se smažou funkcí Delete Face. Teď se již může nad touto geometrií vystavět objem. A to tak, že se rozkopírují hraniční body ve směru osy z do určité výšky, aby tvořily rovinu. K tomu slouží funkce Move/Copy Vertices. Následně se všechny body pospojují. Z jednotlivých spojnic se vytvoří v logickém pořadí plochy nástrojem Create Face from Wireframe. Z těchto ploch se zkonstruuje objem pomocí funkce Stich Faces (Obr. 9.6).
Obr. 9.6 Výpočetní objem
65
Za přeliv byl připojen další objem, kvádr který nahradil zbývající část koryta. Jelikož se ukázalo během testování výpočtů, že je výhodné, když úsek koryta blízko spodní okrajové podmínce je nezvlněný a pravidelný. Pokud tomu tak není nastávají značné problémy s konvergencí a funkcí okrajové podmínky. Výpočet buď po určité době spadne nebo se okrajová podmínka začne chovat jako nepropustná. Druhý objem se vytvoří a připojí dle postupů uvedených výše. Dalším krokem je nastavení okrajových podmínek. GAMBIT spolu s FLUENTem umožňují nastavit několik typů okrajových podmínek. Jejich použití záleží na typu úlohy a na informacích, které má uživatel ohledně vstupu, výstupu, teploty apod. V případě této úlohy jsou k dispozici okrajové podmínky pro vstup Velocity Inlet a Mass Flow Inlet. Vybrána byla Mass Flow Inlet (modře), jelikož je znám průtok. Na výstupu se naskýtá použití podmínek Outflow a Pressure Outlet. Byla použita Pressure Outlet (červeně), protože je citlivější a v kombinaci s Mass Flow Inlet dokáže jako jediná definovat sklon dna. Všechny ostatní plochy byly definovány jako Wall (šedě a zeleně) viz. (Obr. 9.7). Okrajové podmínky je samozřejmě možné definovat ve FLUENTu. Více k tématu ve FLUENT User’s Guide.
Obr. 9.7 Kompletní model koryta
66
Předposledním krokem je tvorba sítě. Správné vytvoření sítě je jednou z nejdůležitějších věcí celé úlohy. Jelikož správnost a přesnost výpočtu, ale i doba jeho trvání závisí na hustotě a celkovém uspořádání sítě. Síť byla vytvořena operací Mesh Volumes, které předcházelo vytvoření Size function dle (Obr. 9.8), kde hrany přelivu byly zahuštěny na 2 cm a postupně se hustota směrem od přelivu zmenšovala faktorem 1,1 až do finální velikosti 10 cm.
Obr. 9.8 Postup vytvoření Size function
Obr. 9.9 Vytvořená síť
67
Vzhledem ke komplikovanosti geometrie byly pro síťování zvoleny čtyřstěny, které lépe vyplní prostor (Obr. 9.9). Jejich nevýhodou je, že jich muselo být použito při stejné konfiguraci 4x více než v případě hexa prvků. Celkový počet činil 180 000 buněk. Což v případě takového výpočtu není mnoho. Ale vzhledem k nedostatku času i výpočetní kapacity a hlavně licencí, kdy se výpočet dal paralelizovat pouze na 2 části je to dost. Bohužel kvůli tomu byla snížena ostrost fázového rozraní tekutin (hladiny). Posledním krokem je zhodnocení kvality sítě. Jak již bylo uvedeno, neměla by míra zkosení buňky přesáhnou hodnoty 0.9 a v žádném případě hodnotu 0.97. Ke zkoumání kvality sítě slouží nástroj Examine Mesh. Ten umožňuje zobrazení kvality buněk v barevné škále a prostorové prohlédnutí sítě ve všech třech směrech, viz. (Obr. 9.10).
Obr. 9.10 Kvalita sítě ve směrech x,y, z.
Nejhorší elementy (červeně) mají zkosení 0.74 a jejich počet je 42 což je 0.02% z celkového počtu. Na základě tohoto faktu se dá síť klasifikovat jako 68
vhodná pro výpočet. Na závěr se síť vyexportuje, aby se dala použít pro výpočet ve FLUENTu. Provede se tak přes File – Export – Mesh.
9.3 Vlastní výpočet Výpočet probíhal v programu Fluent. Úloha byla řešena trojrozměrně a nestacionárně. Pro výpočet byl použit model vícefázového proudění (VOF). Testovací výpočty probíhaly na počítači se 4GB operační pamětí s procesorem Intel Core 2 Quad Q9300, což je 4 jádrový procesor s taktovací frekvencí 2,5 GHz, 6MB L2 Cache a frekvencí sběrnice 1333MHz. Následně se naskytla možnost využití stroje se 2GB operační paměti a s procesorem Intel Xeon X5550 2.67GHz, který má 8 jader, 8MB L2 Cache a FSB 1066 MHz a všechy další výpočty probíhaly zde. Tento procesor byl o poznaní rychlejší a jeden výpočet na něm trval necelých 72 hodin při napočítaném čase 600 sekund. Celkem se počítaly 3 průtoky, které byly vybrány z měření v roce 2008. Dva průtoky byly kalibrační, sloužily pro určení výšky hladiny v koruně přelivu. Tzn., jaký poměr fáze voda – vzduch tvoří skutečnou hladinu. Tyto průtoky byly 53.07 l.s-1 z 22.04.2008 v 16:00 a 100.59 l.s-1 z 1.03.2008 ve 14:00. Průtokem, který bylo třeba ověřit, bylo doposud maximálně naměřených 277 l.s-1 ze 7.08.2008 v 0:00. Průtoku jsou měřeny pomocí thomsonova přelivu, který má délku ramena 74.1 cm, maximální měřitelnou výšku h1 od koruny 52.4 cm a výšku koruny ode dna p = 45.57 cm a střední šířku B = 337.3 cm. Celková délka modelovaného úseku je 12.63 m a průměrný sklon dna činí 11%.
9.3.1 Nastavení výpočtu Pře File – Read – Case se nahraje do FLUENTu vytvořená síť *.msh. Poté je vhodné síť zkontrolovat Grid – Check. Dalším krokem je už vlastní nastavení modelů a solveru přes Define – Model viz. (Obr. 9.11, 9.12, 9.13). Schéma pro výpočet gradientu bylo zvoleno Least Square Cell Based popsané v kapitole 7.3.4.2. VOF schéma lze pro Open Channel použít pouze
69
implicitní – to je však výhodné, jelikož umožňuje použití většího časového kroku a je stabilnější na výpočet. .
Obr. 9.11 Nastavení modelů
70
Obr 9.12 Nastavení okrajových podmínek
71
Obr 9.13 Nastavení Solveru
Jako model turbulence byl použit standardní model k – ε se standardní stěnovou funkcí, který je popsán v kap. 7.3.2. Proudící médium – voda bylo vybráno z databáze fluentu. Hodnoty jsou platné pro standardní teplotu 20 °C. Drsnost povrchu se nastavuje dost netradičně přes Wall – Roughness ve vlastnostech stěny. Ve FLUENT User’s Guide doporučují pro koryta použít hodnotu 0.5 a nastavit jakousi požadovanou výšku této drsnosti. Po dohodě s Jirkou Pavláskem bylo nastavena výška 1cm. Nastavení na (Obr. 9.13) odpovídá použitým výpočetních algoritmů, interpolačních schémat a relaxačních faktorů, které jsou vysvětleny v kap. 7.3. Výpočet se začne inicializací Solve – Initialize – Initialize, kde se „nastřelí“ počáteční hodnoty parametrů pro celou oblast. Rychlosti na hodnotu 0, přetlak na 0, turbulentní kinetická energie a její disipace na 1 a podíl vody taky na 1. To nám zajistí vyplnění celého výpočetního prostoru vodou. Následně se provede adaptace regionu Adapt – Region, která umožní definovat prostor, který má být vyplněn vzduchem. Bylo zadáváno od výšky 0.5m po celém výpočetním prostoru. Vybrané buňky se označí Mark a následně se vybraná oblast přepíše hodnotou nula pro voda volume fraction přes Solve – Initialize – Patch, dle (Obr. 9.14). Tzn. že vybraná oblast bude obsahovat pouze vzduch. Grafické zobrazení výpočtu v čase v t = 0 ukazuje (Obr. 9.15). 72
Obr. 9.14 Označení regionu a nastavení požadované hodnoty
Obr. 9.15 Počátek výpočtu v čase t = 0
Paralelizace úlohy se vytvoří přes Parallel – Partition, kde se vybere metoda rozdělení a počet částí, které mají vzniknout. Následně se tento projekt uloží do souboru *.cas. Pro zpuštění paralelního výpočtu je nutné FLUENT spouštět z příkazového řádku s parametry 3D a –t2 (fluent.exe 3D – t2). Po spuštění se načte soubor *.cas a výpočet se spustí přes Solve – Iterate. Výpočet musí začínat na velmi malém časovém kroku, aby bylo dosaženo konvergence. Konverguje-li výpočet může se postupně časový
73
krok zvyšovat. Výhodné je když výpočet konverguje tak po 5 až 10 iteracích. Míra konvergence se sleduje pomocí reziduálů jednotlivých veličin, ty se zobrazí v Solve – Monitor – Residuals.
Obr. 9.16 Nastavení iterací a monitoring reziduálů
Výpočet byl spuštěn s časovým krokem 1.10-6 s a průběžně pomalu zvyšován až na krok 0,05 s. Počet iterací během výpočtu byl udržován v rozmezí 5 – 10 na jeden časový krok.
10. Výsledky a diskuze Jak již bylo předesláno, výpočty jsou pro 3 průtoky Q1 = 53.07 l.s-1, Q2 = 100.59 l.s-1 a Q3 = 277 l.s-1. První dva jsou kalibrační, kdy bylo potřeba zjistit jaký poměr fází voda – vzduch vlastně tvoří hladinu. Jestli to je 20% vody v buňce nebo 15% ? A na základě známé hladiny potom ověřit jestli průtok 277 l.s-1 je změřený správně či nikoli popř. s jakou chybou. K této jednoduché analýze poslouží grafy závislosti poměru fází ku přepadové výšce vody na přelivu. Grafy zachycuje (Obr. 10.1). Tabulka (Tab. 10.1) ukazuje jednotlivé průtoky, kde h1 je přepadová výška, Q je průtok a S jsou srážky. date 22.4.2008 1.3.2008 7.8.2008
time 16:00 14:00 0:00
h1 (mm) Q (l/s) 0.272 53.07654 0.3515 100.5943 0.5275 277.0086
S (mm) 0 73.6
Tab. 10.1 Tabulka průtoků
74
poměr fáze
1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.45
0.55
0.6
0.65
přepadová výška [m]
0.5
0.7
0.75
0.8
0.85
Závislost poměru fází na přepadové výšce
0.4
0.9
0.95
1
1.05
1.1
1.15
Q1 53.07 l.s-1
Q2 100.59 l.s-1
Q3 277 l.s-1
Obr. 10.1 Grafy závislostí poměru fází na přepadové výšce
Z grafu je vidět že pro průtok Q1 s přepadovou výškou 0.272m je poměr fází
0.1435 a pro průtok Q2 s přepadovou výškou 0.3515m činí hodnota 0.1516
75
(přesné hodnoty byly získány lineární interpolací). Pro získání skutečné přepadové výšky průtoku Q3mod zprůměrujeme hodnoty 0.1435 a 0.1516. Výsledná hodnota je 0.1475. Zprůměrováním vznikla chyba ± 1mm na přepadové výšce pro Q3mod, která činí 0.5019 m. Rozdíl oproti původní výšce 0.5275 m měřeného průtoku je tedy 0.0256 m. Souhrn výsledků uvádí (Tab. 10.2.). frakce prep. vyska prutok
Q1 Q2 Q3simul Q3mer 0.1435 0.1516 0.1475 0.2720 0.3515 0.5019 0.5275 53.0746 100.5907 277.0000 277.0000
Tab. 10.2 Souhrn kalibračních a modelovaných průtoků
Výpočet průtoku je realizován dle rovnice (6.3) (s koeficientem Ce = 0.578), která pro přepadovou výšku 0.5275 m udává průtok 277 l.s-1. To je ale v rozporu pro modelovaný průtok Q3mod 277 l.s-1.pro který je přepadová výška 0.5019 m. Když se do vzorce (6.3) dosadí výška 0.5019 vychází průtok 244.27 l.s-1.Z toho vyplívá, že při použití rovnice (6.3) pro výpočet průtoku přes Thompsonův přeliv na povodí Modrava2, dochází k podhodnocení větších a extrémních průtoků. V konkrétním případě průtok 277 l.s-1 podhodnotí o 32.73 l.s-1, což je významný rozdíl. Znamená to tedy, že při změřené přepadové výšce 0.5275 muselo téci přes přeliv určitě větší množství než 277 l.s-1. Při stávajícím trendu konzumpční křivky (Obr. 10.2) lze odhadnout průtok v rozmezí 300 – 310 l.s-1. Tento rozdíl je zřejmě způsoben vlivem přítokové rychlosti, kterou rovnice (6.3) neuvažuje. Nádrž před přelivem je totiž relativně malá, přičemž platí zásada, čím menší nádrž tím větší vliv přítokové rychlosti. Na obrázcích v příloze 1 lze pozorovat zvýšení rychlosti v blízkosti přelivu zejména v případě průtoku 277l.s-1. Dále jsou v příloze obsaženy grafické výstupy sledovaných veličin v korytě a na přelivu – rychlost, tlak, poloha hladiny.
76
h1 [m]
0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0
25
50
75
125
150
175
200
225
Konzumpční křivka měrného přelivu
100
-1
Q [l.s ]
250
275
300
Měřená
Kalibrační
Modelovaná
Obr. 10.2 Konzumpční křivka měrného přelivu na Modravě 2
77
11 Závěr Práce si kladla za cíl využití metody CFD analýzy na praktickém příkladě proudění tekutin. CFD modelování je určitě silný nástroj pro vyšetřování proudění. Ikdyž příprava veškerých podkladů pro uskutečnění výpočtu je celkem náročná. Výpočet vyžaduje hodně trpělivosti a času, zejména má-li řešitel malé zkušenosti z této oblasti. Přesnost výpočtu je silně závislá na počtu výpočetních elementů. Avšak při dobře sestaveném modelu a dostatečném výpočetním výkonu se může dosáhnout výborné přesnosti. I tak ale výpočet zabere řádově dny. Metoda CFD
je
použitelná
v případě,
máme-li
k dispozici
nějaká
kalibrační
(srovnávací) data ať z fyzikálního modelování nebo ze skutečného měření. V opačném případě je CFD modelování těchto typů úloh nevyhovující. Nicméně tato forma simulace může poskytnout náhled do procesů v proudění, které se jinak těžko odhalují Pro tento typ úlohy se tato metoda osvědčila nejen v této práci, ale např. i v (Kantor 2007). Takže bych tuto metodu doporučil i při navrhování způsobů měření průtoků, kde může být užitečným pomocníkem.
78
12 Použitá literatura BOOR B. (ed.), 1968: Hydraulika pro vodohospodářské stavby, Praha. BOS M.G. (ed.), 1976: Discharge Measurement Structures, Wageningen. BAŠTA P., 2008: Digitální model terénu povodí Modrava 2, Nepublikováno. ČÁBELKA J. et GABRIEL P., 1987: Matematické a fyzikální modelování v hydrotechnice, Praha. DRÁBKOVÁ S. (ed.), 2007: Mechanika tekutin, Ostrava. HAVLÍK V. (ed.), 1992: Matematické modelování neustáleného proudění, Praha. KANTOR M., 2007: Hydraulika bezpečnostních přelivů vodních děl za extrémních průtoků, Nepublikováno. KOLÁŘ V. (ed.), 1966: Hydraulika, Praha. KOZUBKOVÁ M., 2008: Modelování proudění tekutin, Ostrava. KUNŠTÁTSKÝ J. et PATOČKA C., 1971: Základy hydrauliky a hydrologie pro inženýrské konstrukce a dopravní stavby, Praha. ROUB R. et Pech P., 2003: Hydraulika příklady, Praha. ŠOUKAL J. et SEDLÁŘ M., 2009: CFD analýza článkových čerpadel v turbínovém režimu, Nepublikováno. STURM T. W., 2001: Open Channel Hydraulics, New York. Fluent Inc., 2006: FLUENT User’s guide. Fluent Inc., 2006: Gambit User’s guide.
Internetové zdroje www.techsoft-eng.cz www.wikipedia.org www.kvhem.cz www.navajo.cz www.cfd-online.com
79
13 Příloha 1 Sledované veličiny pro průtok 53.07 l.s-1 v čase t = 100s
Zobrazení průběhu hladiny v podélném profilu a na přelivu. Stupnice ukazuje podíl fází.
Rychlost [m.s-1] a tlak [Pa] v podélném profilu
Rozdělení vektorů rychlosti [m.s-1] v blízkosti přelivu u a hladina s rychlostmi [m.s-1]
80
Histogram rychlostí a tlaků v celé oblasti
V čase t = 300s
Zobrazení průběhu hladiny v podélném profilu a na přelivu. Stupnice ukazuje podíl fází.
Obr. Rychlost [m.s-1] a tlak [Pa] v podélném profilu
81
Rozdělení vektorů rychlosti [m.s-1] v blízkosti přelivu u a hladina s rychlostmi [m.s-1]
Histogram (rozdělení) rychlostí a tlaků v celé oblasti
V čase t = 600s
Zobrazení průběhu hladiny v podélném profilu a na přelivu. Stupnice ukazuje podíl fází.
82
Rychlost [m.s-1] a tlak [Pa] v podélném profilu
Rozdělení vektorů rychlosti [m.s-1] v blízkosti přelivu u a hladina s rychlostmi [m.s-1]
Histogram (rozdělení) rychlostí a tlaků v celé oblasti
83
Sledované veličiny pro průtok 277 l.s-1 v čase t = 100s
Zobrazení průběhu hladiny v podélném profilu a na přelivu. Stupnice ukazuje podíl fází.
Rychlost [m.s-1] a tlak [Pa] v podélném profilu
Rozdělení vektorů rychlosti [m.s-1] v blízkosti přelivu u a hladina s rychlostmi [m.s-1]
84
Histogram (rozdělení) rychlostí a tlaků v celé oblasti
V čase t = 300s
Zobrazení průběhu hladiny v podélném profilu a na přelivu. Stupnice ukazuje podíl fází
Rychlost [m.s-1] a tlak [Pa] v podélném profilu
85
Rozdělení vektorů rychlosti [m.s-1] v blízkosti přelivu u a hladina s rychlostmi [m.s-1]
Histogram (rozdělení) rychlostí a tlaků v celé oblasti
V čase t = 600s
Zobrazení průběhu hladiny v podélném profilu a na přelivu. Stupnice ukazuje podíl fází
86
Rychlost [m.s-1] a tlak [Pa] v podélném profilu
Rozdělení vektorů rychlosti [m.s-1] v blízkosti přelivu u a hladina s rychlostmi [m.s-1]
Histogram (rozdělení) rychlostí a tlaků v celé oblasti
87
