TECHNICKÁ UNIVERZITA V LIBERCI
Fakulta mechatroniky, informatiky a mezioborových studií
Porovnávání numerických modelů podzemního proudění v rozpukaných horninách
Diplomová práce
2012 Prokop Žižka
TECHNICKÁ UNIVERZITA V LIBERCI
Fakulta mechatroniky, informatiky a mezioborových studií
Porovnávání numerických modelů podzemního proudění v rozpukaných horninách
Prokop Žižka
Studijní program: B 2612 – Elektrotechnika a informatika Studijní obor: Automatické řízení a inženýrská
informatika
Pracoviště: Ústav nových technologií a aplikované informatiky
Fakulta mechatroniky, informatiky a mezioborových studií
Technická univerzita v Liberci Studentská 2, 461 17 Liberec Školitel: doc. Ing. Otto Severýn, Ph.D.
Rozsah diplomové práce Počet stran:
Počet obrázků:
Sem přijde zadání od školy (jen info)
Prohlášení:
Byl jsem seznámen s tím, že na mou bakalářskou práci se plně vztahuje zákon č.
121/2000 o právu autorském, zejména § 60 (školní dílo).
Beru na vědomí, že Technická univerzita v Liberci má právo na uzavření licenční smlouvy o užití mé práce a prohlašuji, že souhlasím s případným užitím mé bakalářské práce (prodej, zapůjčení apod.).
Jsem si vědom toho, že užití své bakalářské práce či poskytnout licenci k jejímu využití mohu jen se souhlasem Technické univerzity v Liberci, která má právo požadovat ode mne přiměřený příspěvek na úhradu nákladům vynaložených univerzitou na vytvoření díla (až do jejich skutečné výše).
Bakalářskou práci jsem vypracoval samostatně s použití uvedené literatury, pod vedením školitele.
V Liberci dne 20. 12. 2011 ………....
Abstrakt:
Tato práce se zabývá srovnáváním vlastností numerických modelů proudění podzemních vod v rozpukaném prostředí. Nejprve provedeme stručný výklad základních pojmů hydrogeologie a matematického popisu podzemního proudění. Další část je věnujeme simulačním programům Flow123D a FEFLOW. Následuje podrobný popis zadání shodné úlohy do obou programů, abychom mohly porovnávat výsledky testů. Dále se seznámíme s opensource programem ParaView, který slouží pro analýzu a vizualizaci dat. Ten využívat k porovnání výsledků námi nasimulovaných testů.
Abychom mohli srovnávat a zobrazovat výsledky z obou programů (Flow123D a FEFLOW), zvolíme si strukturu datového souboru typu VTK. O struktuře pojednává další část práce, která kromě popsání složení dat, ukáže i zacházení s jednoduchým programem na převod výsledů z FEFLOW do souboru VTK. V poslední části se věnujeme řadě ukázkových testů, na kterých demonstrujeme výsledky této práce.
Klíčová slova: Flow123D, FEFLOW, ParaView, proudění podzemní vody, srovnání numerických modelů, metoda konečných prvků
Abstract:
The diploma thesis deals with the comparison of the characteristics of numerical models of subterranean water circulation in cracked surface. First of all, a brief explanation of fundamental concepts of hydrogeology and mathematical specification of subterranean circulation is presented. The following section of the thesis focuses on the simulation programmes Flow123D and FEFLOW. This chapter gives details on these programmes, their control and environment. The description of identical assignment instructions for both programmes is provided. Afterward, open source programme ParaView, which enables an analysis and visualisation of data, is introduced. This programme is used for comparing simulated test results and their graphical representation. In order to enable comparison and projection of the results from both programmes (Flow123D a FEFLOW), the procedure of the format transmission to the structure of the data file, type VTK, has been chosen. This data organization is described in chapter five. For the transmission of the results from FEFLOW to the VTK file, a simple transfer programme, which has been created to smooth lengthy manual work, is used. Finally, on the basis of a range of sample tests that are presented in the last section of the diploma thesis, the results of the investigation are evaluated and displayed.
Key words: Flow123D, FELOW, ParaView, groundwater flow, comparison of numerical models, finite element method
8
Obsah:
Seznam použitých symbolů ... 11
Seznam použitých zkratek ... 12
Seznam obrázků ... 13
Úvod ... 15
1 Základní pojmy hydrogeologie ... 16
1.1Podzemní vody ... 16
1.2 Zvodeň ... 18
1.3 Porozita ... 19
1.4 Propustnost ... 20
2 Proudění podzemní vody ... 21
2.2 Omezení Darcyho zákona ... 24
2.3 Okrajové podmínky ... 25
2.4 Proudění v rozpukaném prostředí ... 26
3 Metoda konečných prvků ... 28
4. Simulační programy ... 29
4.1 Flow123D ... 29
4.1.1. Struktury souborů pro Flow123D ... 30
4.1.1.1 Síť *.msh ... 30
4.1.1.2 Materiály *.mtr ... 32
4.1.1.3 Okrajové podmínky *.bcd ... 32
4.1.1.4 Sousednosti *.ngh ... 34
4.1.1.5 Zdroje *.src ... 34
4.2 FEFLOW ... 35
4.2.1 Modelování ve FEFLOW ... 36
4.2.1.1 Síť ve FEFLOW ... 36
4.2.1.2 Propustnost FEFLOW ... 36
4.2.1.3 Okrajové podmínky Feflow ... 37
4.2.1.4 Zdroje FEFLOW ... 37
9
4.2.1.5 Samotná simulace FEFLOW ... 37
5 Způsob srovnávání výsledků ... 38
5.1Struktura VTK... 39
5.2Převod z FEFLOW do VTK ... 41
5.3 ParaView ... 42
5.3.1Zobrazení v ParaView ... 42
5.4 Postup porovnání v ParaView ... 43
6 Nastavení stejné úlohy v obou systémech ... 45
6.1 Zadání příkladu ... 45
6.1 Síť ... 46
6.1.1 Síť v Flow123D ... 46
6.1.2 Síť ve FEFLOW ... 47
6.2 Propustnost ... 49
6.2.1 Propustnost Flow123D ... 49
6.2.2 Propustnost FEFLOW... 49
6.3 Okrajové podmínky ... 50
6.3.1 Okrajové podmínky Flow123D ... 50
6.3.2 Okrajové podmínky FEFLOW ... 50
6.5 Vlastní simulace úlohy ... 51
6.5.1 Simulace ve Flow123D ... 51
6.5.2 Simulace ve FEFLOW ... 533
6.6 Srovnání ... 54
7 Úlohy ... 55
7.1 Úloha 1 ... 55
7.2 Úloha 2 ... 56
7.3 Úloha 3 ... 57
7.3.1 Úloha 3 – zadání OKP ... 58
7.3.2 Úloha 3 – jiné OKP ... 59
10 Závěr ... 60 Literatura ... 61
11
Seznam použitých symbolů
gravitační zrychlení k propustnost horniny
K propustnost
tlak
tlaková výška
piezometrická výška
hustota zdrojů rychlost proudění
polohová souřadnice výšky
hustota
dynamická viskozita
oblast
hranice oblasti
pórovitost
objem pórů
celkový objem zeminy
filtrační rychlost koeficient filtrace
12
Seznam použitých zkratek
OKP okrajové podmínky
MKP metoda konečných prvků
13
Seznam obrázků
Obr.1.: Rozdělení vody, převzato z [4] ... 17
Obr. 2: Freatická zvodeň, převzato z [5] ... 18
Obr. 3: Artézská zvodeň, převzato z [5] ... 19
Obr.4.: Pórovitost různě uspořádaných vrstev kulovitých částic, převzato z [4] ... 20
Obr.5.: Příklady hodnot koeficientu filtrace pro různé zeminy, převzato z [6] ... 21
Obr. 6: Piezometrická výška hladiny ... 22
Obr.7.: Závislost rychlosti proudění na gradientu potenciálu [4] ... 25
Obr.8: Idealizace puklinového systému, převzato z [8] ... 27
Obr. 9: Schematické znázornění Flow123D ... 30
Obr. 10: Prostředí FEFLOW - vlevo Classic, vpravo Standart ... 35
Obr.11: FEFLOW-generování sítě ... 36
Obr.12: FEFLOW - nastavení OKP ... 37
Obr.13: FEFLOW - zobrazení výsledků ... 38
Obr.14: Struktura VTK souboru ... 40
Obr.15.: Program na převod z FEFLOW do vtk ... 41
Obr.16: Nastavení zobrazení proudění ... 42
Obr.17: Zobrazení proudění ve ParaView ... 42
Obr.18: Nastavení ParaView pro srovnání 1 ... 43
Obr.19: Výsledek odečtení ploch od sebe ... 44
Obr.20.: Porovnání hodnot v řezu ... 44
Obr.22: Vygenerovaná síť v GMSH ... 47
Obr.23: Nadefinovaný prostor ... 48
Obr.24: FEFLOW-vygenerovaná síť ... 48
Obr.25.: FEFLOW – zvýraznění pukliny v oblasti ... 49
Obr.26: FEFLOW - nastavení OKP ... 51
Obr.27: FEFLOW - zobrazení výsledku ... 53
Obr.28: Výsledný rozdíl tlaků u sítí s 50 e1ementy ... 54
Obr.29: Výsledný rozdíl tlaků u sítí s 300 elementy ... 54
Obr.30: Zadání první úlohy ... 55
Obr.31: Porovnání první úlohy ... 56
Obr.32: Výsedek druhé úlohy ... 57
14
Obr.33: Zadání třetí úlohy ... 57
Obr.34: Výsledek třetí úlohy s velkou vodivostí puklin ... 58
Obr.35: Výsledek třetí úlohy s malou vodivostí puklin ... 58
Obr.36: Výsledek třetí úlohy s velkou vodivostí puklin OKP na puklině ... 59
Obr.37: Výsledek třetí úlohy s malou vodivostí puklin OKP na puklině ... 59
15
Úvod
V dnešní moderní době se vyvíjejí stále výkonnější počítače pro matematicko- fyzikální výpočty. Jejich důležitou aplikací jsou simulace, které se ve velkém využívají například v průmyslu či stavebnictví. Výhodou simulací je, že dokážeme vytvořit model, který je sice zjednodušením našeho skutečného systému, ale má podobné vlastnosti. A v něm si můžeme vyzkoušet, jak se zachová při určitých okolnostech.
V našem případě se budeme zabývat simulací proudění podzemních vod v rozpukaném prostředí.
Začátek diplomové práce je věnován základnímu seznámení s hydrogeologií a fyzikálními popisy proudění podzemní vody. Následuje postup modelování, ve kterém dochází ke zjednodušení reálného problému v důsledku použití numerických metod při jeho řešení. Jednou z nejpoužívanějších metod numerických metod pro diskretizaci původního spojitého problému a následné výpočty hodnot zkoumané veličiny je metoda konečných prvků. Při diskretizaci modelu záleží, jak moc potřebuje podrobný model.
Protože simulací velkých oblastí nebo použití velmi jemné sítě nám narůstá počet dat, která zabírají více výpočetního času, proto se neustále klade požadavek na zvyšování výkonnosti hardwaru, na kterém jsou výpočty a vizualizace spuštěny. Důležité je i vyvíjet software, který bude stabilně pracovat s takto rozsáhlými daty.
Software je možno rozdělit do několika základních skupin, i když tyto skupiny se mohou různě prolínat. Komerční software, akademický software a freeware.
Akademické softwary mají velkou výhodu, že jsou často zdarma ve formě opensource balíků a jsou zcela zdarma. Opensource, že jsou volně dostupné i zdrojové kódy softwaru. Tudíž je možné software dle potřeb vyvíjet nebo upravovat. Jak již z názvu plyne, tak velké oblibě se tento druh softwaru těší v akademické vrstvě. V diplomové práci používám akademický software Flow123D, komerční FEFLOW, který je použit v demo verzi, kde je maximální počet elementů sítě 500 kusů a opensource program pro vizualizaci Paraview.
Cílem diplomové práce je srovnávání totožné úlohy ve dvou různých simulačních systémech. Nejprve je potřeba námi zadat požadovanou úlohu tak, aby se co nejméně lišila v obou programech. Protože každý používá jinou strukturu sítě, je vhodné si převést výsledky z jednoho programu na síť z druhého. Pro kontrolu,
16 že vše je dobře zadané a nastavené, je vhodné provést si několik jednodušších simulačních úloh, na kterých si tyto vlastnosti ověříme.
1 Základní pojmy hydrogeologie
Hydrogeologie je vědní obor, který se zabývá studiem podzemních vod.
Hydrogeologie zkoumá původ vody, výskyt, zákony pohybu, režim vod, fyzikální a chemické vlastnosti, jejich vzájemné působení s horninami, jejich vztah k atmosférické a povrchové vodě a její hospodářský význam. Historie hydrogeologie je velmi dlouhá a je spojována se zkušenostmi pří získávání vody po dobu několika tisíců let. Lidstvo své nabyté znalosti využívalo převážně k zavlažování půdy a k napájení zvířat. Kromě vyhledávání zdrojů řeší hydrogeologie také záležitosti okolo znečištění a ochrany podzemních vod. V dnešní době hydrogeologové vyhledávají nové zdroje podzemních vod včetně stanovení kvality zdrojů a starají se o ochranu stávajících zdrojů podzemní vody, případně řeší sanaci znečištěných podzemních vod. Dále provádějí hydrogeologický průzkum. Ten je nutný v případech, kde je nutné zjistit hloubku podzemní vody, propustnost hornin či směr a rychlost proudění podzemní vody. Hydrogeologii lze rozdělit do různých odvětví podle specializace jako je například všeobecné hydrogeologie, kontaminační hydrogeologie, ložisková hydrogeologie atd.
1.1 Podzemní vody
Všechna voda, která se nachází pod povrchem země, se většinou označuje jako podpovrchová voda. Dále je to voda, která dopadá na povrch země, částečně prostupuje vlivem gravitace směrem dolů do podloží a hromadí se nad nepropustným podložím.
Podzemní voda proudí prostředím, které je tvořeno pevnou fází a volnými prostory.
Jako volné prostory můžeme uvažovat póry, pukliny apod. Podle toho jak jsou vyplněny relativně skuliny a póry, můžeme rozdělit vodu podpovrchovou na několik horizontálních zón. (viz. Obr.1).
17
Obr.1.: Rozdělení vody, převzato z [4]
Podle toho jak jsou póry zaplněny vodou, můžeme zónu rozdělit na nasycenou a nenasycenou. Nasycené prostředí je, když jsou všechny póry zcela zaplněny vodou, Nenasycenou znamená, že obsahují póry provzdušnění.
Nasycená zóna zahrnuje tři pásma. První je pásmo půdní vody, druhé přechodové pásmo a třetí částečné pásmo kapilární vody.
Pohyb kapilární vody je ovlivňován převážně kapilárními silami, což jsou síly v kapilárách na fázovém rozhraní mezi plynnou, kapalnou a pevnou fází. Při plném nasycení zanikají.
Podzemní voda je ze shora ohraničena volnou hladinou podzemí vody a ze zdola nepropustným podložím. Na volné hladině je hydraulický tlak roven tlaku atmosférickému.
Z hlediska energetického rozdělení půdní vody se dělí na absorpční, kapilární a gravitační. Voda absorpční zahrnuje molekuly vody poutané k povrchu pevných částic adsorpčními a osmotickými silami. Nachází se v horninách při nízké vlhkosti, kde je v kapalném stavu zcela nepohyblivá. Dále se pak dělí na kapilární vodu, u které je pohyb závislý na kapilárních silách. Jsou to síly vznikající při vytváření menisku spojeného s adhezními silami mezi vodou a povrchem pevných částic a kohezních sil mezi molekulami vody. Poslední voda je gravitační voda, její pohyb je určován gravitací. Rozděluje se podle původu. Voda juvenilní, která je magmatického původu (pochází z hloubky), voda meteorická pocházející z dešťových srážek a dostávající se do země vsakováním. Konátní je voda, v níž byl sediment uložen a ona se v něm v části může zachovat. Při povrchu leží pásmo provzdušnění, kde část pórů je vyplněna vzduchem. Voda pohybující se v této části se nazývá voda vadózní. Níže leží pásmo
18 nasycení, kde jsou všechny póry zaplněny vodou. Svrchní hranicí tohoto pásma je hladina podzemní vody a v podstatě odpovídá horní hranici zvodně.
1.2 Zvodeň
Zvodeň je akumulace gravitační podzemní vody a to buď volné, nebo napjaté.
Volná voda (freatická voda) je v malé hloubce pod povrchem, snadno dostupná k čerpání (Obr. 2) a pohybuje se především díky vlivem gravitačních sil.
Obr. 2: Freatická zvodeň, převzato z [5]
Artézská zvodeň se vyznačuje hydrostatickým tlakem větším než tlak atmosférický vlivem nepropustného nadloží (viz Obr.3). Hornina, v níž se voda hromadí, se nazývá kolektor a jeho fyzikální vlastnosti určují, kolik vody je schopna pojmout. Hlavními vlastnostmi jsou pórovitost a propustnost, tedy množství pórů v hornině, jejich velikost a vzájemné spojení o rozpukání hornin.
19
Obr. 3: Artézská zvodeň, převzato z [5]
1.3 Porozita
Porozita udává přítomnost jakýkoliv pórů a dutin v horninách. Vyznačuje se jako podíl objemu pórů na celkovém objemu horniny. Takto spočítaná pórovitost se nazývá koeficient pórovitosti.
Vztah pro výpočet pórovitosti
(1)
Kde je pórovitost je objem pórů je celkový objem zeminy
Pórovitost může ovlivňovat například i vzájemné uložení zrn. Na Obr.4 je znázorněno několik typických příkladů pravidelného uspořádání kulových částic ve vrstvě. U volně nasypaných vrstev složených z běžných zrnitých materiálů hornin se pórovitost pohybuje v rozmezí 0,05 – 0,25.
20
Obr.4.: Pórovitost různě uspořádaných vrstev kulovitých částic, převzato z [4]
Efektivní pórovitost, též zvaná účinná pórovitost, vyjadřuje objem gravitační vody, který vyteče z plně nasyceného vzorku zeminy. Se zvyšující se velikostí zrna zeminy pórovitost obvykle klesá a efektivní pórovitost stoupá (obr. 4). Efektivní pórovitost například u štěrku je 25%, písku 20%, pískovce 10% a jílu 3%.
Homogenní porézní prostředí je dané ke své vlastnosti (např. hydraulické vodivosti), pokud je ve všech bodech tato jeho vlastnost stejná. Pokud se vlastnost mění v závislosti na poloze ve vyznačené oblasti, tak se pak jedná o nehomogenní prostředí.
1.4 Propustnost
Propustnost je schopnost horniny propouštět kapaliny (vodu, ropu) nebo plyn v horninovém pórovitém prostředí (pískovec, štěrk) nebo propouštět médium podél puklin apod., její hodnota může být různá v každém směru. Propustnost je vztažena k proudění kapaliny o určitých vlastnostech (podzemní voda), vyjadřuje se pomocí koeficientu filtrace k (viz. Obr.5). V přírodě neexistuje hornina, která by byla absolutně nepropustná. Některé horniny však mají tak malou propustnost, že se považují za nepropustné (např. mastné jíly). Především je ovlivňována typem horniny. Jednotka propustnosti se nazývá Darcy (budeme se jí ještě dále věnovat). Nejlépe propustné jsou hrubé sedimenty bez jílovité příměsi (písky, sutě, vápence atd.)
21
Obr.5.: Příklady hodnot koeficientu filtrace pro různé zeminy, převzato z [6]
2 Proudění podzemní vody
Pohyb tekutin, při kterém se částice tekutiny pohybují neuspořádaným pohybem a zároveň se posunují v určitém směru, se nazývá proudění. Než se dostaneme k rovnicím, které nám budou popisovat proudění, popíšeme si nejprve použité veličiny.
V hydraulice podzemní vody se pracuje s piezometrickou výškou značí to nadmořskou výšku hladiny vody ve vrtu. Ta se dá spočítat jako součet nadmořské výšky bodu měření a tlakové výšky v tomto bodě (viz. Obr. 6).
(2)
kde (3)
je gravitační zrychlení a - hustota vody a výsledný vzorec je:
(4)
22
Obr. 6: Piezometrická výška hladiny
Další veličinou je filtrační rychlost . Filtrační rychlost je definována jako objem podzemní vody, který proteče jednotkou plochy.
(5)
Filtrační rychlost ovlivňují vlastnosti filtrace jako je vlastnost vody (viskozita) a tlakový spád (stlačitelnost).
Koeficient filtrace se číselně rovná filtrační rychlosti při jednotkovém hydraulickém gradientu I. Má rozměr rychlosti a představuje směrnici lineární závislosti proudění na hydraulickém gradientu v Darcyho zákoně. Q je objemový průtok a A je průřez oblasti.
Nejprve proudění vody v podzemí budeme předpokládat jako ustálené proudění s napjatou hladinou. Tento děj se dá popsat rovnicí kontinuity ve tvaru:
na (6)
kde q v je to hustota zdrojů, je filtrační rychlost a je oblast, ve které zkoumáme proudění.
Dále budeme potřebovat Darcyho vztah, který popisuje proudění kapaliny v porézním prostředí. Jde o závislost mezi filtrační rychlostí tekutiny a hydraulickým gradientem.
na (7)
23 rychlost proudění se udává v , protože často jsou průtoky velmi malé,
používá se většinou jednotek , . K představuje tenzor
propustnosti.
na (8)
Většinou používáme propustnost ve tvaru , kde je číslo a je jednotková matice. Propustnost K obsahuje v sobě jak vlastnosti horniny, tak i vlastnosti proudící kapaliny, proto je vhodné ji používat jen v případech, kde se vlastnosti tekutiny výrazně nemění, což je splněno ve většině úloh proudění podzemních vod. V případě, kdy jsou vlastnosti kapaliny výrazně závislé na tlaku a teplotě, jako např. při proudění uhlovodíků (ropa), je výhodné používat propustnost závislou pouze na vlastnostech horniny k. Mezi K a k je vztah:
(9)
Kde pro k se používá jednotka Darcy , fyzikální rozměr má , je dynamická viskozita kapaliny, , hustota tekutiny , gravitační zrychlení
.
Jednotka Darcy je definována jako propustnost pro tekutinu s viskozitou 1 cP (centipoisson), kdy objemový tok 1 laminárně profiltruje při gradientu 1 atm
profilem o ploše 1 . Platí, že , jednotka je to příliš
velká, proto se v praxi používá .
Propustnost K se zadává v , což je pro nás někdy příliš velká jednotka, je proto vhodnější si ji zadávat v jednotkách . Proto je hodné přepočítat si vztah mezi k a K:
24
dosazení do rovnice:
po přepočtu nám vyjde:
2.2 Omezení Darcyho zákona
Darcyho zákon, který jsme si popsali v kapitole 2, je lineární zákon, který vyjadřuje lineární závislost filtrační rychlosti na hydraulickém gradientu. Tato lineární závislost neplatí na celém rozmezí hodnot, ale je omezena dolní a horní limitní hodnotou gradientu (viz Obr. 7).
25
Obr.7.: Závislost rychlosti proudění na gradientu potenciálu [4]
Při průsaku velmi jemnozrnným materiálem s nízkou propustností existuje limitní hodnota hydraulického gradientu, při které ustává pohyb kapaliny. Druhé omezení Darcyho zákona je při průsaku velmi hrubozrnným materiálem, při kterém dochází k nelineární závislosti mezi růstem gradientu potenciálu a růstem rychlosti. Horní limit může být překročen při proudění v krasových vápencích a dolomitech a ve vulkanických horninách s kavernami. Pro běžné případy proudění podzemní vody je lineární Darcyho zákon aplikovatelný.
2.3 Okrajové podmínky
Na hranici vymezené oblasti , označované jako , musíme definovat okrajové podmínky (OKP). Existují tři typy podmínek. První je Dirichletova okrajová podmínka, která zadává piezometrickou výšku
na (10)
26 Druhá je Neumannova okrajová podmínka udávající množství kapaliny, která protekla přes hranici oblasti
na (11)
Často se používají tvz. homogenní Neumannovy okrajové podmínky, kde je přetok roven nule. To v praxi znamená, že je hranice okolí nepropustná.
Newtonova okrajová podmínka je zadaná lineární kombinací Dirichletovy a Neumanovy okrajové podmínky.
na (12)
kde vyjadřuje vlastnosti hranice.
2.4 Proudění v rozpukaném prostředí
Protože je velice složité a prakticky nemožné vystihnout přesnou strukturu prostředí ve všech detailech hrající roli v procesech proudění (od mikrotrhlin až k velkým zlomům), je potřeba simulované prostředí do jisté míry zjednodušit a nahradit modelem vystihujícím podstatné věci pro dané účely a v daném měřítku.
V praxi se uplatňuje několik koncepcí.
Náhrada ekvivalentním porézním médiem (označována jako EPM, Equivalent Porous Medium [7]) představuje „průměrováním“ vlastností reálného horninového prostředí v modelu proudění. To znamená, že všechny individuální hydraulické charakteristiky vyskytující se v systému propojených puklin jsou nahrazeny efektivní hydraulickou vodivostí. Výhodou této metody je zpravidla to, že jeho výpočetní realizovatelnost je méně náročná a dají se díky tomu vytvářet rozsáhlé oblasti.
Vlastnosti propojených propustných pukliny systému slepých póru jsou průměrovány v celém simulovaném bloku. Náhrada ekvivalentním porézním médiem se používá zpravidla pro simulaci globálního charakteru toku v dané oblasti.
27 Dual porozity a Dual permeability je pokročilá funkce, která převede požadovanou oblast na matici bloků (viz Obr. 8) mezi nimiž jsou pukliny. Koncept dvojí porozity a permeability je založen na idealizovaném modelu, ve kterém existují dvě překrývající kontinua, která se ovlivňují. V každém bodě je známa permeabilita a porozita .
V jednom kontinuu platí
(13)
(14)
v druhém
(15)
(16)
a mezi těmito kontinui probíhá výměna pomocí členu výměny.
Obr.8: Idealizace puklinového systému, převzato z [8]
Další metodou je koncepce založená na metodě diskrétních puklinových sítí (zkratkou DFN, Discrete Fracture Network [7]) kde je simulované prostředí nahrazeno sítí puklin. Proudění potom probíhá pouze po síti propojených puklin. Hornina mezi puklinami
se v modelu může i nemusí zanedbávat. Potřebná puklinová síť je zpravidla generována stochasticky na základě známé pravděpodobnostní distribuci puklin v reálném horninovém prostředí. Tyto modely zpravidla umožňují zahrnout hydrogeologicky významné pukliny. Hydraulické vlastnosti jsou v modelu průměrovány na úrovni
28 jednotlivých puklin. Problémy tohoto přístupu se zpravidla vztahují k výpočetní náročnosti výsledného modelu. Počet puklin významně narůstá s velikostí simulované oblasti. Koncepce diskrétních puklinových sítí se používá pro simulace proudění a transportu v menším měřítku.
3 Metoda konečných prvků
Metoda konečných prvků se používá při řešení fyzikálních problémů v analýze návrhu nových zařízení. Funguje na principu, že oblast rozdělíme na drobné elementy (síť) a převedeme řešení parciálních diferenciálních rovnic na diskrétní problém (tj. řešení soustavy lineárních rovnic). Při modelování fyzikálního problému musíme brát v úvahu jisté předpoklady, které výsledný model sice znepřesní, ale umožní nám ho vyřešit. Po zidealizování problému dostaneme matematický model, který se skládá z diferenciálních rovnic. Parciální diferenciální rovnice jsou rovnice, v nichž hledáme funkci u na základě daného vztahu mezi jejími parciálními derivacemi. Řešení takové úlohy nebude stoprocentně odpovídat realitě, protože parciální diferenciální rovnice většinou řešíme uvnitř nějaké otevřené množiny, přičemž na hranici této množiny zpravidla uvažujeme dodatečné tzv. okrajové podmínky. Dále nedokážeme popsat plně přírodní děj a také při řešení matematického modelu dochází k jistým chybám.
Výsledek můžeme zpřesnit více způsoby. Zaprvé upravení přímo samotného řešení metody konečných prvků a to tak, že zjemníme síť. Další možností je upravit matematický model. Můžeme změnit geometrii, kinematiku, okrajové podmínky, vlastnosti materiálů.
Hlavním krokem je volba matematického modelu, který vybíráme na základě zkoumaných vlastností objektu. Potřebujeme-li získat co nejpřesnější a nejefektivnější model…
Ve světě se stále častěji uplatňuje přístup, kdy bloky hornin reprezentované porézním médiem obsahují deterministické pukliny reprezentující hydraulicky propustné pukliny, či rozsáhlé tektonické zóny.
29
4. Simulační programy
Simulačních programů pro proudění kapalin v prostředí je velmi mnoho.
Ať už komerční nebo volně dostupné. V práci jsme použili programy Flow123D, který je vyvíjený na Technické univerzitě v Liberci a komerční FEFLOW. Podrobněji si popíšeme tyto programy v dalších kapitolách.
4.1 Flow123D
V případě simulačního programu Flow123D, který je vyvíjen na ústavu NTI v Liberci, byla přijata koncepce umožňující kombinaci porézního modelu s diskrétní puklinovou sítí. Tato koncepce umožňuje individuálně u každé realizace modelu konkrétní oblasti zvážit, do jaké míry bude provedeno zjednodušení simulovaného reálného prostředí na jedné straně, a do jaké míry bude zohledněna síť hydrogeologicky významných puklin nebo regionálních poruch na straně druhé. Software potom umožňuje využít ryze porézní model, popřípadě ryze puklinový přístup. Pro aproximaci je zde použita metoda konečných prvků založená na smíšené hybridní formulaci úlohy proudění. Software umožňuje výpočet ustáleného i neustáleného proudění a je doplněn o transport rozpuštěných látek zahrnující advekci, sorpci a difuzi do imobilní zóny.
Program Flow123D se nyní nachází ve verzi 1.6, které podporuje už i export do různých výstupních souborů. Jedná se o konzolovou aplikaci, která se spouští výhradně z příkazového řádku příkazem:
flow123 [parametr] ini_file
význam parametrů s – spočítá úlohu proudění, příp. transportu c – konverze vstupních souborů pro vizualizaci n – generování souboru sousednosti NGH (v novějších verzí programu se generují NGH samostatným programem
30 ini_file je konfigurační soubor, ve kterém si nastavujeme vstupní soubory (síť, okrajové podmínky, materiály a další), průběh výpočtu a výstupní soubory.
Znázornění všech potřebných souborů je na Obr. 9. Výstupní soubor *.log obsahuje vypsaný celý průběh výpočtů a druhý soubor *.pos obsahuje výsledky výpočtů ve formátu pro zobrazení programem GMSH. Nebo místo souboru *.pos můžeme provést výpočet do souboru*.vtu, který lze otevřít v jiných programech jako ParaView, což se nám velice hodí. Protože program ParaView kromě zobrazení umí i různé filtry a porovnávání více výsledků z rozdílných programů, které použijeme.
4.1.1. Struktury souborů pro Flow123D
Nyní si popíšeme obecně konkrétní typy souborů potřebné pro výpočet a jejich strukturu.
4.1.1.1 Síť *.msh
První soubor, který si popíšeme je *.msh. Obsahuje celou síť, jeho struktura je:
$MeshFormat
2.0 file-type data-size
$EndMeshFormat
$Nodes
number-of-nodes
node-number x-coord y-coord z-coord ...
$EndNodes
Flow123D
*.msh *.mtr *.bcd *.ngh *.src *.ini
*.pos *.vtu *.log
Obr. 9: Schematické znázornění Flow123D
31 Number-of-nodes obsahuje celkový počet uzlů. Další řádky obsahují číslo uzlu a jeho souřadnice v x.y a z. Další část obsahuje už vytvořené elementy. Opět jsou očíslovány od jedničky (po vygenerování GMSH, ale není to podmínka) a na začátku je uveden celkový počet elementů. Ve druhém sloupci jsou definovány geometrické typy elementů (1-linka definovaná 2 uzly, 2 – trojúhelník definovaný 3 uzly atd. více v manuálu k Flow123D [3]). Dalším důležitým sloupcem je čtvrtý sloupec. Ten udává označení materiálu. Podrobnější vlastnosti materiálů si nadefinujeme v souboru s materiály.
Poslední sloupce se věnují místu, kde bude element zobrazen. Udávají čísla uzlů, ve kterých tvoří daný geometrický element.
$Elements
number-of-elements
elm-number elm-type number-of-tags <tags> node-number-list ...
$EndElements
Pro ukázku jednoduchý příklad:
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes -hlavička začátku uzlů
4 -celkový počet uzlů
1 0 0 0 -souřadnice uzlů (x=0, y=0, z=0)
2 10 0 0
3 0 5 0
4 10 5 0
$EndNodes -konec uzlů
$Elements -hlavička začátku elementů
2 -celkový počet elementů
1 1 2 1008 7 0 1 2
- 1D element, materiál 1008 a je mezi uzly 1 a 2
2 2 2 1009 7 1 2 3
- 2D element, materiál 1009 a je mezi uzly 1, 2 a 3
$EndElements -konec elementů
32
4.1.1.2 Materiály *.mtr
Mtr obsahuje vlastnosti materiálů, které budeme potřebovat. Dá se zde definovat typ materiálu, geometrie, sorpční vlastnosti atd.
$MaterialFormat 1.0 file-type data-size
$EndMaterialFormat
$Materials
number-of-materials
material-number material-type <material-type-speci_c-data> [text]
$EndMaterials
$Geometry
material-number geometry-type <geometry-type-specific-coeficient> [text]
$EndGeometry
$Sorption
material-number substance-id sorption-type <sorption-type-specific-data> [text]
$EndSorption
Nás bude hlavně zajímat Materials. Nejprve určíme celkový počet materiálů.
Pak si je označíme material-number, například 1008 (na který jsme se odkazovali při výrobě sítě). Pak material-type, typ materiálu, například 21 znamená, že , ostatní typy jsou uvedeny v manuálu ke Flow123D [3]. Poslední sloupec udává propustnost.
4.1.1.3 Okrajové podmínky *.bcd
Okrajové podmínky jsou opět definovány hlavičkou a pak daty OKP podmínek.
$BoundaryFormat 1.0 file-type data-size
$EndBoundaryFormat
$BoundaryConditions
condition-number type <type-specific-data> where <where-data>
number-of-tags <tags> [text]
33 type – typ okrajové podmínky : 1 - podmínka typu Dirichlet
2 - podmínka typu Neumann 3 - podmínka typu Newton
<type-specific-data> - zadání hodnot podle vybrané podmínky, je-li vybrána 1 – Dirichlet, hodnoty značí tlakovou výšku
<where> - elementy na kterých se nachází tato OKP
<tags> - definovaná podmínka, 1 = typ voda
Abychom nemuseli jednotlivé body vypisovat ručně, používá se program BCD a z něho je výstupní soubor*.fbc. K nastavení programu se používá ini_soubou se strukturou:
$Input
Mesh_file = ./input/test.msh - cesta k síti
Ngh_file = ./input/test.ngh - cesta k sousednostem N_subst = 1
$Output
Output_digits = 12 Accuracy = 1E-16 Write_all_BC = NO
FBC_file = ./input/test.fbc - kam se uloží výsledné OKP
$FBC
2 1 0 0 0 1 0 0 -1 100 1 1
– definování OKP na celé levé straně, pomocí Dirichleta a tlakové výšky 100m
2 1 0 0 -10 1 0 0 -1 120 1 2
– definování OKP na celé pravé straně vzdálené 10metrů od levé, pomocí Dirichleta a tlakové výšky 120m
$EndFBC
Podrobněji popsáno v kapitole 6.3.1.
34
4.1.1.4 Sousednosti *.ngh
Soubor obsahující sousednosti je přímo generován programem nacházející se u Flow123D. Výsledná struktura souboru je:
$NeighbourFormat 1.0 file-type data-size
$EndNeighbourFormat
$Neighbours
number-of-neighbours
neighbour-number type <type-specific-data>
$EndNeighbours
Udává počet sousedností a typy sousedností (hranice se společnými uzly, hranice se specifickými uzly atd.)
4.1.1.5 Zdroje *.src
Zde se definují zdroje. Která vyznačují oblast, kde se čerpá nebo vtláčí kapalina.
A kolik kapaliny skrz oblast prochází.
$SourceFormat 1.0 _le-type data-size
$EndSourceFormat
$Sources
number-of-sources eid density
...
$EndSource
eid – číslo elementu, na kterém se nachází
destiny – hustota zdroje v , kladné hodnoty jsou pro vtláčení, záporné hodnoty pro čerpání.
35
4.2 FEFLOW
FEFLOW se řadí mezi profesionální software pro modelování proudění a transport rozpuštěných látek nebo tepla v podpovrchových procesech. Tento program byl vyvinut v Německé akademii věd v Berlíně na Ústavu mechaniky. Nyní je jako komerční software vyvíjený skupinou DHI, což je skupina zabývající se technologickým rozvojem v oblasti vodního hospodářství a životního prostředí.
Program pracuje na základech metody konečných prvků. Jeho hlavní výhodou je, že se jedná o integrovaný systém, který od návrhu simulace v grafickém rozhraní umožňuje i prohlížení výsledků. Program je používán výzkumnými ústavy, vysokými školami, firmami a vládními organizaci po celém světě. Dnes obsahuje už dvě verze a to FEFLOW Classic, která obsahuje starší uživatelské prostředí pocházející ze starších verzí pro Linux. A novou verzi FEFLOW Standart, která má přepracované ovládání (viz. Obr.10). Testovací úlohy se vytvářely v prostředí Classic, kvůli přehlednější obsluze.
Obr. 10: Prostředí FEFLOW - vlevo Classic, vpravo Standart
36
4.2.1 Modelování ve FEFLOW 4.2.1.1 Síť ve FEFLOW
V sekci Edit- supercelement mesh si navrhneme náš 2D prostor. Protože je síť generována automaticky, nemůžeme ovlivnit, kde a jak se zobrazí elementy. Například pokud máme použity dvě různé látky (viz Obr.11) nebo pukliny. Je nutné už při návrhu polygonů počítat s jejich rozmístěním. Pak přejdeme ke generování sítě. Zde se nastavuje použitý typ elementů (např. trojúhelníky, čtverce). Pokud stačí defaultní nastavení, tak pro vytvoření sítě potřebujeme pouze jeden parametr, a to počet elementů. Tím ovlivňujeme jemnost sítě.
Obr.11: FEFLOW-generování sítě
4.2.1.2 Propustnost FEFLOW
Nastavení vlastností materiálů se nachází v sekci Edit v Problem attributes.
Zde se nastavuje kromě zmíněných věcí i problém, který máme simulovat. V našem případě (se jedná o) podzemní proudění bez transportu látek. Poslední věcí, kterou zde můžeme nastavit, jsou vlastnosti materiálů. (Flow materials). Zde si nastavíme propustnosti materiálů, buď pro celou plochu najednou, nebo každému elementu zvlášť.
37
4.2.1.3 Okrajové podmínky Feflow
Nastavení okrajových podmínek se provádí ve Flow data. Zde se nacházejí okrajové podmínky (boundaries). Podmínky jsou vypsány jak podle typu, jak je uvedeno v kapitole o okrajových podmínkách, tak podle veličiny, kterou zadáváme (např. Head (1st kind). Což je podmínka Dirichletova, při které se zadává tlaková výška hladiny. Po vybrání podmínky a zadání hodnoty si myší vybereme hrany elementů, na kterých má podmínka platit (viz Obr.12).
Obr.12: FEFLOW - nastavení OKP
4.2.1.4 Zdroje FEFLOW
Zdroje u FEFLOW se zadávají u vlastností materiálů. Zde se pouze zadá, zda se voda bude vhánět či čerpat v jednotkách [10-4/d] a označí se požadovaný element.
4.2.1.5 Samotná simulace FEFLOW
Samotná simulace probíhá v Run – start simulation. Zde je důležité si před samotným spuštěním nastavit nahrávání simulace do datového souboru. V Control output vybereme Record komplete data ASCII a místo souboru kam se má uložit. Nyní můžeme pustit samotnou simulaci. FEFLOW rovnou ukazuje průběh simulace. Můžeme si dokonce zvolit i body, ve kterých nám program ukáže průběh veličiny v daném bodě.
Po skončení simulace se nám zobrazí výsledek. Ten pomocí uživatelského prostředí
38 můžeme prohlížet a zkoumat. Pro další porovnávání je vhodné si po simulaci vyexportovat soubor s hodnotami sledovaných veličin v uzlech elementů.
Obr.13: FEFLOW - zobrazení výsledků
5 Způsob srovnávání výsledků
Abychom mohli porovnat výsledky obou programů, bylo by vhodné najít nějaké řešení grafické či početní. Jelikož síť v každém programu má jinou strukturu, jsou výstupní data z nich také odlišná. Největší rozdíl je v umístění porovnávaných bodů.
Tím, že každá síť je jiná, výsledky v uzlech se nachází na jiných souřadnicích a proto je potřeba provést aproximaci, která nám pomůže porovnat dvě úplně odlišné sítě.
Nakonec byl zvolen program ParaView, který umí přepočítávat data z jedné sítě na
druhou a navíc
je umí i porovnat mezi sebou. Ovšem podporu pro program ParaView má pouze Flow123D (verze 1.6). Takže se objevil problém jak dostat data z FEFLOW. Byl zvolen
39 soubor VTK, který má jednoduchou strukturu a do něj se převáděly data z vyexportovaných souborů.
5.1 Struktura VTK
Pro volbu formátu *.vtk byl hlavní důvod pro vytvoření jednotného schéma pro různě popsaná data. Existují dva základní typy souboru. První je jednoduchý, určený pro přehledné zorientování. Druhý je složitější XML formát, který poskytuje více možností konfigurace. Nám postačí základní formát. Ukázka výsledného souboru (viz. Obr.14)
Standartní VTK je rozdělen na pět částí.
1. Typ souboru a verze VTK, která musí být přesně definována pro komunikaci a kompatibilitu se softwarem.
2. Hlavička. Slouží pro naši orientaci, co následující data popisují. Její délka musí být maximálně 256 znaků, typu string.
3. Formát souboru. Máme dvě možnosti. Buď ASCII nebo binární kódování. My použijeme typ ASCII.
4. Popis geometrie objektu. Nejprve je nutné definovat
datový set (strukturovaná, nestrukturovaná síť, atd.). Poté souřadnice bodů, spojované body a tvar, jež spojované body představují.
5. V poslední části popisujeme konkrétní datové hodnoty pro jednotlivé body (např. vektory, skaláry, tensory, atd.)
40
Obr.14: Struktura VTK souboru
41
5.2 Převod z FEFLOW do VTK
Na převod z FEFLOW budou potřeba dva vyexportované soubory, které získáme přímo z programu. Jedná se o soubor *.dac, který se vytvoří během průběhu simulace, z něho náš nejvíce zajímají vrcholy elementů (NODE). Ty se nejprve musí upravit a to tak, že všechna čísla se musí zmenšit o hodnotu 1, protože FEFLOW čísluje od 1, ale v souboru VTK je potřeba číslovat od 0. Když máme takto upravené hodnoty, můžeme je nahrát do souboru *.vtk do části CELLS. Další soubor, který potřebujeme, má koncovku *.dat a obsahuje souřadnice (x,y,z) a výsledné hodnoty (např. tlaku).
Souřadnice nahrajeme do DATASET UNSTRUCTURED_GRID a hodnoty tlaku do POINT_DATA. Chybí nám už jen informace o CELL_TYPES, které si sami vytvoříme. Jedná se o tvar elementů, v našem případě jde o trojúhelník a ten má hodnotu 5 a výskyt tohoto čísla je dán počtem uzlů. Abychom po každé simulaci nemuseli tento proces opakovat, vytvořil jsem jednoduchý program, který tento postup udělá za nás. Je velice jednoduchý na ovládání. Pomocí tlačítek přiřadíme cestu k požadovaným souborům, pak vygenerujeme vtk soubor.
Obr.15.: Program na převod z FEFLOW do vtk
42
5.3 ParaView
Jedná se o open-source vizualizační program, který je vhodný pro prohlížení a úpravu dat. V programu si pomocí filtrů uděláme postup pro srovnání vybraných ploch. Princip spočívá v tom, že si data z jedné sítě převedeme na stejnou síť pomocí interpolace. Poté, kdy máme data na stejné síti, můžeme je od sebe odečíst. Tím uvidíme, jak velké rozdíly se nacházejí na ploše a kde se přesně nacházejí.
5.3.1 Zobrazení v ParaView
Zobrazení v ParaView probíhá tak, že načteme naše nasimulované soubory z Flow123D a FEFLOW (po převodu na vtk) a v Object Inspector klikneme na Apply.
Tím se nám zobrazí naše síť. Pokud chceme zobrazit tlakovou výšku, musíme vybrat vrstvu, která ji obsahuje. Pokud ale chceme zobrazení vektorů proudění, je potřeba aplikovat dva filtry. Nejprve Cell Centers a poté Gyph (viz Obr. 16).
Obr.16: Nastavení zobrazení proudění
Díky těmto filtrům se nám zobrazí proudění v naší oblasti. (viz Obr.17)
Obr.17: Zobrazení proudění ve ParaView
43
5.4 Postup porovnání v ParaView
V této kapitole si ukážeme a popíšeme úkony, které jsou nutné udělat, pro porovnání.
Výsledkem je grafické znázornění rozdílů hodnot na vybrané ploše.
1 Načteme soubory z obou programů pomocí Calculatoru. Tím si vyextrahujeme pouze námi porovnávanou veličinu, se kterou budeme dále pracovat.
2 Teď použijeme filtr ResampleWithDataset. Tím docílíme přepočtu dat z jedné sítě na druhou. Do kolonky input vybereme jemnější síť a do source síť hrubší.
3 Dále použijeme filtr Apend Atributes, který nám vypočtené hodnoty z filtru ResampleWithDataser spojí dohromady. Výsledkem je jedna síť, která obsahuje data z obou programů.
4 Na konec pomocí Calculatoru můžeme různě pracovat s hodnotami v síti. Pro zobrazení rozdílu ploch je vhodné tyto hodnoty od sebe odečíst (viz Obr.19).
Při vytváření filtrů se nám v Pipeline Browseru budou zobrazovat veškeré naše kroky.
Výsledek by měl vypadat přibližně takhle:
Obr.18: Nastavení ParaView pro srovnání 1
44
Obr.19: Výsledek odečtení ploch od sebe
Další věcí, kterou můžeme v ParaView použít pro porovnání, je Plot Over Line, který nám zobrazí hodnotu např. tlaku v řezu. Na obrázku jsme udělali průřez naší oblastí podél osy x.
Obr.20.: Porovnání hodnot v řezu
45
6 Nastavení stejné úlohy v obou systémech
Abychom mohli porovnávat naše nasimulované výsledky, bude potřeba, aby se úloha do obou programů zadala co nejvíce shodná v obou programech. Aplikace Flow123D pro svůj výpočet, jak jsme si již řekli, potřebuje nadefinovat několik vstupních souborů a správně napsaný ini soubor. Navíc pro vytvoření sítě se využívá program GMSH. FEFLOW je komplexní systém, který obsahuje vše pro tvorbu, simulaci a zobrazení výsledků. Proto se liší kromě tvorby sítě i v zadávání dalších veličin, jako jsou okrajové podmínky, vlastnosti materiálů či definování vrtů. Proto se v této kapitole budeme věnovat tomu, jak nejlépe tyto parametry zadat tak, aby si byly úlohy co nejvíce podobné. A aby vše bylo vidět, ukážeme si to na konkrétním jednoduchém příkladu s puklinami. (konkrétní popis struktury souborů u Flow123D je popsán v kapitolách o Flow123D).
6.1 Zadání příkladu
Je zadán kvádr 10x5m. Na levé a pravé straně je zadaná Dirichletova podmínka.
Vlevo je tlaková výška 100m a vpravo 120m. Pak jsou zde dvě pukliny. Spodní puklina má 100x větší propustnost než okolní prostředí a vrchní má 100x menší propustnost.
10m 4m
3m 5m
5 3m m
2 m
1 m
2 m
Obr.21.: Zadání úlohy
46
6.1 Síť
Nejprve si vytvoříme síť pro naši oblast. Pro tento příklad si připravíme dvě sítě, jednu hrubší o 50 elementech (pro tu budeme popisovat nastavení), druhou jemnější o 300 elementech.
6.1.1 Síť v Flow123D
Síť jsme vytvářeli pomocí programu GMSH z důvodu, že chceme strukturovanou síť, kterou program vyžaduje pro svůj výpočet. Můžeme si ji graficky vytvořit nebo lépe naprogramovat pomocí souboru *.geo. V něm použijeme funkci extrude (vytažení), která nám hodně usnadní práci. Nejprve si vytvoříme bod, který pomocí této funkce vytáhneme v daném směru (například ve směru osy x). Se vzniklou čárou provedeme obdobný proces v jiné ose a nakonec to samé s plochou. Zde je ukázán v geo programu pro naši úlohu)
Extrude Point { 1, { 10.0, 0.0, 0.0 }} { Layers { 10, 1.0 };
Recombine;
};
Extrude Line {21, { 0.0, 5, 0.0 }} { Layers { 5, 1.0 };
};
Extrude Surface { 3, { 0.0, 0.0, 0.0 }} { Layers { 1, 1.0 };
};
Při otevření souboru *.geo v textovém editoru vidíme vypsané příkazy a snadno si zde můžeme měnit dělení sítě nebo její celkové rozměry. Po vytvoření tohoto objektu stačí přepnout v GMSH na Mesh a vybrat vytvoření 2D sítě a program si ji sám vygeneruje (Obr.22).
47
Obr.22: Vygenerovaná síť v GMSH
Pokud máme vytvořenou síť, musíme do ní dodefinovat pukliny. Editujeme soubor
*.msh a do sekce Elements vytvoříme nové 2D, které nám budou popisovat pukliny. Pro příklad ukázka zapsání jedné pukliny mezi dvěma uzly 12 a 133 a typ materiálu 1008.
Původní 2D elementy mají 1009 .
1 1 3 1008 7 0 12 133 - 1D elementy
2 2 3 1009 7 2 12 133 - 2D elementy
6.1.2 Síť ve FEFLOW
Navrhneme si náš prostor, viz kapitola 4.2.1.1. Nesmíme zapomenou při návrhu zakreslit pukliny (viz Obr.23). Pak přejdeme ke generování sítě, kde si zvolíme počet elementů sítě.
48
Obr.23: Nadefinovaný prostor
Obr.24: FEFLOW-vygenerovaná síť
49
6.2 Propustnost
6.2.1 Propustnost Flow123D
Zadávání propustnosti K se nastavuje u materiálů v souboru *.mtr, kde definuje vlastnost materiálů a propustnost. Zde si nadefinujeme i propustnost puklin. Výsledek je:
$Materials 3
1007 11 0.01 - 100x propustnější puklina 1008 11 100 - 100x propustnější puklina 1009 21 1 - propustnost materiálu v naší oblasti
6.2.2 Propustnost FEFLOW
Nastavení vlastností materiálů se nachází v sekci Edit v Problem attributes. Zde si nastavíme propustnosti horniny. Pro nastavení propustností puklin, využijeme sekci pro diskrétní elementy (Discrete feature elelements). Tam zadáme nejprve hodnotu propustnosti a potom označíme požadované 1D elementy.
Obr.25.: FEFLOW – zvýraznění pukliny v oblasti
50
6.3 Okrajové podmínky
6.3.1 Okrajové podmínky Flow123D
Okrajové podmínky se mohou nastavit přímo editací souboru *.bcd, který má strukturu popsanou v kapitole 4.1.1.3 Okrajové podmínky. Nevýhoda tohoto způsobu je, že pro větší nebo složitější modely je velmi pracné vypisovat všechny stěny s okrajovými podmínkami, proto byl použit program BCD.
Program BCD je opět konzolovou aplikací, využívající konfigurační ini_soubor, kde se vkládá název souborů sítě a sousedností, přesnost výpočtu a název výstupního souboru s příponou *.fbc. V prvním sloupci je číslo podmínky, ve druhém sloupci je plocha, na které má být podmínka definována ve tvaru ax+by+cz+d=0, ve třetím sloupci je opět definovaný typ okrajové podmínky 1.Dirrichlet, 2.Neumann, 3.Newton. Ve čtvrtém sloupci je hodnota tlakové výšky a v posledním se vybírá číslo stěny elementu, ke které bude hodnota přiřazena.
$FBC
2 0 1 0 0 1 0 0 -1 100 1 1
2 0 1 0 -10 1 0 0 -1 120 1 2
$EndFBC
6.3.2 Okrajové podmínky FEFLOW
Nastavení okrajových podmínek provedeme v sekci Flow data. Vybereme typ Head (1st kind) a zadáme tlakovou výšku hladiny, potom myší vybereme hrany elementů, na kterých mají podmínky platit (viz Obr.26).
51
Obr.26: FEFLOW - nastavení OKP
6.5 Vlastní simulace úlohy 6.5.1 Simulace ve Flow123D
Program Flow123D (viz kapitola 4.1) se spouští příkazem z příkazového řádku.
Pro správné proběhnutí programu je potřeba mít v souboru ini. nastaveny všechny potřené parametry a uvedeny správné cesty k souborům se sítí, OKP, zdrojem a sousednostmi. V souboru ini je důležité uvést hlavně cestu k potřebným souborům a nastavení výpočtů.
[Global]
Problem_type = 1 Description = test Stop_time = 1 //1275 Save_step = 0.010 Density_steps = 10 Density_on = No
[Input] - cesty ke vstupním souborům File_type = 1
Mesh = ./input/test.msh
52 //Source = ./input/test.src
Material = ./input/test.mtr Boundary = ./input/test.fbc Neighbouring = ./input/test.ngh
[Run]
Log_file = pokus.log - log soubor, zde je uložen popis procesu simulace Screen_verbosity = 7
Log_verbosity = 7 Pause_after_run = No Checkpoints = -9999
[Output] - nadefinování formátu výstupního souboru a místa uložení POS_format = VTK_SERIAL_ASCII
Write_output_file = Yes
Output_file = ./output/test.vtu -náš výstupní soubor Output_file_2 = ./output/test.pos
Output_digits = 12 Output_file_type = 1
53
6.5.2 Simulace ve FEFLOW
Výpočet probíhá v Run – start simulation. Zde je důležité si před samotným spuštěním nastavit nahrávání simulace do datového souboru (viz. Kapitola 4.2.1.5).
Po skončení simulace se nám zobrazí výsledek. Pomocí uživatelského prostředí ho můžeme prohlížet a zkoumat. Nesmíme zapomenout vyexportovat soubor s hodnotami sledovaných veličin v uzlech elementů. Provedeme to v Halt and view results a uložíme hodnoty Head.
Obr.27: FEFLOW - zobrazení výsledku
Ze souborů Head.dat a ze záznamu dat v ASCII.dac si vytvoříme soubor se strukturou VTK. Postup převodu je popsán v kapitole 5.2.
54
6.6 Porovnání
Nyní máme dva soubory, které lze otevřít v ParaView. Pro porovnání použijeme již popsaný postup. Zobrazíme si už jen výsledný rozdíl tlaků ploch z obou programů.
Na prvním obrázku je síť z 50 elementů, na druhém z 300 elementů.
Obr.28: Výsledný rozdíl tlaků u sítí s 50 e1ementy, absolutní rozdíl je 2%
Vidíme, že největší rozdíl je asi o 2%.
Obr.29: Výsledný rozdíl tlaků u sítí s 300 elementy, absolutní rozdíl je 2,16%
Z obrázků je patrné, že výsledky z obou programů se docela shodují. Rozdíl je oproti hrubé síti větší, a to o 2,16%. Vidíme, že u 100x propustnější pukliny oproti okolí, je větší rozdíl výsledků než u méně vodivé pukliny. Dále pak je zřetelné, jak moc nám ovlivňuje výsledek jemnost sítě. Při jemnější síti se na většině plochy ukazuje menší rozdíl výsledků (více je zobrazen modrý odstín barvy).
55
7 Úlohy
7.1 Úloha 1
Jako první ukázkovou porovnávací úlohu máme vytvořenou na čtvercové oblasti 10x10metrů. Máme v ní dvě kvádrovité oblasti 2x8metru s nižší propustností prostředí a OKP jsou zadány pouze v rozích a nejsou zde žádné pukliny.
Obr.30: Zadání první úlohy
p=100
k=1
k=0,01
k=0,01
10m
2m 2m 2m 2m
2m
8 m
8 1 m
0 m
p=120
56
Obr.31: Porovnání první úlohy, absolutní rozdíl je 0,9%
Z obrázku je patrné, že největší rozdíly jsou na hranách prostředí s různou propustností o 0,9% a pak hlavně v levé části oblasti. Jinak jsou výsledky úloh shodné.
7.2 Úloha 2
Druhá úloha byla vybrána na ukázku toho, co se stane, když se zadá puklina mimo hranice elementů. Přes naši oblast byly vytvořeny pukliny a OKP byly zadány přímo do pukliny. Dole byla zadána tlaková výška 120metrů a nahoře 100metrů.
Z obrázku je patrné, že spodní OKP program v pořádku přijal a počítá s ní. Ale u pukliny nahoře, která nekopíruje hranice elementů je evidentně vidět, že má problém se zadanou okrajovou podmínkou i s vypočítanými hodnotami. Proto pokud bychom chtěli použít takto rozestavené pukliny, musíme použít jinou strukturu sítě, tak aby puklina procházela skrz hrany elementů.
57
Obr.32: Výsledek druhé úlohy
7.3 Úloha 3
Jako další úloha byla vytvořena oblast s puklinami. Využili jsme ji pro několik různých simulací, kde se měnila místa zadání OKP a propustnost puklin, kde za prvé mají pukliny propustnost 100x větší než ostatní materiál a za druhé 100x menší.
Obr.33: Zadání třetí úlohy
10m 5
m
2m
[3,5;1,5]
[4,5;1,5]
[6,5;4]
[3,5;3]
[5;4,5]
[6;4,5]
6,5m
58
7.3.1 Úloha 3 – zadání OKP
Nejprve byly zadány OKP tlaková výška 120metrů na celé pravé straně a 100 metrů na levé. Vodivost pukliny v prvním případě je k=100.
Obr.34: Výsledek třetí úlohy s velkou vodivostí puklin, absolutní rozdíl je asi 3%
Obr.35: Výsledek třetí úlohy s malou vodivostí puklin, absolutní rozdíl je 2,3%
Z výsledků víme, že největší rozdíly jsou na hraně pukliny, kde je změna proudění.
Vetší rozdíl je pak v situaci, kdy je puklina vodivější a svádí proudění v oblasti do sebe.
Dále pak je vidět, že oblasti, kde jsou nastaveny OKP, jsou skoro shodné.
59
7.3.2 Úloha 3 – jiné OKP
Nyní na stejnou úlohu aplikujeme jiné OKP a to tak, že je přímo zadáme na hranici našich oblastí do místa, kde začíná puklina.
Obr.36: Výsledek třetí úlohy s velkou vodivostí puklin OKP na puklině, absolutní rozdíl je 4,9%
Obr.37: Výsledek třetí úlohy s malou vodivostí puklin OKP na puklině, absolutní rozdíl je 10,4%
Při použití OKP přímo v bodě na hraně oblasti a pukliny vidíme, že největší rozdíl je právě tam, kde jsou podmínky zadány. Je to způsobeno tím, že každý program různě rozpočítává OKP na další elementy. Tím pádem trošku lepšího výsledku můžeme dosáhnout jemnější sítí v okolí OKP. Ale jinak výsledek celé oblasti je velmi podobný.
60
Závěr
Hlavním výsledkem práce je metodika nastavení shodné úlohy u stacionárního proudění podzemní vody v puklinovém prostředí v programech Flow123D a FEFLOW.
Nejprve bylo nutné seznámit se programy. Důležité bylo zjistit, kde, jak a v jakých veličinách se zadávají vstupní parametry. Pak byl vybrán vizualizační program ParaVIEW, kromě podpory velkého množství formátu dat je i velmi přehledný. Abych mohl v ParaVIEW data porovnávat, tak bylo potřeba prostudovat výstupní data z obou simulačních programů. Program Flow123D lze rovnou otevřít bez dalších úprav.
Program FEFLOW sice umožňuje uložit jeden výstupní soubor, ale s ním umí pracovat pouze programy od stejného výrobce. Ale umí vyexportovat samostatnou strukturu sítě a pak zvlášť výsledná data. Výstupní data jsem prozkoumal. Podíval jsem se na podporu dat ve vizualizačním programu. Jako nejvhodnější mi přišlo, použít strukturu VTK.
V tomto uspořádání dat, jsou přehledně zapsány uzly, nody a hodnoty dat. Pak jsem nalezl postup a úpravy dat, potřebné k převodu. Abych zjednodušil a zrychlil porovnávání, tak jsem vytvořil program, který snadno a rychle udělá tuto zdlouhavou práci za nás. Dále v ParaVIEW bylo potřeba pomocí různých filtrů navrhnout postup, který by umožňoval srovnat data z různých sítí. Pomocí několika kroků byla sestavena posloupnost kroků, která nám zaručovala porovnání výsledků, ze stejně velkých oblastí, ale s rozdílnou sítí. Tato posloupnost kroků lze v programu uložit a pojmenovat jako nový filtr. Díky tomu při dalším porovnávání stačí vybrat pouze námi vytvořený filtr.
Předposlední částí práce bylo ověření vymyšleného postupu na jednoduchých úlohách.
Když bylo ověřeno, že úlohy jsou zadány shodně a rozdíl výsledků je nulový, tak se přešlo ke složitějším úlohám. Tyto demonstrační úlohy byly zvoleny, tak abych na nich mohl zadat různé situace, ve kterých se předvídalo, že bude vidět nějaký rozdíl.
Z výsledů je patrné, že rozdíly byli maximálně v řádu pár procent a nacházeli se na hranicích puklin nebo na hranicích různých propustností materiálů. Největší rozdíl pak byl u úlohy, kde se OKP zadaly přímo do pukliny. Je to způsobeno tím, že každý program rozpočítává od bodu OKP hodnoty na další elementy trochu jinak a proto se dá částečně eliminovat větší jemností sítě.
Navržená metodika může posloužit k porovnání přesností výsledků z programu Flow123D (vyvíjeného na TUL) vůči komerčnímu systému. Nebo lze práci použít, pro seznámení se s problematikou ustáleného proudění v puklinovém prostředí.
61
Literatura
[1] Český hydrometeorologický ústav oddělení povrchových vod. Oběh vody: The Water Cycle [online]. 2011 [cit. 2009-05-13]. Dostupný z WWW:
<http://ga.water.usgs.gov/edu/watercycleczech.html#freshstorage>.
[2] Hladiny podzemní vody [online]. 29.6.2008 , 16.5.2011 [cit. 2009-05-25].
Dostupný z WWW: < http://cs.wikipedia.org/wiki/Zvode%C5%88>.
[3] Uživatelská dokumentace k programům Flow123D a Eclipse.
[4] Česká zemědělská univerzita v Praze. Modelování odtoku z povodí pomocí Boussinesqovy rovnice. [online]. 2011 [cit. 2011-10-22]. Dostupný z WWW:
< http://www.kvhem.cz/wp-content/uploads/2009/02/Kajtarova.pdf
[5] Hladiny podzemní vody [online]. 29.6.2008 , 16.5.2009 [cit. 2011-11-15].
Dostupný z WWW: < http://cs.wikipedia.org/wiki/Zvode%C5%88>.
[6] Podzemní voda [online]. 29.6.2008 , 16.10.2011 Dostupný z WWW:
< http://geotech.fce.vutbr.cz/studium/geologie/skripta/podvoda.htm>.
[7] KRÁLOVCOVÁ, Jiřina. Proudění a transport látek v různých typech hostitelské horniny s různou geologickou stavbou. 2009. Liberec, Dostupný z WWW:
<http://artec.tul.cz/?content=upload/SIMONA09_kralovcova.pdf&lang=cs>.
[8] ULEBERG and Kleppe, Dual Porosity, Dual Permeability Formulation for Fractured Reservoir Simulation. 1996. Savanger. Dostupný z WWW: <
http://www.ipt.ntnu.no/~kleppe/TPG4150/fracturedpaper.pdf>.
