ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE – FAKULTA STAVEBNÍ Katedra hydrauliky a hydrologie
Diplomová práce Hydraulika bezpečnostních přelivů vodních děl za extremních průtoků
MARTIN KANTOR
Školitel: Ing. Petr SKLENÁŘ, Ph.D.
Praha, leden 2007
Martin Kantor
Diplomová práce
Prohlášení : Prohlašuji, že tuto práci jsem vypracoval samostatně pouze za odborného vedení vedoucího diplomové práce Ing. Petra Sklenáře, Ph.D. a konzultantů: Ing. Zdeňka Cháry, CSc., Doc. Ing. Jiřího Béma, CSc. Dále prohlašuji, že veškeré podklady, ze kterých jsem čerpal, jsou uvedeny v seznamu použité literatury.
V Praze 8.1.2007
Martin Kantor
Martin Kantor
Diplomová práce
Poděkování Rád bych touto cestou vyjádřil svůj dík Ing. Petru Sklenářovi, Ph.D. za jeho cenné připomínky, trpělivost a ochotu při vedení mé diplomové práce. Rovněž bych chtěl poděkovat Ing. Zdeňku Chárovi, CSc. a Doc. Ing. Jiřímu Bémovi, CSc., kteří mi vyšli maximálně vstříc a umožnili mi přístup ke všem potřebným informacím.
Tyto výsledky jsou součástí grantového úkolu GAČR 103/04/1328 "Nejistoty hydraulických výpočtů na vodních tocích pro extrémní hydraulické jevy", řešeného na katedře hydrauliky a hydrologie, Fakultě stavební, ČVUT v Praze, a aktivit výzkumného centra CIDEAS v rámci projektu 1M6840770001 MŠMT ČR.
Martin Kantor
Diplomová práce
Obsah OBSAH........................................................................................................................................................................... 3 1
ÚVOD .................................................................................................................................................................... 5
2
BEZPEČNOSTNÍ PŘELIVY VODNÍCH DĚL ............................................................................................. 6 2.1
3
4
ZÁKLADNÍ TYPY BEZPEČNOSTNÍCH PŘELIVŮ ............................................................................................... 6
2.1.1
Korunový přeliv ....................................................................................................................................... 6
2.1.2
Boční přeliv.............................................................................................................................................. 8
2.1.3
Šachtový přeliv......................................................................................................................................... 9
2.2
KONZUMPČNÍ KŘIVKA V OBLASTI EXTRÉMNÍCH PRŮTOKŮ .......................................................................10
2.3
POSUZOVÁNÍ BEZPEČNOSTI PŘEHRAD ZA POVODNÍ ...................................................................................12
MATEMATICKÉ MODELOVÁNÍ...............................................................................................................13 3.1
CFD ANALÝZA ............................................................................................................................................13
3.2
ZÁKLADY PROGRAMU FLUENT...................................................................................................................14
3.2.1
Preprocesing ..........................................................................................................................................14
3.2.2
Rovnice proudění a turbulence.............................................................................................................15
3.2.3
Vícefázové proudění ..............................................................................................................................17
3.2.4
Solvery a nastavení výpočtu..................................................................................................................18
3.2.5
Paralelizace výpočtu..............................................................................................................................19
3.2.6
Postprocesing.........................................................................................................................................20
PRAKTICKÁ APLIKACE..............................................................................................................................22 ÚČEL A POPIS VD ORLÍK ............................................................................................................................23
4.1 4.1.1
Korunový přeliv .....................................................................................................................................23
4.2
SITUACE NA VD PŘI POVODNI 2002 [18] ...................................................................................................25
4.3
FYZIKÁLNÍ MODELOVÁNÍ............................................................................................................................26
4.3.1
Stavba fyzikálního modelu.....................................................................................................................27
4.3.2
Měření.....................................................................................................................................................30
4.3.3
Vyhodnocení...........................................................................................................................................33
4.4
CFD MODELOVÁNÍ......................................................................................................................................35
4.4.1
4.4.1.1
Výpočetní síť, okrajové a počáteční podmínky...................................................................................36
4.4.1.2
Zadání a spuštění výpočtu .................................................................................................................37
4.4.1.3
Vyhodnocení ....................................................................................................................................41
4.4.2
4.5
2D model ................................................................................................................................................35
3D model ................................................................................................................................................43
4.4.2.1
Výpočetní síť, okrajové a počáteční podmínky...................................................................................43
4.4.2.2
Zadání a spuštění výpočtu .................................................................................................................46
4.4.2.3
Vyhodnocení ....................................................................................................................................47
ZHODNOCENÍ ...............................................................................................................................................48
-3-
Martin Kantor 5
Diplomová práce
ZÁVĚR................................................................................................................................................................53 5.1
HYDRAULIKA BEZPEČNOSTNÍCH PŘELIVŮ ZA POVODNÍ ............................................................................53
5.2
CFD MODELOVÁNÍ......................................................................................................................................53
5.3
APLIKACE METOD NA VD ORLÍK ...............................................................................................................53
6
POUŽITÉ ZDROJE..........................................................................................................................................54
7
PŘÍLOHY...........................................................................................................................................................55
-4-
Martin Kantor
Diplomová práce
1 Úvod Po povodních, které se udály v nedávné minulosti, došlo k přehodnocení návrhových parametrů vodních děl. Vzhledem ke zvýšení návrhové kapacity bezpečnostních objektů je nutné ověřit kapacity stávajících bezpečnostních přelivů vodních děl. Diplomová práce si klade za cíl zhodnocení hydraulických poměrů vybraných typů bezpečnostních přelivů vodních děl (VD) za extremních průtoků a snaží se upozornit na hlavní hydraulické problémy těchto objektů. Posouzení hydraulických problémů vybraných typů přelivů VD je provedeno několika metodami. Základní metodou posouzení je využití dostupných empirických pozntaků. Jako rozšiřující přístupy jsou použity metody fyzikálního a matematického modelování. V mnoha odvětvích se setkáváme s moderními numerickými metodami (CFD), a proto si tato práce klade za cíl inplementovat metody CFD do oblasti hydraulických výpočtů v hydrotechnické praxi. Celá práce je rozdělena do několika kapitol a jejich řazení za sebou odpovídá postupu řešení diplomové práce. V první kapitole je pojednáno o jednotlivých typech bezpečnostních přelivů a jejich vztahu k bezpečnosti přehrad za povodní. Druhá kapitola popisuje základní možnosti matematického modelování pomocí CFD. Těžiště celé práce je popsáno v poslední kapitole, která se zabývá praktickým využitím fyzikálního a matematického modelování aplikovaného na bezpečnostní přeliv VD Orlík. V závěru jsou shrnuty nejdůležitější poznatky a doporučení.
-5-
Martin Kantor
Diplomová práce
2 Bezpečnostní přelivy vodních děl Bezpečnostní přeliv vodního díla patří k funkčním objektům přehrad. Bezpečnostní přeliv plní funkci pojistného zařízení přehrad, které zabezpečuje přehradní těleso proti přelití za povodně a zajišťuje neškodné převedení povodňového průtoku do koryta pod přehradou. Podle umístění a konstrukčního uspořádání lze rozlišit přelivy [3]: -
korunové, které jsou součástí koruny přehrady a voda se jimi převádí přímo přes přehradní těleso;
-
postranní, vybudované v prodloužení koruny přehrady v boku údolí, přes které voda přepadá souběžně se směrem vodního toku;
-
boční, vybudovaná mimo přehradní těleso v boku nádrže, přes které voda přepadá převážně příčně na směr vodního toku;
-
šachtové, přes které voda z nádrže přepadá do svislé šachty;
-
aj. Dále lze přelivy rozdělit [3]:
-
nehrazené (volné) přelivy, jejichž průtočnost závisí při daném uspořádání jen na stavu hladiny vody v nádrži;
-
hrazené přelivy,jejichž průtočnost lze řídit pomocí uzávěrů. Bezpečnostní přeliv má zásadní vliv na převedení extrémních průtoků za povodní přes
přehradní těleso do koryta pod hrází. Typ přelivu a jeho kapacita ovlivní bezpečnost přehrady jako celku, proto je při návrhu a posouzení nutno přelivu věnovat odpovídající pozornost.
2.1 Základní typy bezpečnostních přelivů V jednotlivých kapitolách jsou popsány základní typy bezpečnostních přelivů, jejich základní hydraulický výpočet a citlivost zařízení z hlediska překročení návrhové kapacity.
2.1.1 Korunový přeliv Korunový přeliv převádí vodu přes přehradní těleso. Přelivná hrana je nejčastěji přímá (betonové tížní přehrady) nebo zakřivená (klenbové přehrady). Přelivná plocha je navržena tak, aby se dosáhlo vysokého součinitele přepadu m při zajištění stability přepadového paprsku. K výpočetu přepadu Q se nejčastěji používá vztah:
Q = mb0 2 gh3/ 2 kde m je součinitel přepadu , platí m = 2/3μ, b0 – účinná šířka přelivu, h – přepadová výška včetně zahrnutí přítokové rychlosti.
-6-
(1)
Martin Kantor
Diplomová práce
Účinná šířka přelivu je celková šířka b přelivných polí, zmenšená o boční kontrakci:
b0 = b − 0.1ξ nk h
(2)
kde nk je počet míst kontrakce, ξ – součinitel závislý na tvaru zhlaví pilíře. Přelivná plocha se konstruuje na principu proudnicové plochy, může být tlaková, beztlaková nebo podtlaková. Tvar proudnicové plochy byl odvozen ze spodní obálky přepadového paprsku Bazinova přelivu. Mezi základní beztlakové plochy patří plocha Scimemiho. Souřadnice přelivné plochy na vzdušní straně přelivu se určí dle: 1.85
x y = 0.5 ha ha
(3)
kde ha je návrhová přepadová výška. Pro návrhovou přepadovou výšku ha má součinitel přepadu Scimemiho beztlakové plochy hodnotu ma=0.504, pro jiné hodnoty přepadové výšky se součinitel určí podle Engeze a Brundenella:
h m = ma ha
0.12
,
platí pro 0.2 ≤ h / ha ≤ 1.3
(4)
Mírně podtlakovou plochou je Smetanova plocha odvozena z plochy Scimemiho. Souřadnice plochy na vzdušní straně přelivu se určí dle: 1.85
x y = 0.461 ha ha
(5)
Návrhová hodnota součinitele přepadu se uvažuje ma=0.5 a pro h < ha platí:
h m = 0.5 0.63 + 0.37 ha
(6)
Při překročení návrhového průtoku se na bezpečnostním přelivu projeví podtlakový režim proudění, přepadový paprsek se přisaje k přelivné ploše a zvýší se součinitel přepadu. Z hlediska stability konstrukce musíme brát zřetel na vznikající podtlaky na přelivné ploše, jež mohou mít vliv na kavitační jevy a sání částeček materiálů z přelivné konstrukce do proudu. Problematická mohou být ruzná přemostění lávky, vyhrazené hradící konstrukce a konstrukce předsazené na bezpečnostním přelivu. Při zachycení proudu o tyto konstrukce dojde ke změně charakteru proudění, půjde o kombinaci přepadu a výtoku velkým otvorem (viz Obr. 1).
-7-
Martin Kantor
Diplomová práce h
Ovlivněný přepad ∆h
Ovlivněný ∆Q1 ∆Q2
přepad Prostý přepad Prostý přepad
Koruna přelivu Q Obr. 1 Korunový přeliv
2.1.2 Boční přeliv Přes boční přeliv přepadá voda převážně kolmo na směr toku do spadiště, z něhož je odváděna rovnoběžně s přelivnou hranou skluzem do prostoru pod přehradou. Přesný výpočet všech jevů souvisejících s bočním přelivem je velmi obtížný, protože charakter proudění je prostorový. Hydraulický výpočet přepadu vody přes boční přeliv je obdobný jako u výpočtu korunového přelivu. Spadiště bočního přelivu se navrhne tak, aby přepad přes přelivné těleso byl i pro největší průtok dokonaly, tj. aby nejvyšší hladina ve spadišti nesahala výše, jak do poloviny přepadového paprsku viz Obr. 2. Pro výpočet platí:
Q = σ mb0 2 gh3/ 2
(7)
kde m je součinitel přepadu, σ – součinitel zatopení, b0 – účinná šířka přelivu, h – přepadová výška včetně zahrnutí přítokové rychlosti. Hydraulický výpočet spadiště bočního přelivu se provádí dle Komorova grafu [2]. Max. hladina h/2 Koruna přelivu
Skluz yk
Spadiště
Obr. 2 Schéma bočního přelivu
-8-
Martin Kantor
Diplomová práce
Na kapacitu bočního přelivu při zatížení extrémním průtokem má vliv zatopení přepadového paprsku vodou ze spadiště (viz Obr. 2), jež se projeví v součiniteli zatopení σ. Dalším omezujícím objektem na přelivu jsou přemostění a lávky (viz Obr. 3), pro jejichž hydraulický výpočet platí rovnice výtoku velkým otvorem:
Q=
(
2 µvb 2 g h23/ 2 − h13/ 2 3
)
(8)
kde μv je součinitel výtoku otvorem, b – šířka otvoru, b2 – hloubka horní hrany otvoru pod hladinou, b1 – hloubka dolní hrany otvoru pod hladinou.
Výtok velkým otvorem
h Zatopený přepad
Přepad
∆h
Nezatopený
∆Q1 ∆Q2
přepad
Koruna přelivu
Q Výtok velkým otvorem Q
Obr. 3 Boční přeliv
2.1.3 Šachtový přeliv Šachtovým přelivem voda přepadá z prostoru nádrže do svislé šachty. Je vhodný pro menší návrhové průtoky a u přehrad, kde je obtížné vytvořit korunový nebo boční přeliv. Šachtový přeliv se skládá z vtokové části, přechodové části, šachty, kolena a odpadní štoly. Pro šachtové přelivy jsou charakteristické dva režimy proudění a to nezahlcený a zahlcený režim viz Obr. 4. Pro nezahlcený režim platí:
Q = mb0 2 gh3/ 2 kde m je součinitel přepadu, b0 – účinná šířka přelivu, h – přepadová výška včetně zahrnutí přítokové rychlosti.
-9-
(9)
Martin Kantor
Diplomová práce
Průtok přelivem při zahlceném režimu zjistíme:
Q=µ
π d2 2 gH 1/ 2 4
(10)
kde d je kritický průměr odpadní šachty, H – tlaková výška vztažena ke kritickému profilu odpadní štoly, μ – výtokový součinitel. Kritérií pro určení hranice mezi zahlceným a nezahlceným režimem je řada od různých autorů. Šachtovými přelivy se zabývali např. Haindl, Wagner, Masiar, Kamenský, Kiselev a další. Při návrhu šachtového přelivu se vychází z předpokladu, že návrhový průtok musím být převeden, aniž by došlo k zahlcení šachty. Dojde-li k zahlcení, šachtový přeliv již není schopen převést větší průtok i když dojde ke zvýšení hladiny v nádrži (viz Obr. 4). Zahlcený režim
h
Nezahlcený režim
Nezahlcený režim
∆h Q
Koruna ∆Q1
přelivu
∆Q2
Zahlcený režim
Q Obr. 4 Šachtový přeliv
2.2 Konzumpční křivka v oblasti extrémních průtoků Stanovení konzumpční křivky bezpečnostního přelivu pro průtoky překračující návrhový může v některých případech být komplikované. Při překročení návrhového průtoku dochází k jevům, jež lze obtížně popsat základními hydraulickými výpočty, proto stanovení konzumpční křivky bezpečnostního přelivu pro průtoky větší jak návrhové, je zatíženo určitou chybou. Abychom minimalizovali velikost chyby, musíme použít vhodný koncepční přístup. Ke stanovení a ověření platnosti konzumpční křivky v oblasti extrémních průtoků máme k dispozici několik metod (Obr. 5): •
matematická extrapolace,
•
hydraulický výpočet,
•
CFD (Computational Fluid Dynamics),
•
fyzikální modelování. - 10 -
Martin Kantor h
Diplomová práce Konzumpční křivka
Známa hladina
Skutečné VD
?
Q=?
Q
Fyzikální model
CFD
Extrapolace
Hydraulický výpočet
Q=f(h)
Obr. 5 Metody stanovení konzumpční křivky v oblasti extrémních průtoků
Matematická extrapolace spočívá v proložení stávající konzumpční křivky polynomem ntého stupně a v následné extrapolaci hodnot do oblasti průtoků větších jak návrhových. Výhodou metody je její nenáročnost a je vhodná v případech, kdy je proudění na přelivu jasně charakterizováno a neovlivněno vnějšímí jevy (proudění okolo překážek jako jsou mostovky, lávky, předsunuté konstrukce a jiné). Hydraulický výpočet spočívá ve výpočtu rovnice přepadu (jak je uvedeno v kapitole 2.1) zohledňující jen základní charakteristiky proudu (součinitel přepadu, součinitel kontrakce, součinitel zatopení). Tato metoda je nenáročná a použitelná u jednoduchých geometrií. CFD (Computational Fluid Dynamics) jako metoda matematického modelování (založená na metodě konečných prvků nebo konečných objemů) se stále víc uplnatňuje v ruzných odvětvích. S rostoucí výpočetní kapacitou stolních počítačů je možné tuto metodu aplikovat pro potřeby hydrotechnických výpočtů. Nutností využití metody je potřeba výpočetního zázemí a znalosti z oboru
- 11 -
Martin Kantor
Diplomová práce
CFD. Náročnost metody se odvíjí od charakteru modelu a použitých omezení (rozměrových, tvarových, popisných). Fyzikální modelování je metodou používanou v minulosti i současnosti k ověření charakteru proudění na bezpečnostním přelivu ve zmenšeném měřítku. Metoda vychází z předpokladu mechanické podobnosti
dvou fyzikálně stejnorodých jevů,a to na prototypu a jeho zmenšeném
modelu, a umožňuje následnou extrapolaci výsledků získaných výzkumem na modelu do skutečnoti. Metoda je náročná na materialní a prostorové zázemí laboratoří.
2.3 Posuzování bezpečnosti přehrad za povodní Posouzením bezpečnosti vodních děl za povodní se zabývá oborová norma TNV 752935, zákon 254/2001 Sb., Metodický pokyn MŽP k posuzování přehrad za povodní (Věstník, duben 1999, ročník IX,částka 4). Tyto předpisy hodnotí bezpečnost vodního díla z pohledu průchodu kontrolní povodňové vlny (KPV), jež může negativně působit na těleso hráze ve smyslu destrukce nebo havarie a následnému vzniku průlomové vlny. Extrémní hydrologický jev je charakterizován KPV, extremita jevu (n-letost opakování) je vztažena k třídě VD (závisí na charakteru možného rizika protržením VD). Předpisy zavádí pojem mezní bezpečné hladiny (MBH), která je definována jako hladina, při niž je v dané lokalitě zaručena bezpečnost a stabilita vodního díla. Vzhledem k extremálnosti jevu a současně krátké době jeho trvání se obecně nevylučuje vznik drobných škod a počítá se sníženým stupněm bezpečnosti. Dále je zaveden pojem kontrolní maximální hladiny (KMH), jež je definována jako úroveň maximální hladiny vody v nádrži při posuzované KPV. Její určení je výstupem vodohospodářské úlohy transformace povodňové vlny retenčním účinkem nádrže za předem přijatých předpokladů a podmínek, které se určí individuálně podle konkrétního konstrukčního typu hráze, použitých funkčních objektů, provozních a dalších souvisejících faktorů. Výrazným faktorem při určení KMH je tedy kapacita přelivu.
- 12 -
Martin Kantor
Diplomová práce
3 Matematické 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), jež se dá charakterizovat jako matematické modelování proudící tekutiny. V současnosti je to oblast, která se výrazně dynamicky rozvíjí v oblasti návrhu virtuálních protoptypů. Dle typu numerického řešení turbulence rozlišujeme:
-
Direct numerical simulation (DNS) – přímá numerická simulace,
-
Large numerical simulation (LES),
-
Reynolds average Navier-Stokes equations (RANS). Jak ukazuje Obr. 6, jednotlivé metody počítají (resolved) vírovou kaskádu až do určité mezní
velikosti ∆xxx, víry menší něž tato velikost jsou již domodelovány. Na metodě závisí také hustota výpočetní sítě, čím menší víry chceme počítat, tím jemnější síť musíme vytvořit a tím se stává výpočet náročnější.
vírová kaskáda
l
η = l/ReL
3/4
Direct numerical simulation (DNS)
Large eddy simulation (LES)
Reynolds averaged Navier-Stokes equations (RANS) Obr. 6 Typy numerického řešení
3.1 CFD analýza Analýza CFD je založena na softwarovém prostředí, které v sobě obsahuje jednotlivé prostředky umožňující: preprocesing (vytvoření geometrie a její vysíťování), vlastní výpočet (nastavení počátečních, okrajových podmínek, spuštění výpočtu) a postprocesing (vyhodnocení výsledku).
- 13 -
Martin Kantor
Diplomová práce
Základní kroky při CFD analýze: -
Definice cílů,
-
Stanovení modelované oblasti,
-
Výběr správného řešiče,
-
Vytvoření výpočetní sítě,
-
Nastavení numerického modelu,
-
Řešení,
-
Zkonvergování řešení,
-
Prohlížení výsledků,
-
Adaptace výpočetní sítě,
-
Revize modelu,
-
Využití výsledku analýzy, jejich interpretace. Trh se softwarem CFD v dnešní době ovládá společnos ANSYS, která disponuje produkty
Fluent a CFX, a tím pokrývá víc jak 60 % trhu s CFD. Mezi dálší produkty patří STAR-CD, FLOW3D a mnoho dalších komerčních, nekomerčních a univerzitních programů. V rámci své diplomové práce jsem se seznámil se softwarem Fluent, protože v rámci ČVUT bylo možno využít kampusovou licenci.
3.2 Základy programu Fluent Fluent je softwarový prostředek umožňující komplexní řešení proudění tekutin, přenosu tepla a hmoty, včetně chemických reakcí, spalování, aerodynamiky a akustiky. Díky jeho multifunkčnímu uplatnění jej lze využít v oblastech energetiky, turbínářství, aerodynamiky, vodohospodářství, automobilovém průmyslu, vzduchotechnice, metalurgii a atd. Vzhledem ke zaměření své práce na aplikaci CFD v hydraulice přelivů jsem se ve Fluentu podrobněji seznámil s problematikou tvorby sítě, okrajovými podmínkami proudění, turbulencí, vícefázovým prouděním ve formě volné hladiny, solverů a jejich nastavení, možností paralelizace výpočtu a možnostmi postprocesingu. Následující kapitoly by měly jednotlivé problémy jednoduše přiblížit.
3.2.1 Preprocesing Základním úkolem je vytvoření výpočetní sítě. K tomuto slouží preprocesor GAMBIT viz Obr. 7 (jež je standartně dodáván k FLUENTU), nebo můžeme využít řadu jiných softwarů (T/Grid, GRIDPRO), popřípadě importovat síť z CAD a jiných CFD programů.
- 14 -
Martin Kantor
Diplomová práce
Obr. 7 Preprocesor GAMBIT
Fluent je založen na metodě konečných objemů a z toho vychází možnosti použití různých výpočetních elementů, ze kterých je vystavěna analyzovaná geometrie. V Obr. 8 jsou uvedeny typy jednotlivých elementů, jak pro 2D tak 3D oblasti, jež FLUENT podporuje. Při tvorbě geometrie můžeme všechny typy elementů vzájemně kombinovat. Pokud je to možné, je výhodné používat qaud(2D) hexa- (3D) elementy, protože při stejné velikosti hrany jsou prostorově úspornější než ostatní elementy a při výpočtu je s nimi dosažena vyšší stabilita,přesnost a rychlejší konvergence. V přápadech složitých geometrií se používají tri- (2D) a tetra- (3D) elementy, pomocí nichž jsme je schopni vysíťovat. Od verze FLUENT 6.3 (prosinec 2006) máme také k dispozici polyhedry (3D), jež svými vlastnostmi tvoří přechod mezi hexa- a tetra- elementy.
Obr. 8 Základní typy elementů
3.2.2 Rovnice proudění a turbulence S laminárním prouděním se v běžné praxi téměř nesetkáme, většina jevů probíhajících kolem nás je charakteristická turbulentním prouděním. Ale z hlediska výpočtu je modelování laminárního proudění o mnoho jednodušší než turbulentního, výpočet je méně náročný na výpočetní čas, rychleji - 15 -
Martin Kantor
Diplomová práce
konverguje a je stabilnější. Z tohoto hlediska je použití laminarního modelu v případech jednoduchých výpočtů dostačující a v případech, kdy potřebujeme hrubě charakterizovat proudění a otestovat okrajové podmínky, nám poslouží pro stanovení počátečních podmínek složitějších výpočtů. V rámci svých výpočtů jsem využíval numerické řešení RANS rovnic pro turbulentní proudění, jde o Reynoldsem upravené Navier-Stockesovy rovnice zavedením středních hodnot parametrů proudění:
ρ
Dui ∂p ∂ ∂ui ∂u j 2 ∂ui ∂ =− + + − δ ij − ρ u 'i u ' j µ + Dt ∂xi ∂xi ∂x j ∂xi 3 ∂xi ∂x j
(
)
(11)
kde Rij = − ρ u 'i u ' j je Reynoldsovo napětí. Modelování turbulence rozdělujeme na modelovaní
v jádře proudu a na modelování
v blízkosti stěny. Modely trbulence jádra proudu založené na RANS rovnici se rozdělují podle počtu doplňujících rovnic na jendorovnicové, dvourovnicové a vícerovnicové. Ve svých výpočtech jsem využil dvourovnicový model k-epsilon, k-omega a sedmirovnicový model Reynolds-Stress (RSM). Samotná rovnice RANS je neřešitelná a musí být doplněna další rovnicí odvozenou Boussinesq hypotézou:
∂u ∂u j 2 ∂u − ρ u 'i u ' j = µt i + − ρ k + µt i δij ∂x ∂xi j ∂xi 3
(12)
kde µt turbulentní viskozita. Dvourovnicové modely turbulence k-epsilon vztahují turbulentní viskozitu k turbulentní kinetické energii (k) a k disipaci energie (ε), k-omega vztahuje turbulentní viskozitu k turbulentní kinetické energii (k) a k disipaci energie (ω). Sedmirovnicový model RSM řeší transportní rovnice přenosu turbulentní energie pro každý člen tenzoru Reynoldsova napětí a je oproti ostatním modelům turbulence robustnější, ale také náročnější na výpočet. V rámci řešení proudění v blízkosti pevné stěny máme několik možností výpočtu turbulence viz. Obr 9. Jako základní se nám nabízí metoda standartní stěnové funkce, která je vhodná v případě hrubé sítě v blízkosti stěny. Tato metoda je založena na modelování proudu v oblasti viskozní a přechodové vrstvy, je doporučeno, aby střed buňky přiléhající ke stěně náležel do inetrvalu 50≤y+≤500, kde y+ je bezrozměrná hodnota vysvětlena rovnicí (13). V případě podrobnějšího řešení proudění v blízkosti stěny přistupujeme k metodě, která je založena na přímém numerickém řešení proudění ve viskozní a přechodové vrstvě, nutností ovšem je zahuštění sítě v blízkosti stěny tak, abychom ve viskozní vrstvě měli alespoň 4 výpočetní buňky tedy y+≈1.
- 16 -
Martin Kantor
Diplomová práce
Přímé řešení
Stěnová funkce
Obr. 9 Možnosti řešení proudění v blízkosti stěny
Závislost mezi rychlostí a vzdáleností od pevné stěny udává Obr. 10. Z obrázku je patrné rozdělení oblastí v blízkosti stěny na jednotlivé vrstvy. Hodnota y+ je definovaná následovně:
y+ =
ū/uτ
ρ uτ ∆y µ
(13)
ū/uτ= 2.5 ln(uτy/v)+5.45
ū/uτ= uτy/v
ln(uτy/ν) Obr. 10 Graf – závislost rychlosti na vzdálenosti od stěny
3.2.3 Vícefázové proudění FLUENT umožňuje pomocí modelu Multiphase Flows řešit vícefázové proudění kapalin, plynů a tuhých látek a jejich vzájemnou interakci. Dle charakteru fází, vzájemného poměru fází a interakcí FLUENT nabízí několk modelů. Pro proudění o volné hladině je speciálně určen model Volume of Fluid (VOF), který řeší vícefázové proudění o velkém vzájemném rozhraní jednotlivých fází.
- 17 -
Martin Kantor
Diplomová práce
Metoda VOF je založena na výpočtu 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 plná dané fáze
0 < αq < 1
výpočetní buňka obsahuje rozhraní fází
Rozhraní fází je řešeno rovnicí kontinuity pro jednotlivé fáze:
∂α q ∂t
+ uj
∂α q ∂xi
=0
(14)
3.2.4 Solvery a nastavení výpočtu Dle řešené úlohy a typu použitého modelu vybíráme ze základních nastavení (segregovaný nebo caplovaný řešič, explicitní nebo implicitní schéma, ustálené nebo neustálené časové řešení). Nastavením
solveru
upravujeme
způsob
numerického
řešení
soustavy
parciálních
diferenciálních rovnic( v Obr. 11 – Equatins). Máme také možnost upravovat relaxační faktory (v Obr. 11 – Under-Relaxation Faktors) a tím rychlost konvergence a stabilitu výpočtu. Stabilitu a rychlost kovergence lze významně ovlivnit vhodným nastavením počátečních a okrajových podmínek a u neustáleného proudění velikostí časového kroku. Fluent, jako řešič založený na metodě konečných objemů, počítá jednotlivé proměnné vždy ve středu buněk a proměnné na hranici elementu dopočítává pomocí extrapolace. K dispozici je několik metod, které se liší přesností a tím i výpočetní náročností (v Obr. 11 – Discretization).
Obr. 11 Nastavení solveru
- 18 -
Martin Kantor
Diplomová práce
Velmi důležitou informací při výpočtu je určení konvergenční meze, tj. kritéria kterým určíme, že výpočet dokonvergoval a průběh simulace je možno považovat za dokončený. Ve FLUENTU nám k tomu slouží průběh residuí (viz. Obr. 12), jež jsou kvantitativním vyjádřením sumy rozdílů jednotlivých proměnných mezi dvěma iteracemi. Jako implicitní hodnota pro residua je nastavena hodnota 1e-03, po dosažení této meze všemi proměnnými je výpočet ukončen (v případě ustáleného proudění), nebo následuje výpočet dalšího časového kroku (u neustáleného proudění). V případě složitějších úloh a neustáleného proudění je však nutno sledovat řadu dalších hodnot, např. sledovat průběh některé z neznámých, integrální síly (odpor, vztlak) nebo integrální bilance toků (hmotnostních, tepelných).
Konvergenční mez
Obr. 12 Graf průběhu residui v průběhu výpočtu
3.2.5 Paralelizace výpočtu Vzhledem k náročnosti výpočtu a možnosti využít výpočetního clusteru na katedře hydrauliky a hydrologie, bylo využito paralelizace výpočtu ve FLUENTU. Princip singl a paralelního výpočtu je znázorněn na Obr. 13. V případě singl výpočtu je celá úloha načtena z pevného disku do operační paměti RAM a poté zpracovávána jedním procesorem (jádrem), důležitým a omezujícím parametrem v tomto případě je velikost operační paměti.
- 19 -
Martin Kantor
Diplomová práce Paralelní výpočet
Singl výpočet
Obr. 13 Schéma singl a paralelního výpočtu
Paralelní výpočet je založen na rozdělení výpočetní domény na více částí (partition) jak ukazuje Obr. 14 a jejich nezávislém zpracování na více procesorech (jádrech). Protože je nutné aby si mezi jednotlivými iteracemi procesory (jádra) předávaly informace o krajových podmínkách mezi jednotlivými partition, je výhodné provést rozdělení tak, aby délka hranice byla co nejkratší (2D), nebo o nejmenší ploše (3D). Komunikace mezi jednotlivými procesory (jádry) je zajišťována protokolem MPI (Message Passing Interface) a to tak, že jeden procesor (jádro) má funci host a rozděluje výpočet na ostatní procesory (jádra) s funkcí nodes. Obdobně jako u singl výpočtu je celá úloha načtena z pevného disku a distribuována do operační paměti, odkud je zpracovávána procesory (jádry). Omezujícími parametry v případě paralelního výpočtu je opět velikost operační paměti a propustnost a kapacita přenosu dat mezi jednotlivými procesory (jádry). Před rozdělením na partition
Po rozdělení na 2 partition
Obr. 14 Schéma rozdělení na partition
3.2.6 Postprocesing K vyhodnocení výsledku simulace postačuje samotný FLUENT, jež obsahuje řadu funkcí pro grafické nebo numerické vyhodnocení. V případě složitějších úloh je možné data vyexportovat do různých formátů a následně je zpracovat v jiných aplikacích (např. MAT-LAB apod). Mezi základní grafické nástroje patří vyhodnocení proudového pole pomocí kontur (vrstevnic), vektorů a trajektorií (viz. Obr 15). Mezi numerické nástroje patří vyhodnocení rovnováhy toků, plošné a objemové itegrály, nebo vyhodnocení silových a momentových účinků. - 20 -
Martin Kantor
Diplomová práce Kontury
Vektory
Obr. 15 Možnosti vizualizace proudění
Po vyhodnocení výsledku simulace a porovnání s daty získanými měřením na skutečném objektu nebo fyzikálním modelu, přistoupíme k revizi modelu. Pokud je ve výsledcích patrný rozdíl, snažíme se najít příčiny těchto rozdílů.
- 21 -
Martin Kantor
Diplomová práce
4 Praktická aplikace Pro praktickou aplikaci metod bylo vybráno vodní dílo Orlík. Při povodni v roce 2002 došlo k překročení návrhové maximální hladiny na VD o více jak 1.5 metru, mající za následek nekontrolovaný odtok z nádrže. Obr. 16 znázorňuje stav na přelivech zhruba při kulminaci povodňové vlny v 8:00 dne 14.8.2002.
Obr. 16 VD Orlík za povodně v roce 2002
Během povodně a bezprostředně po ní vyvstala otázka, jaký byl odtok z nádrže. Vzhledem k neexistenci konzumpční křivky pro oblast tak vysokých průtoků bylo přistoupeno k její extrapolaci autory Brožou a Zezulákem viz [4], [11]. Určení celkového odtoku z nádrže a odtoků přes přelivy, autory Brožou a Zezulákem, se zdá pro daný vodní stav (povoděň 2002) správné. Motivací, proč vybrat toto VD, bylo pro mě zhlédnutí video záznamu [16], na kterém jsem si všiml zvláštního povrchového útvaru v oblasti před přelivy (viz Obr. 16 – zvýrazněné oblasti). Po bližším prozkoumaní vyšlo najevo, že jde o povrchový válec, zachycený konstrukcí mostu. Začalo mě tedy zajímat, jak tato konstrukce může ovlivnit charakter
- 22 -
Martin Kantor
Diplomová práce
proudění na bezpečnostním přelivu VD. Z povahy jevu bylo přistoupeno k modelování fyzikálnímu a matematickému (formou CFD), viz následující kapitoly.
4.1 Účel a popis VD Orlík VD Orlík je nejobjemnější nádrží v ČR a je nejdůležitější částí Vltavské kaskády. Bylo postaveno v letech 1954 až 1961 na Vltavě u Solenic na Příbramsku. Jde o betonovou tížnou přehradu dosahující v koruně výšky 91 m s délkou koruny 450 metrů. Přehrada obsahuje korunový přeliv o třech polích, spodní výpusti, vodní elektrárnu a v pravé části tělesa hráze dvě plavební zařízení. Hlavním účelem VD je především výroba elektrické energie, akumulace vody pro nadlepšení průtoků na spodní části Vltavy a Labe, částečná ochrana území pod přehradou a také Prahy před velkými vodami. Vedle toho slouží nádrž k rekreaci, provozování plavby a vodních sportů a využívá se i k rybímu hospodářství. Základní hydrologické údaje pro přehradní profil VD Orlík: 2
-
plocha povodí
-
průměrný roční úhrn srážek
705
mm
-
průměrný dlouhodobý roční průtok
83,5
m3s-1
12 105,96
km
M-denní průtoky (Qm) m3s-1 M Qm
30
60
90
150
120
180
210
240
270
300
330
355
364
178,0 127,1 100,0 83,4 70,7 60,7 52,4 45,0 38,3 31,9 25,0 17,3 11,3 N-leté průtoky (Qn) m3s-1
N
1
2
Qn
498
688
5
10
20
50
100
966 1190 1140 1770 2050
4.1.1 Korunový přeliv Korunový přeliv o třech polích světlé šířky 15 metrů je hrazen segmentovými uzávěry výšky 8 metrů (viz Obr. 17 - a). Přelivy jsou zakončeny betonovými rozražeči, pod nimiž je umístěn vývar.
- 23 -
Martin Kantor
Diplomová práce
a
b
c
d
Obr. 17 Bezpečnostní přeliv VD Orlík
Pro potřeby pojezdu portálových jeřábů po koruně je před přelivnou hranou umístěna konstrukce mostu (viz Obr. 17, b a c zobrazují detail uložení mostní konstrukce v bočních pilířích, d ukazuje konstrukci samotné mostovky s vodící kolejnicí). Vzhledem k omezené dostupnosti projektové dokumentace k přelivům, bylo nutné pro další práci stanovit přesný tvar přelivné plochy pro potřeby fyzikálního a matematického modelování. Tvar přelivné plochy byl určen dle Scimemiho souřadnic pro beztlakovou plochu. Dle vzorce uvedeného v kapitole 2.1.1 byly spočítány a vykresleny body přelivné plochy pro různé návrhové přepadové výšky. Výsledným porovnáním tvaru vypočteného se skutečným tvarem odečteným ze zapůjčených situací [13], byla návrhová přepadová výška stanovena na hodnotu ha = 9 metrů (viz Obr. 18).
- 24 -
Martin Kantor
Diplomová práce
0
0 1 2 3 4 5 6 7 8 9 10
2
4
6
8
10
12 x (m)
14
odměřeno Scimemi ha=9 m Scimemi ha=9.5 m Scimemi ha=8.5 m y (m)
Obr. 18 Graf – souřadnice přelivné plochy
Přelivná plocha pro souřadnice x<0 je tvořena dvěma oblouky o poloměrech 0.81 a 6.00 metru [13]. Návodní líc má sklon 1:0.066 a vzdušní líc 1:0.728. Situační a výškové uspořádání přelivného pole viz Obr. 19. Koruna hráze 361.00 mn.m.
Měřítko 0
Max. hladina
10 m
354.18
354.60
353.60 mn.m.
Segmentový uzávěr Koruna přelivu
x
345.60 mn.m.
y Obr. 19 Schéma bezpečnostního přelivu VD Orlík
4.2 Situace na VD při povodni 2002 [18] Vlivem nepříznivé metorologické situace v srpnu 2002, postihly převážnou část Jihočeského, Plzeňského a Středočeského kraje dvě srážkové epizody (první 6. až 7. srpna a druhá 11. až 13. srpna), které vyvolaly katastrofální povodňovou situaci.
- 25 -
Martin Kantor
Diplomová práce
Na dolní Vltavě byla první vlna manipulacemi na VD Orlík transformována tak, že v Praze nebyl překročen průtok odpovídající třetímu stupni povodňové aktivity (tj. stavu ohrožení). Na nádrži Orlík se při nástupu druhé povodňové vlny manipulovalo tak, aby byl pozdržen nástup vlny v obcích pod posledním stupněm kaskády a v Praze. Tím byl získán čas k provedení potřebných povodňových opatření (evakuace, stavba protipovodňové stěny). Vzhledem k rychlému nárůstu přítoku se volný prostor nádrže rychle zaplnil a po úplném otevření všech přelivů se odtok z nádrže stal dále neovladatelný. Přítok do nádrže kulminoval 13.8. v poledne na hodnotě zhruba 3900 m3.s-1. P řibližně v té době došlo k havarijnímu přerušení provozu vodní elektrárny a tím ke zmenšení kapacity zařízení odvádějících vodu z nádrže o přibližně 600 m3.s-1. Potom ani kapacita plně otevřených přelivů a výpustí nestačila na převedení kulminujícího přítoku přes hráz a došlo ke stoupnutí hladiny až na kótu 355.17 mn.m., tj. 1.57 m nad maximální povolenou hladinu (viz Obr.20). Maximální odtok z nádrže činil 3100 m3.s-1 (přes přelivy 2950 m3.s-1 [4]), takže kulminace povodňové vlny byla v nádrži Orlík snížena cca o 800 m3.s-1 a zpožděna o 18 hodin.
Hl. povodeň 2002 355.17 mn.m.
354.60
Max. hladina 353.60 mn.m.
Koruna přelivu 345.60 mn.m.
Obr. 20 Bezpečnostní přeliv VD Orlík při povodni 2002
4.3 Fyzikální modelování Fyzikální modelování proběhlo ve školním hydraulickém žlabu ve vodohospodářské hale Fakulty stavební ČVUT v Praze. Základním omezením pro návrh typu modelu je šířka hydraulického žlabu a rozsah dosažitelných průtoků. Na základě těchto omezení byl navržen výsekový 2D model o určitém měřítku zmenšení. Bližší informace k návrhu modelu jsou uvedeny v následujících kapitolách. Před započetím stavby modelu muselo být rozhodnuto, jaké veličiny se budou na modelu měřit, jaké přístrojové vybavení bude použito a jaké doprovodné jevy se budou zkoumat. Vzhledem
- 26 -
Martin Kantor
Diplomová práce
k omezeným znalostem autora v oblasti fyzikálního modelování, bylo rozhodnuto, že se na modelu bude měřit pouze poloha hladiny a bude využito digitální fotografické techniky k zachycení doprovodných jevů. Varianty fyzikálního experimentu: -
geometrie, bude vytvořena základní geometrie (pro zkoumání prostého přepadového paprsku) jež bude v druhé fázi doplněna o model mostní konstrukce, umožňující pojezd portálového jeřábu (umožní zkoumání vlivu mostovky na charakter proudění).
-
řada průtoků (od minimálního, návrhového po maximální ) .
4.3.1 Stavba fyzikálního modelu Vzhledem k omezením, která vyplývají z použití daného hydraulického žlabu a podmínek mechanické podobnosti, bylo stanoveno měřítko modelu M=65 (1:65). Hydraulický žlab, do něhož bude umístěn výsekový model přelivu, je zobrazen na Obr. 21. Oběh vody v hydraulickém žlabu je uzavřený. Voda je čerpána dvěma čerpadly ze zásobní nádrže do potrubí, na kterém jsou umístěny regulační armatury (velké šoupě pro hrubou regulaci a malý ventil pro jemnější regulaci). Z potrubí voda vytéká do uklidňovacího prostoru, kde jsou umístěny uklidňovače a voštinový usměrňovač. Poté již voda protéká skleněným žlabem, který má na začátku výšku 1 m a po 2 metrech je snížen na 0.5 m. Voda před opuštěním žlabu protéká žaluziovým uzávěrem, kterým můžeme regulovat výšku vody ve žlabu, do nádrže, kde přepadá přes Thomsonův přeliv (umožňuje měřit průtok) do zásobní nádrže. BEZPEČNOSTNÍ VOŠTINOVÝ USMĚRŇOVAČ PŘELIV
SCHEMATICKÝ ŘEZ A - A
HROTOVÉ MĚŘÍTKO NA POJEZDU
ADV SONDA NA POJEZDU
ŽALUZIOVÝ UZÁVĚR
1,0 m
B UKLIDŇOVAČ HLADINY
UKLIDŇOVAČE (DĚROVANÝ PLECH)
5,0 m
1,3 m
B
0,5 m
KOLEJNICE POJEZDU
VÁLEC S HROTOVÝM MĚŘÍTKEM
POTRUBÍ S UZÁVĚREM PRO JEMNĚJŠÍ REGULACI PRŮTOKU
SÍTO
OBTOK PRO JEMNĚJŠÍ REGULACI PRŮTOKU
ČERPADLO
ZÁSOBNÍ NÁDRŽ
OBSLUŽNÁ LÁVKA
PŘEPÁŽKA S THOMSONOVÝM PŘELIVEM VOŠTINOVÝ USMĚRŇOVAČ
BEZPEČNOSTNÍ PŘELIV
SCHEMATICKÝ POHLED B - B OBSLUŽNÁ LÁVKA
UKLIDŇOVAČ HLADINY
ADV SONDA NA POJEZDU
VÁLEC S HROTOVÝM MĚŘÍTKEM ŽALUZIOVÝ UZÁVĚR SÍTO
A 0,25 m
ČERPADLA OBSLUŽNÁ LÁVKA
Obr. 21 Schéma hydraulického žlabu
- 27 -
PŘEPÁŽKA S THOMSONOVÝM PŘELIVEM
A
Martin Kantor
Diplomová práce
Omezení vyplývající z použití hydraulického žlabu: -
rozsah dosažitelných průtoků 0 – 44 l.s-1,
-
šířka modelu max. 0.25 m,
-
max. výška modelu s přepadovým paprskem 0.9 m. V tomto případě modelování se vycházelo z Froudova zákona mechanické podobnosti [5],
který vyjadřuje podmínku dynamické podobnosti hydrodynamických jevů za výhradního působení gravitačních sil. Za předpokladu vhodně zvoleného měřítka modelu je vhodný zejména u proudění o volné hladině, např. při modelování proudění na přelivech, vtoků nebo výtoků, větších povrchových vln a při proudění v krátkých úsecích otevřených koryt (vodní skok apod.). Mezní podmínky vyplývající z mechanické podobnosti [5]: -
musí být dosaženo mechanické a dynamické podobnosti,
-
proudění by mělo být v kvadratické oblasti odporů, aby se neuplatnil vliv odporových sil,
-
proudění na modelu by nemělo být ovlivněno kapilárními silami, které způsobují povrchové napětí, musí být splněny následující podmínky:
-
•
přepadová výška musí být h ≥ 50 mm,
•
povrchová rychlost proudu na objektech má být u ≥ 230 mm.s-1,
•
hloubka vodního proudu na modelu musí být h ≥ 15 mm,
šířka přelivného pole má být b0 ≥ 60 mm a šířka hydraulického žlabu, do nichž se umisťuje výsekový model má být b ≥ 200 mm, aby se při výzkumu nepříznivě neuplatnil vliv drsnosti stěn. Na základě těchto omezení a měřítka byl stanoven základní rozměr modelu viz Obr. 22. Výška
modelu je 0.46 m a jeho délka v patě modelu je 0.50 m, model bude zaplňovat celou šířku žlabu, která je 0.25 m. Při stanovení výšky modelu byl brán zřetel na možnost jeho největšího zmenšení, ale za předpokladu výrazného neovlivnění charakteru proudění. Proto jsem se držel předpokladu, že proudnicový tvar přepadového paprsku není ovlivněn, když s/h ≥ 2.8 [9], kde s je v mém případě výška modelu a h je max. přepadová výška. Model bude instalován do žlabu v místě, kde se mění výška žlabu viz Obr. 22, a to proto, že v daném místě je žlab opatřen zpevňujícími svislými výztuhami, které zabrání předpokládanému rozevření žlabu vlivem zvýšeného hydrostatického tlaku (předpokládaný vodní sloupec 0.5 až 0.7 m).
- 28 -
Martin Kantor
Diplomová práce
Max. hladina
x
Min. hladina 0.46 m
y 0.50 m
Obr. 22 Schéma umístění modelu ve žlabu
Model přelivu je vyroben z plastu viz Obr. 23. Základ tvoří tři svislá žebra o tloušťce 6 mm, která jsou ve vodorovném směru spojena v několika úrovních špalíky. Žebra jsou na vzdušní straně potažena plastovou fólií o tloušťce 1 mm a na návodní straně ocelovým plechem o tloušťce 1 mm (z důvodu předpokládaného účinku hydrostatické síly). Dno modelu není uzavřeno. Mostní konstrukce je vyrobena ze dvou kusů plastových špalíků slepených k sobě (rozměr 2.1 x 2.9 cm).
Obr. 23 Pohled na výsekový model přelivu
Vzhledem k předpokládaným velkým silám působícím na model přelivu (max. síla na model ve vodorovném směru bude 53 kg), bylo zvažováno několik možností uchycení modelu ve žlabu. Jako nejednoduší možnost se jevilo uchytit model silikonovým lepidlem, které bude aplikováno po všech okrajích modelu a zároveň bude tvořit těsnění. Tato varianta byla nakonec realizována s předpokladem, že pokud nebude dostačující, provedou se dodatečná opatření. Mostní konstrukce je taktéž vlepena pomocí silikonu do prostoru mezi skly hydraulického žlabu, silové působení se zde nepředpokládá, a tak by silikonové uchycení mělo být dostačující. Souřadnice dolního rohu na vzdušní straně přelivu je x = -0.0431 m, y = -0.1320 m dle souřadného systému v Obr. 22. Osazení bylo provedeno s přesností na 0.1 mm ve svislém směru a 0.5 mm ve vodorovném směru, referenčním bodem byla koruna přelivu.
- 29 -
Martin Kantor
Diplomová práce
Po prvním testovacím spustění žlabu byly zjištěny následující nedostatky: vzhledem k velkému hydrostatickému tlaku došlo k rozevření žlabu a následnému vnikání vody přes silikonem zalepené spáry do vnitřního prostoru modelu, bylo také zjištěno, že voda vniká přes netěsnosti mezi jednotlivými skly žlabu, zdálo se, že vzdušní líc je daleko těsnější než líc návodní a docházelo k tlakování konstrukce modelu zevnitř. Za těchto okolností hrozilo, že dojde buď k nadzvednutí celého modelu vlivem vztlakové síly, nebo k utržení vzdušního líce. Nápravná opatření měla tedy za úkol dotěsnit spáry mezi modelem a sklem tak, aby nedocházelo ke vnikání vody do modelu a pokud by nebylo dotěsnění dostačující, přišlo by do úvahy dodatečné přichycení modelu šrouby ke dnu žlabu. Po opětovném dotěsnění silikonem se ukázalo, že voda stále vniká do modelu. Přistoupil jsem tedy k riskantnímu experimentu otestovat stabilitu uchycení při maximálním průtoku. Konstrukce modelu vydržela opakované maximální možné zatížení (průtok Q ≈ 44 l.s-1) a shledal jsem ji vhodnou k měření.
4.3.2 Měření Jak již bylo uvedeno v úvodu kapitoly, hlavní měřenou veličinou na hydraulickém modelu je poloha hladiny pro řadu průtoků. V prvopočátku se počítalo také s měřením tlaku na přelivné ploše, ale od tohoto záměru bylo upuštěno z důvodu větší složitosti modelu. Základním předpokladem pro měření hladiny bylo nastavení požadovaného průtoku z průtokové řady. Nastavení průtoku se dělo prostřednictvím škrcení na regulačním šoupátku a ventilu viz Obr. 24, průtok se určoval pomocí přepadu přes Thomsonův přeliv viz Obr. 24 a jako měřící nástroj bylo použito hrotové měřítko viz Obr. 24 s možností pojezdu po kolejnicích v horní části žlabu a s možností odečtu hladiny na 0.1 mm. Regulační armatury
Thomsonův přeliv
Hrotové měřítko
Obr. 24 Součásti hydraulického žlabu
Poloha hladiny pro jmenovitý průtok (viz Tabulka 1) byla zaměřena v řadě profilů, hustota profilů byla zvolena tak, aby oblast v okolí přelivné hrany byla zahuštěna. Průběh hladiny pro jednotlivé průtoky je zobrazen na Obr. 25 a Obr. 26. - 30 -
Martin Kantor
Diplomová práce
Nejprve proběhlo měření na modelu bez osazené mostní konstrukce (geometrie G1), poté na modelu doplněném o mostní konstrukci (geometri G2).
Tabulka 1 – Zaměření polohy hladiny pro jednotlivé průtoky Q (l.s-1)
40.4
38.7
36.9
35.7
34.7
34.0
32.9
30.2
27.6
22.7
18.1
14.0
10.3
7.1
G1
X
X
X
X
X
X
X
X
X
X
X
X
X
X
G2
X
X
X
X
X
X
-
-
-
-
-
-
-
-
Poznámka: X měřeno, - neměřeno
Y (cm) 45 40 35 30 25 20 15 10 160
180
200
220
240
260
280
300 X (cm) 320
Obr. 25 Průběh hladiny pro geometrii G1a průtoky 40.4 – 7.1 l.s-1 Y (cm) 45 40 35 30 25 20 15 10 160
180
200
220
240
260
280
300 X (cm)320
Obr. 26 Průběh hladiny pro geometrie G2 a průtoky 40.4 – 34.0 l.s-1
- 31 -
Martin Kantor
Diplomová práce -1
Pro půtoky 7.1 – 34 l.s se přepadový paprsek choval stabilně, hladina významně neoscilovala a polohu hladiny bylo možné zaměřit jedním měřením hrotového měřítka. Hladina přepadového paprsku pro průtoky nad 34 l.s-1 již silně pulzovala a byla zvlněna povrchovými vlnkami (rozměrů 1 – 2 mm) a polohu hladiny bylo nutné určit jako průměr ze tří měření hrotového měřítka. Průtok byl měřen prostřednictvím Thomsonova přelivu. Poloha hladiny v odměrném válci byla měřena hrotovým měřítkem s přesností na 0.1 mm náhodně na začátku, v průběhu a na konci měření. Výsledný průtok byl vypočítan z průměru odečtených hodnot prosřednictvím měrné křivky [17] přelivu. Měření na modelu proběhlo v rozmezí dvou týdnu. Vždy po napuštění modelu vodou (po předcházejícím odvodnění modelu) bylo nutné zavodnit prostor uvnitř vlastního modelu přelivu. Zavodnění probíhalo samovolně vlivem netěsnosti modelu a výhodné bylo model zatížit max. možným průtokem (44 l.s-1), protože se model zavodnil v kratším časovém úseku (cca. 5 minut). V okamžiku, kdy byl model plný vzduchu, docházelo vlivem netěsnosti v oblasti přelivné hrany, k nasávání vzduchu z prostoru uvnitř modelu a vzniku vzduchové kapsy na přelivné ploše viz Obr. 27. Přepadový paprsek se v tomto okamžiku choval jako paprsek volný (tj. nepodepřený přepadovou plochou) s působením atmosférického tlaku na horní i spodní obalovou plochu. V oblasti konce vzduchové kapsy na přelivné ploše docházelo k pulzaci její polohy a stálému odtrhávání vzduchových bublin do proudu. Působením podtlakového efektu vzduchové kapsy došlo k úplnému zavodnění modelu a poté vzduchová kapsa bez dotace vzduchu postupně zmizela a přepadový paprsek se přisál k přelivné ploše (byl zde pozorovatelný efekt zvýšení souč. přepadu vznikem podtlakové plochy). Pulzace
vzduchové
kapsy
Směr proudění Koruna přelivu
Místo úniku vzduchu Obr. 27 Vzduchová kapsa na přelivné ploše
- 32 -
Martin Kantor
Diplomová práce
4.3.3 Vyhodnocení Naměřená data byla vyhodnocena prostřednictvím programu Excel. Výstupem z měření byla hodnota průtoku, pro niž byla zaměřena poloha hladiny v různých profilech, a fotodokumentace doprovodných jevů. Výsledkem měření na fyzikálním modelu je konzumpční křivka přelivu, charakterizovaná prouděním v hydraulickém žlabu šířky 0.25 m bez uvažování boční kontrakce. Konzumpční křivka je zobrazena na Obr. 28, kde měření bez osazené konstrukce jsou ozančena jako prostý přepad a měření s osazenou mostní konstrukcí jsou označena jako ovlivněný přepad. h (m) 0.18 0.16 0.14 0.12
detail
0.1 0.08
prostý přepad
0.06
sestupná větev - ovlivněný přepad
0.04
vzestupná větev - ovlivněný přepad
0.02
Q (m3/s)
0 0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
detail
h (m) 0.18 0.175
prostý přepad sestupná větev - ovlivněný přepad
0.17
vzestupná větev - ovlivněný přepad
0.165 0.16 0.155 0.15 0.145 0.032
0.034
0.036
0.038
0.04
Obr. 28 Konzumpční křivka přelivu pro fyzikální model
- 33 -
Q (m3/s)
0.042
Martin Kantor
Diplomová práce
V průběhu měření a vyhodnocování byl pozorován zvláštní jev související s prouděním v okolí mostní konstrukce. Aby mohl být tento jev lépe poznán, provedla se řada měření pro interval průtoků 33 – 38 l.s-1. Výsledky měření jsou názorně zobrazeny na Obr. 29, z obrázku jsou patrny dvě větve konzumpční křivky v rozmezí průtoků 34 – 37 l.s-1. h (m) 0.18 prostý přepad 0.175
sestupná větev - ovlivněný přepad vzestupná větev - ovlivněný přepad
0.17
IV
III
0.165 0.16 0.155 0.15
II 0.145 0.032
I
0.034
0.036
0.038
I
II
III
IV
0.04
Obr. 29 Konzumpční křivka přelivu pro fyzikální model
- 34 -
Q (m3/s)
0.042
Martin Kantor
Diplomová práce -1
Jedna větev charakterizuje vzestup hladiny od průtoku 34 l.s , fáze I – na návodní straně mostovky se vytváří drobný povrchový válec, povrchový válec se stále zvětšuje až po fázi II, kdy energie přepadového paprsku už nedokáže válec udržet před konstrukcí a dochází náhle k přelití mostní kon. – fáze III.Druhá větev charakterizuje sestup hladiny od průtoku 37 l.s-1, konstrukce je přelévána, jak ukazuje fáze III, a postupně dochází ke zmenšování průtoku, průtok přes konstrukci se stále zmenšuje až do fáze IV, kdy se vytvoří před mostní kon. vodní válec, poté dojde k náhlému pohlcení válce přepadovým paprskem a dostáváme se opět do fáze I.
4.4 CFD modelování Matematické modelování použitím CFD bylo provedeno v softwarovém prostředí Fluent. Pro základní simulaci byl zvolen 2D model a to z důvodu jednodušší tvorby geometrie a menší výpočetní náročnosti. Pro porovnání charakteru proudění a náročnosti výpočtu byl vytvořen také 3D model.
4.4.1 2D model Model byl vytvořen jako výsekový, tj. ve 2D a ve velikosti odpovídající velikosti fyzikálního modelu (M 1:65), a to z důvodu vzájemného porovnání fyzikálního a CFD modelování. Model bude koncipován pro proudění dvou fází tj. vody a vzduchu (model VOF). Varianty numerického řešení (Tabulka 2): -
geometrie, pro stejnou hustotu výpočetních elementů byly vytvořeny dvě 2D geometrie, jedna pro prostý přepad (geometrie G1), druhá s uvažováním vložené překážky pro ovlivněný přepad (geometrie G2),
-
nastavení charakteru proudění, při výpočtu byly použity 4 modely proudění: laminární, turbulentní k-epsilon (k-ε), k-omega (k-ω) a RSM,
-
řada průtoků, zadána formou okrajových podmínek.
- 35 -
Martin Kantor
Diplomová práce
Tabulka 2 – Modelované varianty pro jednotlivé průtoky Q (l.s-1)
G1
G2
38.7
35.7
34.7
34.0
32.9
30.2
27.6
22.7
7.1
lam
-
X
-
-
X
X
X
X
X
k-ε
-
X
-
-
X
X
X
X
X
k-ω
-
-
-
-
X
-
-
-
-
RSM
-
-
-
-
X
-
-
-
-
lam
-
X
X
X
X
-
-
-
-
k-ε
X
X
X
X
X
-
-
-
-
k-ω
-
-
-
-
-
-
-
-
-
RSM
-
X
-
-
-
-
-
-
-
Poznámka: X modelováno, - nemodelováno
4.4.1.1 Výpočetní síť, okrajové a počáteční podmínky Abych si tvorbu sítě co nejvíce zjednodušil, použil jsem trojúhelníkové elementy (tri-). Schéma modelované oblasti je zobrazeno na Obr. 30, velikost přelivu je shodná s velikostí fyzikálního modelu (výška 0.46 m a šířka ve dně 0.5 m). Výpočetní oblast byla rozšířena do míst nad předpokládanou hladinu (prostor pro proudění druhé fáze – vzduchu) a to tak, aby nad přepadovým paprskem byl prostor o velikosti nejméně jednoho přepadového paprsku. Nátok do modelu byl předsazen do vzdálenosti 2.5 metru před přelivnou konstrukcí. Hustota výpočetních elementů byla navržena s předpokládaným uvážením proudění dvou fází (viz Obr. 30), v oblasti předpokládané hladiny (rozhraní dvou fází) byla síť zahuštěna (v okolí přelivu 2 mm, v oblasti nátoku 5 mm). V místě koruny přelivu a okolí na přelivné ploše (viz Obr. 31) je síť zahustěna na velikost 0.5 mm (předpoklad vzniku mezní vrstvy). Velikost elementů v ostatních oblastech je 5, 10 a 20 mm. Celkový počet výpočetních elementů pro geometrii G1 je 250 tis. a pro geometrii G2 je 250 tis.
Obr. 30 Struktura výpočetní sítě pro 2D model se znázorněním hustoty sítě
- 36 -
Martin Kantor
Diplomová práce
Obr. 31 Schéma výpočetní sítě pro 2D model v oblasti koruny přelivu
Okrajové podmínky jsou definovány takto (viz Obr. 30): -
vstup vody a vzduchu do modelu je definován podmínkou velocity – inlet, do úrovně předpokladane hladiny jako vstup vody, nad předpokládanou hladinou jako vstup vzduchu, vstup fáze je charakterizován normálovou rychlostí na vstupní hraně,
-
odtok vody a vzduchu z modelu je definován podmínkou pressure – outlet v pravé části modelu, výstup fází je charakterizován nulovým podtlakem (výtok do volna),
-
model je z horní části ohraničen podmínkou pressure – inlet, která umožňuje vstup a výstup vzduchu do modelu a je charakterizována nulovým přetlakem (atmosférický tlak),
-
dno modelu a přeliv je definován jako pevná stěna, podmínka wall. Počáteční podmínka pro první výpočet byla definovaná jako vodorovná hladina v celém
modelu (viz Obr. 32 – T = 0 s) s horizontální rychlostí vody a vzduchu uvnitř modelu rovnou hodnotě vstupující vody v okrajové podmínce velocity – inlet. Počáteční podmínka pro další výpočty již vycházela z předešlých výpočtů proudění. T=0s
T = 0.2 s
T = 10 s
vzduch
voda
Obr. 32 Rozhraní fází (hladiny) pro časy T
4.4.1.2 Zadání a spuštění výpočtu Výpočet je řešen jako neustálené proudění (nutnost pro VOF) s výpočetním krokem t (s). Nutností u neustáleného proudění je, aby výpočet zkonvergoval v každém časovém kroku, optimální počet iterací na jeden časový krok je kolem deseti. Pokud se počet iterací na jeden časový krok blíží maximálnímu (lze nastavit, obvykle je 20 it.), je nutné zkrátit časový krok t. Jestliže se naopak počet - 37 -
Martin Kantor
Diplomová práce
iterací blíži dvěma, je vhodné časový krok t prodloužit. Optimální nastavení časového kroku t v průběhu výpočtu lze měnit dle charakteru konvergence. Konvergenční kritérium pro jeden časový krok je sledování rezidui, v tomto případě nastavena defaultní hodnota 0.001 pro všechny proměnné. Nalezení konvergenčního kritéria pro celý výpočet, tj. nalezení podmínky po jejímž splnění bylo možno říct, že charakter proudění je ustálený, bylo o něco složitější a bude popsáno dále. Pro očekávané jednodušší řešení byly nejdříve výpočty provedeny laminárním modelem proudění. Poté byly použity i ostatní modely proudění a to modely turbulence k-epsilon, k-omega a Reynolds stress model (RSM). Nejprve byly výpočty provedeny na geometrii (G1) a poté až na geometrii s vloženou mostní konstrukcí (G2). Laminární model – pro rychlejší řešení byly zvoleny základní nastavení solveru (viz kapitola 3.2.4), diskretizace pro tlak – standart, pro hybnost – First Order Upwind, pro objem fází - First Order Upwind. Turbulentí model k-epsilon – pro přesnější postižení proudění byl vybrán složitější model kepsilon Realizable, jež zavádí zpřesnění vložením doplňujících rovnic. Nastavení konstant turbulentního modelu bylo ponecháno na defaultních hodnotách. Nastavení solveru bylo také změněno ve prospěch přesnosti výpočtu (neprospěch časové náročnosti výpočtu). Diskretizace pro tlak – Body Force Weight, pro hybnost, objem fází a turbulentní členy – Second Order Upwind. Turbulentní model k-omega – pro přesnější postižení proudění byl vybrán složitější model komega SST. Nastavení konstant turbulentního modelu bylo ponecháno na defaultních hodnotách. Nastavení solveru bylo také změněno ve prospěch přesnosti výpočtu (neprospěch časové náročnosti výpočtu). Diskretizace pro tlak – Body Force Weight, pro hybnost, objem fází a turbulentní členy – Second Order Upwind. Turbulentní model Reynolds Stress (RSM) byl ponechán v defaultním nastavení. Nastavení solveru bylo také změněno ve prospěch přesnosti výpočtu (neprospěch časové náročnosti výpočtu). Diskretizace pro tlak – Body Force Weight, pro hybnost, objem fází a turbulentní členy – Second Order Upwind. U všech modelů turbulence bylo provedeno standartní nastavení stěnové funkce - Standart Wall Functions, předpokladem pro její použití je 50≤y+≤500. Po vyhodnocení výsledků numerické simulace bylo zjištěno, že pro oblast přelivné plochy se hodnota y+ pohybuje v rozmezí 10 až 20, což napovídá o příliš jemné síti pro tento typ stěnové funkce. Numerické řešení proudění laminárním modelem probíhalo bez větších obtíží. Počáteční nastavení velikosti časového kroku bylo voleno se zřetelem na stabilitu řešení v rozmezí t = 0.0001 až 0.001 s, v průběhu výpočtu bylo možno časový krok prodloužit na t = 0.002 až 0.01 s a to s přihlédnutím ke stabilitě výpočtu. Řešení za ustálené a zkonvergované jsem předpokládal po
- 38 -
Martin Kantor
Diplomová práce
nasimulované době T = 10 až 50 sekund s přihlédnutím k podmínce, že rozdíl přítoku a odtoku vody do modelu je menší jak 1% a průběh tlaku po přelivné ploše odpovídá předpokladům. Numerické řešení proudění turbulentním k-epsilon modelem probíhalo s menšími obtížemi. Počáteční nastavení velikosti časového kroku bylo voleno se zřetelem na stabilitu řešení v rozmezí t = 0.0001 až 0.001 s, v průběhu výpočtu bylo možno časový krok prodloužit na t = 0.001 až 0.005 s, a to s přihlédnutím ke stabilitě výpočtu. Najít kritérium konvergence pro celý výpočet bylo v tomto případě složitější. Porovnání výsledku simulace laminárního a k-epsilon modelu po určité době výpočtu T ≈ 50 s, kdy výsledek výpočtu laminárního modelu se dal předpokládat za zkonvergovaný (ustálený stav), je zobrazeno na Obr. 33. Z obrázku je jasně patrné, že výsledek simulace k-epsilon modelem za ustálený a zkonvergovaný nelze pokládat. Hlavním a podstatným nedostatkem je výskyt oblasti podtlaku uvnitř přepadového paprsku, což je z fyzikálního hlediska nemožné. Fázové rozhraní
Rozložení podtlaku v přepadovém paprsku
Laminární model proudění
Turbulentní k-epsilon model proudění
Obr. 33 Porovnání výsledku řešení pro T≈ 50 s
Po řadě testování se ukázalo za nejefektivnější sledovat průběh statického tlaku v bodě v průběhu výpočtu. Obr. 34 zobrazuje polohu sledovaného bodu, jež je umístěn na návodní straně přelivu v oblasti předpokládaného maxima podtlaku. Druhá část obrázku ukazuje průběh sledované veličiny v čase. Výpočet je podle tohoto grafu možné považovat za zkonvergovaný v čase T = 180 tis. až 190 tis. s, kde je patrné, že statický tlak osciluje kolem hodnoty -1255 Pa.
- 39 -
Martin Kantor
Diplomová práce
Obr. 34 Monitorování statického tlaku v bodě
Numerické řešení proudění turbulentním k-omega a RSM modelem probíhalo obdobně jako modelem k-epsilon. K výpočtu byl využit katedrový cluster skládající se ze 3 stolních PC (3 x P4 3GHz, 1GB RAM, 100 Mbit síť). Z důvodu spuštění paralelního výpočtu bylo potřeba zadanou geometrii rozdělit na 3 partition. Po testování jednotlivých metod rozdělení a zjištění, že pro tuto geometrii jsou rozdíly jen nepatrné, byla vybrána metoda Principial Axes (rozdělení je zobrazeno na Obr. 35).
Obr. 35 Schéma rozdělení na partition
Náročnost výpočtu byla testována pro výpočet modelem k-epsilon. Potřebný čas na provedení jedné iterace se liší podle možnosti využití paralelizace úlohy. V případě výpočtu na jednom procesoru byla doba výpočtu jedné iterace 2.5 s reálného času. Při využití 2 procesorů 1.4 s reálného času a při využití 3 procesorů 1 s reálného času. Uvážíme-li, že výpočet jednoho časového kroku (např. t = 0.001 s modelovaného času) zabere zhruba 10 iterací a ustálené řešení dostaneme např. po namodelování časového úseku T = 100 s, zabere výpočet jediné varianty 12 dnů při zaměstnání 3 procesorů. Tato hodnota je trošku nadsazená, většinou na jeden časový krok stačilo méně jak 10 iterací, nebo byl zvolen delší časový krok t = 0.005 s. Reálně trval výpočet jedné varianty 2 – 5 dnů (při využití 3 procesorů). Náročnost výpočtu laminárniho modelu byla o něco menší, jak u předchozího k-epsilon. Výpočet lépe konvergoval a potřebný čas na jednu iteraci byl také kratší (asi o 10 %). Reálně trval výpočet jedné varianty 1 – 3 dny. Naopak výpočet RSM modelem byl náročnější než výpočet modelem k-epsilon. Výpočet hůře konvergoval a výpočetní čas jedné iterace byl delší (asi o 10 %). - 40 -
Martin Kantor
Diplomová práce
Reálně trval výpočet jedné varianty 5 dnů. Náročnost výpočtu modelem k-omega byla srovnatelná s modelem k-espilon. 4.4.1.3 Vyhodnocení Vyhodnocení výsledků simulací proběhlo v rámci postprocesinku ve Fluentu a následně v tabulkovém procesoru Excel. Hlavní sledovanou veličinou byla poloha hladiny pro jednotlivé průtoky a geometrie. Hodnota výšky hladiny odečítaná pro potřeby určení konzumční křivky byla odečítaná v dostatečné vzdálenosti před přelivem tak, aby byla splněna podmínka L > 6 h (kde L vzdálenost před přelivem, h je výška přepadového paprsku). Konzumpční křivka pro jednotlivé modelované varianty je zobrazena na Obr. 36, kde výpočty bez osazené konstrukce jsou ozančena jako prostý přepad a výpočty s osazenou mostní konstrukcí jsou označeny jako ovlivněný přepad. Mezi výsledky namodelovanými modely k-epsilon, k-omega a RSM nebyly významné rozdíly, a proto v dalších vyhodnoceních budou porovnávány mezi sebou pouze model laminární s modelem turbulence k-epsilon. h (m)
0.18 0.16 0.14 0.12 0.10 0.08
Fluent lam - prostý přepad 0.06
Fluent lam - ovlivněný přepad
0.04
Fluent k-e - prostý přepad Fluent k-e - ovlivněný přepad
0.02 Q (m3/s) 0.00 0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
Obr. 36 Konzumpční křivka přelivu pro matematický model
Další zajímavou veličinou, která nebyla na fyzikálním modelu měřena, ale kterou z matematického modelu můžeme získat, je průběh tlaku po přelivné ploše. Jak zobrazuje Obr. 37, vypočtený průběh tlaku na přelivné ploše je značně rozkolísaný zvláště v oblasti x > 0. Toto rozkolísání je způsobeno tvarem a nespojitostí přelivné plochy, která byla vytvořena spojením úseček. Maximální výkyvy od vyrovnané hodnoty tlaku jsou právě v místě těchto nespojitostí plochy.
- 41 -
Martin Kantor
Diplomová práce 0.2 0.1
x/Ha
0
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
-0.1 -0.2 -0.3
vypočtený tlak
-0.4
vyrovnaná hodnota tlaku
-0.5 -0.6 hp/Ha
Obr. 37 Průběh tlaku po přelivné ploše pro Q = 27.6 l.s-1
Pro porovnání průběhu tlaku po přelivné ploše musely být souřadnice x a tlak přepočteny na bezrozměrnou hodnotu. Porovnání s hodnotami zjištěnými U.S. Army Engineers Waterways Experiment Station [8] je uvedeno v Obr. 38, kde je patrná dobrá shoda s vypočtenými hodnotami zvláště pro oblast x > 0. Rozdíly v oblasti před korunou přelivu jsou způsobeny rozdílným tvarem geometrie (porovnání geometrií zobrazeno pod grafem). V mém případě vypočítaná oblast vysokých podtlaků v místě před korunou přelivu je způsobená velkým zakřivením proudnic na přelivné ploše v těchto místech.
- 42 -
Martin Kantor
Diplomová práce 0.5
0.25 x/Ha 0 -0.4
0
0.4
0.8
1.2
-0.25 H/Ha=0.5 [8] H/Ha=1 [8]
-0.5
H/Ha=1.33 [8] Q=35.7 l.s-1 (H/Ha=1.12) -0.75
Q=27.6 l.s-1 (H/Ha=0.96) Q=7.1 l.s-1 (H/Ha=0.43)
-1
hp/Ha
model
x
AEWES [8]
y geometrie
Obr. 38 Průběh tlaku po přelivné ploše
4.4.2 3D model Pro potřeby porovnání výsledků a otestování složitosti výpočtu byl vytvořen také 3D model geometrie přelivu VD Orlík. Model byl vytvořen v měřítku M 1:1. Varianty numerického řešení: -
3D geometrie (geometrie přelivu s mostní konstrukcí),
-
nastavení charakteru proudění, při výpočtu byly použity 2 modely proudění: laminární a turbulentní k-epsilon,
-
jeden průtok zadaný formou okrajových podmínek. 4.4.2.1 Výpočetní síť, okrajové a počáteční podmínky Velikostí se řadí přeliv VD Orlík k jedněm z největších v ČR, rozměr každého ze tří polí je 15
m na šířku s přepadovou výškou cca. 9 m. Celý přeliv je schopen převest při maximálním zatížení cca. 3000 m3.s-1. Z těchto informací je zřejmé, že pokud by bylo možno matematický model patřičně zjednodušit, bylo by to výhodné. Rozhodl jsem se pro zjednodušenou formu řešení, a to pouze jedné poloviny přelivu (viz Obr. 40). Model obsahuje polovinu vnitřního pole, jeden vnitřní pilíř, celé krajní - 43 -
Martin Kantor
Diplomová práce
pole a boční pilíř. Model byl ve všech směrech rozměrově ohraničen za podmínky, že ohraničením nedojde k výraznému ovlivnění proudění. Celkový rozměr modelované oblasti je 45 x 85 x 45 m (šířka x délka x výška) o objemu 95 tis. m3. Celá geometrie byla vysíťována v preprocesoru Gambit. V počátku kdy jsem neměl představu o náročnosti 3D výpočtu, jsem vytvořil jednoduchou výpočetní síť tvořenou tetra- prvky (Obr. 39 – geometrie GA). Vycházel jsem z předpokladu, že síť bude v okolí přelivu zahuštěna zhruba na velikost hrany elementu 0.3 metru a všemi směry se síť bude postupně zjemňovat až na rozměr 2 metry pro elementy v oblasti nátoku u dna modelu. Celkový počet elementů této sítě se pohyboval kolem cca. 600 tis. Síť byla importována do Fluentu a proběhl první výpočet, který naznačil, že takovýto typ sítě je zcela nepoužitelný, protože výpočet nebyl schopen zkonvergovat již v prvním časovém kroku. Poté jsem přistoupil k tvorbě modifikované sítě (Obr. 39 – geometrie GB) tvořené kombinací hexa- a tetraprvky. V oblastech, kde to bylo možné, byla vytvořená síť tvořena hexa- prvky, pouze oblast v těsném okolí přelivné hrany, pilířů a mostní konstrukce byla kvůli geometrické složitosti vysíťována tetraprvky. Velikost hrany tetra- prvku byla zvolena 0.2 m pro oblasti u stěn a postupně se velikost zvětšovala až na 0.4 m pro oblasti přechodu v síť tvořenou hexa- prvky. Velikost hrany hexa- prvků se pohybovala v rozmezí 0.3 až 2 metry. Celkový počet všech elementů byl cca. 730 tis. S příchodem nové verze Fluentu s označením 6.3 (v prosinci 2006) jsem využil možnosti výpočtu s prvky tvaru polyheadr. Přechod na polyheadry spočívá v konverzi prvku tetra- v poly- prvky, jak je zobrazeno na Obr. 39, a tak byla vytvořena geometrie GC z geometrie GB. Celkový počet elementů se snížil na 630 tis. Snížením počtu elementů se také snížila výpočetní náročnost.
- 44 -
Martin Kantor
Diplomová práce
GA tetra-prvky
GB tetra- a hexa-prvky
GC hexa- a poly- prvky
Obr. 39 Jednotlivé varianty výpočetní sítě - detail
Okrajové podmínky jsou definovány takto (viz Obr. 40): -
vstup vody a vzduchu do modelu je definován podmínkou velocity – inlet, do úrovně předpokládané hladiny jako vstup vody, nad předpokládanou hladinou jako vstup vzduchu, vstup fáze je charakterizován normálovou rychlostí na vstupní ploše,
-
odtok vody a vzduchu z modelu je definován podmínkou pressure – outlet, výstup fází je charakterizován nulovým podtlakem (výtok do volna),
-
model je z horní části ohraničen podmínkou pressure – inlet, která umožňuje vstup a výstup vzduchu do modelu a je charakterizována nulovým přetlakem (atmosférický tlak),
-
dno modelu a přeliv je definován jako pevná stěna, podmínka wall.
-
osa symetrie je definovaná podmínkou symetry, která je charakterizovaná jako stěna s nulovou drsností.
- 45 -
Martin Kantor
Diplomová práce
pressure - outlet
symetry
wall
velocity - inlet
Obr. 40 Okrajové podmínky
Počáteční podmínka pro výpočet byla definovaná jako vodorovná hladina v celém modelu s horizontální rychlostí vody a vzduchu uvnitř modelu rovnou hodnotě vstupující vody v okrajové podmínce velocity – inlet. 4.4.2.2 Zadání a spuštění výpočtu Výpočet je řešen jako neustálené proudění (nutnost pro VOF) s výpočetním krokem t (s). Pro výpočet bylo využito geometrie GB (tvořené tetra- a hexa- prvky) a jediný průtok Q = 2950 m 3.s-1, který odpovídá průtoku převedeném přelivy při povodni v roce 2002. Pro očekávané jednodušší řešení byly nejdříve výpočty provedeny laminárním modelem proudění. Poté byl použit model turbulence k-epsilon. Laminární model – pro rychlejší řešení bylo zvoleno základní nastavení solveru (viz kapitola 3.2.4), diskretizace pro tlak – standart, pro hybnost – First Order Upwind, pro objem fází - First Order Upwind. Turbulentí model k-epsilon – pro přesnější postižení proudění byl vybrán složitější model kepsilon Realizable, jež zavádí zpřesnění vložením doplňujících rovnic. Nastavení konstant turbulentního modelu bylo ponecháno na defaultních hodnotách. Nastavení solveru bylo také změněno ve prospěch přesnosti výpočtu (neprospěch časové náročnosti výpočtu). Diskretizace pro tlak – Body Force Weight, pro hybnost, objem fází a turbulentní členy – Second Order Upwind. U modelu turbulence k-epsilon bylo provedeno standartní nastavení stěnové funkce - Standart Wall Functions. Numerické řešení proudění laminárního i turbulentního modelu probíhalo bez větších obtíží. Počáteční nastavení velikosti časového kroku bylo voleno se zřetelem na stabilitu řešení v rozmezí t =
- 46 -
Martin Kantor
Diplomová práce
0.0001 až 0.001 s, v průběhu výpočtu bylo možno časový krok prodloužit na t = 0.002 až 0.01 s, a to s přihlednutím ke stabilitě výpočtu. Řešení za ustálené a zkonvergované jsem předpokládal po splnění konvergenčních kritérií (Obr. 41), jedním bylo ustálení průběhu tlaku ve sledovaném bodě na koruně přelivu, druhým kritériem bylo vyrovnání hmotnostního vtoku a výtoku do a z modelu.
Obr. 41 Graf – konvergenčních kritérií
K výpočtu byl využit katedrový cluster skládající se ze 3 stolních PC (3 x P4 3GHz, 1GB RAM, 100 Mbit síť). Z důvodu spuštění paralelního výpočtu bylo potřeba zadanou geometrii rozdělit na 3 partition podle metody Cartesian X-Coordinate. Náročnost výpočtu byla testována pro výpočet modelem laminárním. Potřebný čas na provedení jedné iterace se liší podle možnosti využití paralelizace úlohy. V případě výpočtu na jednom procesoru byla úloha neřešitelná (nedostatečná velikost RAM paměti). Při využití 2 procesorů byla doba výpočtu jedné iterace 9 s reálného času a při využití 3 procesorů 7.2 s reálného času. Reálně trval výpočet jedné varianty 1 – 2 týdny (při využití 3 procesorů). Náročnost výpočtu turbulentním modelem k-epsilon byla oproti modelu laminárnímu o něco náročnější (asi o 10 %). 4.4.2.3 Vyhodnocení Vyhodnocením výsledku simulace byla zjištěna poloha hladiny odpovídající zadanému průtoku. Určení této hladiny nebylo jednoznačné, protože se zde projevuje prostorový jev snížení hladiny směrem z nádrže k přelivu. Jako reprezentativní hodnota hladiny byla vybrána maximální hodnota v oblasti nátoku do modelu. Pro průtok 2950 m3.s-1 odpovídá namodelovaná hodnota přepadového paprsku h = 10.02 m (model k-epsilon). Laminární model dával obdobné výsledky. Využitím 3D modelování jsme byli schopni postihnout prostorové jevy proudění. Obr. 42 např. znázorňuje trajektorie proudící vody pro zkoumaný průtok. Dále je možno získat obdobné výstupy jako byly prezentovány u 2D modelování (např. průběhy tlaků a rychlostí).
- 47 -
Martin Kantor
Diplomová práce
Obr. 42 Trajektorie
Na Obr. 43 jsou znázorněny možnosti vizualizace proudění, jako je zobrazení fázového rozhraní pro hodnotu 0.5, trajektorie a povrchové rychlosti na hladině. Fázové rozhraní (hladina)
Trajektorie
Rychlost na hladině
Obr. 43 Možnosti vizualizace proudění
4.5 Zhodnocení Porovnání výsledků z fyzikálního a numerického řešení 2D na zmenšeném modelu je zobrazeno formou kozumpční křivky na Obr. 44. S fyzikálním modelem dobře koresponduje turbulentní model k-espilon. Laminární model se s fyzikálním modelem shoduje pouze v oblasti
- 48 -
Martin Kantor
Diplomová práce
prostého přepadu. V tomto případě je to způsobeno jednodušší diskretizací u laminárního modelu, která je při použití tri- elementů nedostačující. h (m)
0.18 0.16 0.14 0.12
detail
0.1 0.08
fyzikální modelování
0.06
CFD laminární model
0.04
CFD k-epsilon model
0.02
Q (m3/s)
0 0
0.18
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
detail
h (m) fyzikální modelování
0.175
CFD laminární model 0.17
CFD k-espilon model
0.165 0.16 0.155 0.15 0.145 0.032
0.034
0.036
0.038
0.04
Q (m3/s)
0.042
Obr. 44 Konzumpční křivka zmenšeného modelu přelivu
Porovnání průběhu hladiny pro jeden průtok je zobrazeno na Obr. 45, z obrázku je patrná shoda naměřených dat na fyzikálním modelu s daty vypočítanými turbulentním modelem k-espilon. Laminární model podhodnocuje součinitel přepadu a výsledná hladina je výše, což také způsobuje rozdíly, které byly patrné v konzumpční křivce uvedené v předchozím obrázku.
- 49 -
Martin Kantor
Diplomová práce
detail
detail
Obr. 45 Porovnání polohy hladiny pro Q = 32.9 l.s-1
Dobré shody pro fyzikální model a CFD model bylo také dosaženo v oblasti přelití mostní konstrukce (Obr. 46). Jev popsaný v kapitole 4.3.3 , související s hysterzí konzumpční křivky, pozorované na fyzikálním modelu pro průtoky 34 až 38 l.s-1, nebyl na výsledku numerického výpočtu pozorován. Tento jev, způsobený pravděpodobně povrchovým napětím, nebyl v numerickém řešení zohledněn z několika důvodů: hustota výpočetní sítě v oblasti mostní konstrukce není dostatečně zahuštěna a v simulaci nebylo počítáno s vlivem povrchového napětí na přepadový paprsek.
- 50 -
Martin Kantor
Diplomová práce
Fyzikální model
CFD modelování
Obr. 46 Porovnání charakteru proudění pro Q = 38.1 l.s-1
Předpokládaným výstupem jak fyzikálního, tak CFD modelování by měla být konzumpční křivka skutečného přelivu na VD Orlík. Při přepočtu z modelu (fyzikálního nebo CFD 2D) na skutečnost musíme zohlednit řadu jevů souvisejících s prostorovým prouděním na skutečném objektu. Přepočet přepadové výšky z modelu do skutečnosti byl proveden pouze formou měřítka zmenšení. S přepočtem přepadového průtoku to bylo již složitější, musel být zaveden opravný součinitel zohledňující vliv kontrakce pilířů na přepadový paprsek. Tento součinitel jsem získal úpravou rovnice (2), s použitím součinitele tvaru pilířů ξ v rozmezí hodnot 0.2 až 0.3, podle tloušťky přepadového paprsku (ξ = 0.2 pro h ≈ 6 m až ξ = 0.3 pro h ≈ 10 m). V případě využití výsledku CFD simulace 3D modelem nebylo nutno zavádět opravné součinitele, protože samotný model již postihuje prostorový charakter proudění s uvážením bočního zúžení přepadového paprsku. Konzumpční křivka pro jedno pole bezpečnostního přelivu VD Orlík je zobrazena na Obr. 47. Stávající konzumpční křivka zjištěná z manipulačního řádu [12] je stanovena pro rozsah přepadového paprsku od 0 po 8 m, nad tuto hranici je provedena extrapolace [11] do tloušťky paprsku 9.6 m. Nad hranicí 9.6 m může být přepad vody ovlivněn mostní konstrukcí, a proto konzumpční křivka nad touto hranicí bude určena na základě modelování fyzikálního a numerického. Hodnoty reprezentované CFD modelováním (2D i 3D) jsou výsledkem řešení turbulentního modelu k-espilon (laminární model kvůli nepřesnostem nebyl pro přepočet do skutečnosti vůbec uvažován).
- 51 -
Martin Kantor
Diplomová práce
h (m) 12.00
10.00
8.00
6.00
konzumpční křivka [12] extrapolace [11]
4.00
fyzikální model CFD modelování 2D
2.00
CFD modelování 3D
0.00 0
200
400
600
800
1000
1200 Q (m3/s)
Obr. 47 Konzumpční křivka jednoho pole přelivu na VD Orlík
Z konzumpční křivky je patrné, že výsledky fyzikálního a numerického řešení se pro rozsah přepadového paprsku od 0 po 9.6 m ztotožňují se stávající konzumpční křivkou [12] a její extrapolací [11]. Stanovit konzumpční křivku nad hranici 9.6 m je již problematické a reálně neexistují data s nimiž bychom ji mohli srovnat. Jev pozorovaný při fyzikálním modelování související s hysterzí konzumpční křivky v oblasti průtoků kolem 1000 m3.s-1 se ve skutečném měřítku téměř neprojeví, protože tento jev na fyzikálním modelu pravděpodobně souvisel s povrchovým napětím. Proto se dá očekávat, že konzumpční křivka bude nad hranicí 9.6 m pro přepadový parsek procházet body určenými CFD modelováním (2D) a zároveň bude procházet zprůměrovanou konzumpční křivkou určenou fyzikálním modelováním (výsledná konzumpční křivka bezpečnostního přelivu VD Orlík je zobrazena v příloze 1 a 2). Hodnota reprezentovaná CFD modelováním ve 3D je nepatrně nadsazená, to je způsobeno nedostatečnou hustotou výpočetní sítě v oblasti mostní konstrukce (hrana elementu o velikosti 0.2 m) a v okolí přelivné hrany. Získání přesnějšího výsledku 3D simulace by si vyžádalo zahuštění výpočetních elementů a tomu odpovídající výpočetní kapacity.
- 52 -
Martin Kantor
Diplomová práce
5 Závěr 5.1 Hydraulika bezpečnostních přelivů za povodní Vzhledem k rozsahu diplomové práce nemohlo být zkoumáno více typů bezpečnostních přelivů. Práce prokázala, že problematika stanovení konzumpční křivky (kapacity) bezpečnostních přelivů pro oblasti extrémních průtoků je velmi složitá. Zvláštní pozornost tedy zasluhují především boční a šachtové přelivy.
5.2 CFD modelování Práce si kladla za cíl využití metody CFD analýzy na praktickém příkladě. V tomto bodě nabídla pohled na přístup k dané problematice a jednoduchý postup. Hlavním poznatkem, který bych chtěl zdůraznit, je: metoda CFD analýzy je použitelná v těch případech, kdy máme k dispozici srovnavací data, ať z fyzikálního modelování nebo ze zkutečného provozu. V případech kdy nemáme dostatečné množství dat pro srovnavání, nás výsledky samotného CFD modelování nemohou nikdy dostatečně uspokojit!
5.3 Aplikace metod na VD Orlík Na příkladu bezpečnostního přelivu VD Orlík bylo ukázano použití jednotlivých přístupů, jež by měli vést ke stanovení nebo ověření konzumpční křivky v oblasti extrémních průtoků. Autor se převážně zabýval metodami fyzikálního a numerického modelování. Výsledky experimentu prokázaly, že geometrie mostní konstrukce (předsazená před přelivem) může významně ovlivnit charakter proudění na bezpečnostním přelivu VD. Významné ovlivnění lze očekývat pro průtoky větší jak 1100 m3.s-1 pro jedno pole bezpečnostního přelivu. Pokud se ukáže z praktického pohledu provozovatele vodního díla, že tento problém je dále nutné zkoumat, doporučuje se vytvoření 3D fyzikálního modelu (předpoklad postihnutí výrazného 3D proudění na přelivu VD). Výsledná konzumpční křivka bezpečnostního přelivu VD Orlík je v tabelární formě uvedena v příloze 1. Příloha 2 pak obsahuje její grafickou podobu pro jedno jezové pole.
- 53 -
Martin Kantor
Diplomová práce
6 Použité zdroje [1]
Broža, V., (2005) Přehrady Čech, Moravy a Slezska, Liberec
[2]
Broža, V., Satrapa, L., (2000) Hydrotechnické stavby 10 – Přehrady, Praha
[3]
Broža, V., Kratochvíl, J., Peter, P., Votruba, L., (1987) Přehrady, Praha
[4]
Broža, V., (2002) Výpočet přítoku do nádrže Orlík za mimořádné povodně v srpnu 2002, Praha
[5]
Čábelka, J., Gabriel, P., (1987) Matematické a fyzikální modelování v hydrotechnice 1, Praha
[6]
Fluent Inc., FLUENT User’s Guide
[7]
Fluent Inc., Gambit User’s Guide
[8]
Chow, V. T., (1959) Open-Channel Hydraulics, McGraw Hill
[9]
Kolář, V., Patočka, C., Bém, J., (1983) Hydraulika, Praha
[10]
Králík, M., (2005) Boční přeliv a bezpečnost přehrad, Praha
[11]
Krejčí, J., Zezulák, J., (2003) Vyhodnocení povodně v srpnu 2002 z pohledu průchodu povodňové vlny Vltavskou kaskádou, Praha
[12]
Povodí Vltavy, státní podnik, (2002) Manipulační řád pro vodní dílo Orlík na Vltavě
[13]
Povodí Vltavy, státní podnik, Dokumentace k vodnímu dílu Orlík
[14]
Tesař, V., (1996) Mezní vrstvy a turbulence, Praha
[15]
Metodický pokyn odboru ochrany vod MŽP ČR k posuzování přehrad za povodní, In: Věstník MŽP, duben 1999, roč. IX, částka 4
[16]
Video záznam, (2002) Velká voda na Vltavské kaskádě 9-14.8.2002
[17]
Konzumpční křivka Thomsonova přelivu pro studentský žlab
[18]
Výsledná zpráva o projektu Vyhodnocení katastrofální povodně v roce 2002 a návrhu úpravy systému prevence před povodněmi, (2004) Výzkumný ústav vodohospodářský T.G.Masaryka
- 54 -
Martin Kantor
Diplomová práce
7 Přílohy Příloha č. 1 Tabelární konzumpční křivka bezpečnostního přelivu VD Orlík Příloha č.2 Graf konzumpční křivky bezpečnostního přelivu VD Orlík
- 55 -
Martin Kantor
Diplomová práce
Příloha č.1 Tabelární konzumpční křivka bezpečnostního přelivu VD Orlík
H [mn.m.]
h [m]
Q1pole [mn.m.]
H [mn.m.]
h [m]
Q1pole [mn.m.]
345.60 345.80 346.00 346.20 346.40 346.60 346.80 347.00 347.20 347.40 347.60 347.80 348.00 348.20 348.40 348.60 348.80 349.00 349.20 349.40 349.60 349.80 350.00 350.20 350.40 350.60 350.80 351.00 351.20
0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00 3.20 3.40 3.60 3.80 4.00 4.20 4.40 4.60 4.80 5.00 5.20 5.40 5.60
0 2.04 5.98 11.3 17.7 25.2 33.6 42.9 53.1 64.2 76 88.5 102 116 131 146 162 179 196 214 233 252 272 293 314 336 358 381 404
351.60 351.80 352.00 352.20 352.40 352.60 352.80 353.00 353.20 353.40 353.60 353.80 354.00 354.20 354.40 354.60 354.80 355.00 355.20 355.40 355.60 355.80 356.00 356.20 356.40 356.60 356.80 357.00 357.20
6.00 6.20 6.40 6.60 6.80 7.00 7.20 7.40 7.60 7.80 8.00 8.20 8.40 8.60 8.80 9.00 9.20 9.40 9.60 9.80 10.00 10.20 10.40 10.60 10.80 11.00 11.20 11.40 11.60
453 478 504 530 557 584 611 640 669 698 728 758 789 820 852 884 917 950 984 1020 1048 1079 1094 1107 1120 1130 1151 1175 1200
351.40
5.80
428
- 56 -
Martin Kantor
Diplomová práce
Příloha č.2 Graf konzumpční křivky bezpečnostního přelivu VD Orlík
průtok Q jedním přelivným polem (m3.s-1)
úroveň hladiny (mn.m.)
přepadová výška H (m)
- 57 -
