MASARYKOVA UNIVERZITA V BRNĚ. vody. Mgr. Eva Cvachovcová

MASARYKOVA UNIVERZITA V BRNEˇ PRˇI´RODOVEˇDECKA´ FAKULTA

Matematicky´ model proudeˇnı´ podzemnı´ vody DIPLOMOVA´ PRAĆE

Mgr. Eva Cvachovcova´

Brno, jaro 2009

ii

Prohla´sˇenı´ Prohlasˇuji, zˇe tato diplomova´ praće je my´m pu˚vodnı´m autorsky´m dı´lem, ktere´ jsem vypracovala samostatneˇ. Vsˇechny zdroje, prameny a literaturu, ktere´ jsem prˇi vypracovańı´ pouzˇ´ıvala nebo z nich cˇerpala, v praći rˇa´dneˇ cituji s uvedenı´m u´plne´ho odkazu na prˇ´ıslusˇny´ zdroj.

Vedoucı´ praće: Mgr. Jirˇ´ı Zelinka, Dr. iii

iv

Podeˇkovańı´ Na prvnı´m mı´steˇ bych chteˇla podeˇkovat sve´mu vedoucı´mu, Mgr. Jirˇ´ımu Zelinkovi, Dr., za vsˇechen cˇas i uzˇitecˇne´ rady, ktere´ mi ochotneˇ veˇnoval beˇhem cele´ me´ praće. Da´le bych chteˇla podeˇkovat spolecˇnosti GEOtest Brno, a.s. za spolupraći a poskytnutı´ dat pro mou praći, a sve´mu manzˇelovi za pomoc nejen prˇi hledańı´ gramatickyćh chyb.

v

vi

Obsah ´ vod . . . . . . . . . . . . . . . . . . . . . . . . . . U 1 Vytvorˇenı´ matematicke´ho modelu . . . . . . . 1.1 Definice za´kladnıćh pojmu˚ . . . . . . . . . . 1.2 Koncepcˇnı´ model . . . . . . . . . . . . . . . 1.3 Matematicky´ model . . . . . . . . . . . . . . 1.3.1 Rovnice kontinuity . . . . . . . . . 1.3.2 Darcyho za´kon . . . . . . . . . . . 1.3.3 Rovnice proudeˇnı´ podzemnı´ vody 1.3.4 Okrajove´ a pocˇa´tecˇnı´ podmıńky . 2 Metody rˇesˇenı´ modelu . . . . . . . . . . . . . . 2.1 Metoda konecˇnyćh diferencı´ . . . . . . . . . 2.2 Metoda konecˇnyćh prvku˚ . . . . . . . . . . . 2.3 Srovnańı´ metod . . . . . . . . . . . . . . . . 3 Model oblasti Slavkov u Brna . . . . . . . . . 3.1 Koncepcˇnı´ model . . . . . . . . . . . . . . . 3.2 Rˇesˇenı´ a kalibrace modelu . . . . . . . . . . . Literatura . . . . . . . . . . . . . . . . . . . . . . . .

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 7 7 9 10 11 12 13 16 19 19 24 28 31 34 34 49

1

Seznam obra´zku˚ 1

Rekultivovana´ skla´dka.

1.1 1.2 1.3 1.4

Prˇ´ıklad typu zvodneˇ. 8 Rozdı´l ve vrtu volne´ a napjate´ hladiny. 8 Objemovy´ pru˚tok a filtracˇnı´ rychlost. 9 Darcyho pokus. 12

2.1 2.2

Prˇ´ıklad sı´teˇ s vyznacˇenı´m sousednıćh uzlu˚. 20 Prˇ´ıklad vy´beˇru uzlu˚ ve smeˇru norma´ly pro neprˇesneˇ i prˇesneˇ zadanou hranici. 22 Oznacˇenı´ indexu˚. 23 Prˇ´ıklad triangulace. 27 Ba´zova´ funkce ϕA . 28 Hranicˇnı´ ba´zova´ funkce ϕB . 28

2.3 2.4 2.5 2.6

6

3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14

Odebrane´ vzorky zemin. 32 Vrtna´ souprava. 33 Mapa za´jmove´ oblasti s rozdeˇlenı´m jizˇnı´ a severnı´ cˇa´sti. 36 Privilegovane´ cesty proudeˇnı´ podzemnı´ vody. 37 Mapa nameˇrˇenyćh izohydrohyps. 38 Geologicky´ rˇez za´jmovou oblastı´. 39 Sı´t’modelu s vyznacˇenı´m okrajovyćh podmıńek. 40 Prvnı´ odhad koeficientu˚ propustnosti. 41 Prvnı´ modelovane´ izohydrohypsy. 42 Konecˇny´ stav modelovanyćh izohydrohyps. 43 Konecˇne´ nastavenı´ koeficientu˚ propustnosti. 44 Uka´zka pru˚beˇhu vy´pocˇtu. 45 Vizualizace modelovane´ hladiny podzemnı´ vody. 46 Vizualizace modelovane´ hladiny podzemnı´ vody, nepropustne´ho podlozˇ´ı a pozice vrtu˚. 46 3.15 Vizualizace modelovane´ hladiny podzemnı´ vody a zemske´ho povrchu. Tam, kde voda prˇesahuje povrch, je hladina napjata´. 47 3.16 Modelovany´ rozsah kontaminace. 48

3

´ vod U Podzemnı´ voda zahrnuje vesˇkerou vodu, ktera´ se vyskytuje pod zemsky´m povrchem, tvorˇ´ı jeden z nejvy´znamneˇjsˇ´ıch zdroju˚ pitne´ vody na nasˇ´ı planeteˇ. V Cˇeske´ republice je to asi 44 % za´sob pitne´ vody [La´gner]. V minulosti vsˇak docha´zelo zejmeńa dı´ky nedbale´mu zacha´zenı´ s nebezpecˇny´mi la´tkami prˇi pru˚myslove´ vy´robeˇ cˇasto k jejı´mu znecˇisˇt’ovańı´. Dnes jizˇ nasˇteˇstı´ plneˇ docenˇujeme jejı´ vy´znam ale take´ vy´znam zˇivotnı´ho prostrˇedı´. Dı´ky tomu docha´zı´ k odstranˇovańı´ mnoha staryćh ekologickyćh za´teˇzˇ´ı1 a uva´deˇnı´ prˇ´ırody do pu˚vodnı´ho stavu. Matematicky´ model je zjednodusˇeny´ popis rea´lne´ho syste´mu, ktery´ matematicky´mi prostrˇedky simuluje fungovańı´ tohoto syste´mu. Matematicke´ modelovańı´ proudeˇnı´ vody se vyuzˇ´ıva´ pro simulaci proudeˇnı´ vody v kanalizacˇnıćh syste´mech, vodovodnıćh sı´tıćh, v rˇekaćh (naprˇ´ıklad pro simulaci povodnˇovyćh uda´lostı´), ve vodohospoda´rˇskyćh stavbaćh a v neposlednı´ rˇadeˇ take´ k proudeˇnı´ podzemnı´ vody. Modelovańı´ proudeˇnı´ podzemnı´ vody se vyuzˇ´ıva´ pra´veˇ pro znecˇisˇteˇne´ syste´my. Umozˇnˇuje na´m jednak odhadnout rozsah a budoucı´ postup znecˇisˇteˇnı´, ale take´ navrhnout optima´lnı´ zpu˚soby jeho odstraneˇnı´. Dı´ky matematicke´mu modelu mu˚zˇeme simulovat ru˚zne´ postupy odstraneˇnı´ znecˇisˇteˇnı´ (naprˇ. cˇerpańı´ znecˇisˇteˇne´ vody, stavba podzemnı´ barie´ry, vyuzˇitı´ bakte´riı´ atd.), srovna´vat je a vybrat nejefektivneˇjsˇ´ı. Proto se matematicke´ modelovańı´ stalo nedı´lnou soucˇa´stı´ remedikace2 nejen na nasˇem u´zemı´. Tato praće poda´va´ jak teoreticky´ tak prakticky´ pohled na matematicke´ modelovańı´ proudeˇnı´ podzemnı´ vody. Je koncipovańa tak, aby byla dobrˇe srozumitelna´ i pro cˇtena´rˇe neznale´ho dane´ problematiky at’ jizˇ z matematicke´ho cˇi hydrogeologicke´ho pohledu. Mu˚zˇe slouzˇit jak studentu˚m matematicke´ho modelovańı´ jako prakticka´ uka´zka, tak studentu˚m hydrologie k sezna´menı´ se s matematickou strańkou problematiky. V u´vodnı´ kapitole te´to praće jsou definovańy za´kladnı´ pojmy, ktere´ se vztahujı´ k hydrogeologii a ktere´ jsou da´le pouzˇ´ıvańy. V dalsˇ´ım je popsań postup prˇi vytva´rˇenı´ matematicke´ho modelu. Je zde z rovnice kontinuity a Darcyho za´kona odvozena za´kladnı´ rovnice proudeˇnı´ podzemnı´ vody. Na´sleduje prˇehledna´ tabulka a popis jednotlivyćh typu˚ rovnic, ktere´ se vyuzˇ´ıvajı´ prˇi modelovańı´ proudeˇnı´ podzemnı´ vody spolu s pocˇa´tecˇnı´mi podmıńkami. Druha´ kapitola popisuje dveˇ v soucˇasnosti nejpouzˇ´ıvaneˇjsˇ´ı numericke´ metody pro rˇesˇenı´ modelu proudeˇnı´ podzemnı´ vody, metodu konecˇnyćh diferencı´ a metodu konecˇnyćh prvku˚. Praće se zameˇrˇuje zejmeńa na aspekty prakticke´ aplikace teˇchto metod na model proudeˇnı´ podzemnı´ vody a na jejich srovnańı´ z hlediska pouzˇitı´. Poslednı´ kapitola popisuje pouzˇitı´ modelu proudeˇnı´ podzemnı´ vody v prˇ´ıpadove´ 1. Odborny´ termıń pro znecˇisˇteˇnı´ prˇ´ırodnı´ho prostrˇedı´, ktere´ vzniklo v minulosti. 2. Uvedenı´ prˇ´ırody do pu˚vodnı´ho stavu.

5

studii na oblasti Slavkov u Brna. Tato cˇa´st je doplneˇna graficky´m zpracovańı´m pro lepsˇ´ı prˇedstavu o zkoumane´ problematice.

Obra´zek 1: Rekultivovana´ skla´dka.

6

Kapitola 1

Vytvorˇenı´ matematicke´ho modelu 1.1 Definice za´kladnıćh pojmu˚ V te´to cˇa´sti se sezna´mı´me se za´kladnı´mi pojmy a procesy, ktere´ se ty´kajı´ podzemnı´ vody a jejı´ho proudeˇnı´. Proudeˇnı´m podzemnı´ vody se zaby´va´ samostatny´ obor hydrauliky (cˇa´st technicke´ mechaniky, ktera´ studuje za´kony klidu a pohybu kapalin). Veˇtsˇina podzemnı´ vody vznika´ pru˚sakem sra´zˇkove´ vody pod zemsky´ povrch. Na rychlost a mnozˇstvı´ prosakovańı´ a proudeˇnı´ vody ma´ vliv zejmeńa propustnost horniny. Propustnost je schopnost po´rovite´ho prostrˇedı´ propousˇteˇt vodu. Je-li propustnost vztazˇena k proudeˇnı´ kapaliny o urcˇityćh vlastnostech (podzemnı´ voda), vyjadrˇuje se pomocı´ koeficientu filtrace (k, m s−1 ). V prˇ´ırodeˇ neexistuje hornina, ktera´ by byla absolutneˇ nepropustna´, neˇktere´ horniny vsˇak majı´ tak malou propustnost, zˇe se povazˇujı´ za nepropustne´ (naprˇ. mastne´ jı´ly). Naprˇ´ıklad pro sˇteˇrk se koeficient filtrace pohybuje zhruba od 10−2 do 1, pro jı´l meńeˇ nezˇ 10−7 . Laboratorneˇ se stanovuje propustnost propustomeˇrem, ve vrtech lze stanovit propustnost hydrodynamicky´mi zkousˇkami. K nim patrˇ´ı cˇerpacı´ zkousˇky, kdy se meˇrˇ´ı mnozˇstvı´ vody cˇerpane´ za sekundu a snı´zˇenı´ hladiny vody ve vrtu v za´vislosti na cˇase [PoVo]. Po´rovitost vyjadrˇuje pomeˇr mezi objemem volne´ho prostoru a jednotkovy´m objemem horniny (n) [HoNo]. Zemskou ku˚ru mu˚zˇeme rozdeˇlit do dvou oblastı´: nasycene´ a nenasycene´. Nasycena´ oblast znamena´, zˇe kazˇdy´ dostupny´ prostor je zaplneˇn vodou, zatı´mco v nenasycene´ zońeˇ jsou sta´le jesˇteˇ mı´sta vyplneˇna´ vzduchem, ktery´ mu˚zˇe by´t nahrazen vodou [WikFE]. Oblast pod zemı´, ktera´ je vyplneˇna´ vodou se oznacˇuje jako zvodenˇ (obr. 1.1). Hornı´ povrch zvodneˇ tvorˇ´ı hladinu podzemnı´ vody. V nasycenyćh zońaćh mu˚zˇe vzniknout tzv. napjata´ hladina podzemnı´ vody, kde je tlak vody veˇtsˇ´ı nezˇ atmosfe´ricky´. Pokud navrta´me mı´sto s napjatou hladinou, voda pod tlakem stoupa´ (azˇ do vy´sˇky oznacˇovane´ jako piezometricka´, viz nı´zˇe). Je-li tlak vody v nasycene´ zońeˇ roven atmosfe´ricke´mu, pak jejı´ hladinu oznacˇujeme jako hladinu volnou. (Po navrtańı´ se vy´sˇka hladiny nemeˇnı´, viz obra´zek 1.2.) V nenasycenyćh zońaćh je tlak vody mensˇ´ı nebo roven atmosfe´ricke´mu tlaku. Napjata´ hladina vznika´ obvykle tam, kde se voda dostane mezi dveˇ vrstvy nepropustne´ horniny, prˇ´ıpadneˇ je-li koeficient filtrace dane´ho prostrˇedı´ vy´razneˇ vysˇsˇ´ı nezˇ v okolı´. Vy´sˇka hladiny podzemnı´ vody se meˇrˇ´ı vzhledem ke stanovene´ srovna´vacı´ hladineˇ (obvykle 0 m nadmorˇske´ vy´sˇky). Nejcˇasteˇji se pouzˇ´ıva´ piezometricka´ vy´sˇka (H, m), cozˇ je 7

1. VYTVORˇENI´ MATEMATICKE´HO MODELU hladina podzemní vody zvodenˇ

nepropustné podloží

srovnávací hladina

Obra´zek 1.1: Prˇ´ıklad typu zvodneˇ. vy´sˇka hladiny podzemnı´ vody ve vrtu po usta´lenı´. Pro zna´zorneˇnı´ vy´sˇky a smeˇru proudeˇnı´ proudeˇnı´ podzemnı´ vody pouzˇ´ıva´me izohydrohypsy a proudnice. Izohydrohypsa je krˇivka spojujıćı´ mı´sta se stejnou vy´sˇkou hladiny podzemnı´ vody (je to vlastneˇ „vodnı´ vrstevnice“). Proudnice je trajektorie pohybu jednotlivyćh cˇa´stic kapaliny [WikFE].

volná hladina vrt

napjatá hladina zemský povrch

nenasycené prostředí výška hladiny podzemní vody = piezometrická výška

vrt

zemský povrch piezometrická výška výška hladiny podzemní vody nasycené prostředí

Obra´zek 1.2: Rozdı´l ve vrtu volne´ a napjate´ hladiny.

8

1. VYTVORˇENI´ MATEMATICKE´HO MODELU Podle za´vislosti na cˇase mu˚zˇeme proudeˇnı´ vody rozdeˇlit na usta´lene´ (staciona´rnı´) proudeˇnı´, kdy se rychlost proudeˇnı´ (ale i jine´ velicˇiny, naprˇ´ıklad hustota) v cˇase nemeˇnı´, a na neusta´lene´ (nestaciona´rnı´) proudeˇnı´, kdy se rychlost v cˇase meˇnı´ [WikOE]. Neusta´lene´ proudeˇnı´ je ovlivnˇovańo zejmeńa mnozˇstvı´m sra´zˇek souvisejıćı´m se zmeˇnou rocˇnıćh obdobı´ apod. V prˇ´ıpadeˇ usta´lenyćh proudeˇnı´ jsou tyto zmeˇny male´, nebo zkouma´me kra´tke´ cˇasove´ obdobı´, ve ktere´m se tyto zmeˇny vy´razneˇji neprojevı´. Horninne´ prostrˇedı´, ve ktere´m voda proudı´ vsˇemi smeˇry stejneˇ rychle, nazy´va´me izotropnı´ prostrˇedı´; v opacˇne´m prˇ´ıpadeˇ mluvı´me o anizotropnı´m prostrˇedı´. Abychom zı´skali izotropnı´ prostrˇedı´, museli bychom mı´t naprosto homogennı´ horninu se stejnou hustotou a velikostı´ po´ru˚ v cele´m objemu. Je zrˇejme´, zˇe v prˇ´ırodeˇ je takova´to situace velmi nepravdeˇpodobna´, proto se nejcˇasteˇji vyskytuje prostrˇedı´ anizotropnı´. Koeficient filtrace horniny, ktery´ odpovı´da´ rychlosti proudeˇnı´, se pak lisˇ´ı pro ru˚zne´ smeˇry proudeˇnı´. Mı´sto jednoho cˇ´ısla tedy pouzˇ´ıva´me tenzor filtrace (nejcˇasteˇji symetricka´ nebo diagona´lnı´ matice s rozmeˇry 3×3), ktery´ je slozˇen jednotlivyćh koeficientu˚ filtrace pro ru˚zne´ smeˇry. Popis anizotropnı´ho prostrˇedı´ je ovsˇem dost slozˇity´, a proto se cˇasto nahrazuje izotropnı´m. Rychlost proudeˇnı´ podzemnı´ vody se obvykle vyjadrˇuje filtracˇnı´ (Darcyovskou) rychlostı´ (q, m s−1 ). Skutecˇnou rychlost proudeˇnı´ nelze urcˇit, protozˇe voda neproudı´ v hornineˇ prˇ´ımo, ale obte´ka´ zrna a ru˚zne´ dalsˇ´ı prˇeka´zˇky. Filtracˇnı´ rychlost je tedy fiktivnı´ rychlost, skutecˇna´ dra´ha kapaliny je delsˇ´ı, a proto se filtracˇnı´ rychlost lisˇ´ı od skutecˇne´ rychlosti. Filtracˇnı´ rychlost je definovańa jako podı´l objemove´ho pru˚toku a celkove´ pru˚tocˇne´ plochy horniny [Matousˇek]. Objemovy´ (specificky´) pru˚tok je objem kapaliny, ktery´ protecˇe danou plochou za jednotku cˇasu (Q, m3 s−1 ). skutečná dráha vody průtočná plocha P Q objem vody proteklé za jednu sekundu Q

filtrační rychlost q=Q/P

Obra´zek 1.3: Objemovy´ pru˚tok a filtracˇnı´ rychlost. Meˇrna´ objemova´ za´sobnost (storativita) po´rovite´ho prostrˇedı´ vyjadrˇuje objem vody, ktery´ vytecˇe z jednotkove´ho objemu horniny prˇi jednotkove´m zvy´sˇenı´ piezometricke´ ˇ ´ıha]. Volneˇ rˇecˇeno jde o schopnost dane´ horniny udrzˇet vodu. vy´sˇky hladiny (S, m−1 ) [R

1.2 Koncepcˇnı´ model Schopnost spolehliveˇ prˇedpoveˇdeˇt rychlost a smeˇr proudeˇnı´ podzemnı´ vody je rozhodujıćı´ prˇi plańovańı´ a realizaci cˇisˇteˇnı´ podzemnı´ vody. Model mu˚zˇeme definovat jako 9

1. VYTVORˇENI´ MATEMATICKE´HO MODELU zjednodusˇenou verzi syste´mu existujıćı´ho v rea´lne´m sveˇteˇ (v nasˇem prˇ´ıpadeˇ syste´mu podzemnı´ vody), ktery´ prˇiblizˇneˇ simuluje du˚lezˇite´ vztahy akce a reakce v rea´lne´m sveˇteˇ. Protozˇe skutecˇny´ sveˇt je velmi slozˇity´, je potrˇeba jej zjednodusˇit. Zjednodusˇenı´ je soubor prˇedpokladu˚, ktere´ vyjadrˇujı´ povahu syste´mu a vlastnosti jejichzˇ chovańı´ je rozhodujıćı´ pro zkoumany´ proble´m. Neexistuje jednoznacˇny´ model pro dany´ syste´m podzemnı´ vody. Prvnı´ krok v procesu modelovańı´ je sestavenı´ koncepcˇnı´ho modelu, ktery´ popisuje dany´ syste´m a zjednodusˇujıćı´ prˇedpoklady. Koncepcˇnı´ model obvykle obsahuje: •

popis geometrie hranic zkoumane´ zvodneˇ;

•

druh horniny tvorˇ´ıcı´ zvodenˇ a jejı´ vlastnosti (koeficient filtrace, storativita, atd.);

•

druh proudeˇnı´ ve zvodni (t.j. jednodimenziona´lnı´, horizonta´lnı´ dvoudimenziona´lnı´, nebo trˇ´ıdimenziona´lnı´, usta´lene´ cˇi neusta´lene´);

•

vlastnosti vody (informace o homogeniteˇ, stlacˇitelnosti);

•

zvodenˇ s napjatou cˇi volnou hladinou;

•

zdroje a propady vody ve zkoumane´ domeńeˇ a na jejıćh hranicıćh;

•

pocˇa´tecˇnı´ podmıńky na zkoumane´ oblasti; a

•

okrajove´ podmıńky zkoumane´ oblasti, ktere´ vyjadrˇujı´ vztahy s jejı´m okolnı´m prostrˇedı´m.

Dalsˇ´ı krok v modelovacı´m procesu je vyja´drˇenı´ koncepcˇnı´ho modelu formou matematicke´ho modelu (soustavou rovnic). Matematicky´ model pak rˇesˇ´ıme analyticky´mi nebo numericky´mi metodami. Vy´sledek matematicke´ho modelu prˇina´sˇ´ı pozˇadovane´ prˇedpoveˇdi chovańı´ rea´lneˇ existujıćı´ho syste´mu.

1.3 Matematicky´ model V te´to cˇa´sti uvedeme za´kladnı´ rovnice hydrauliky a odvodı´me jednotlive´ typy modelu˚. Proudeˇnı´ podzemnı´ vody se rˇ´ıdı´ slozˇity´mi fyzika´lnı´mi principy, ktere´ zde ze zrˇejmyćh du˚vodu˚ budou uvedeny pouze strucˇnou a zjednodusˇenou formou. Prˇi tvorbeˇ modelu proudeˇnı´ podzemnı´ vody se vycha´zı´ z principu spojitosti a Darcyho za´kona. Princip spojitosti (kontinua) prˇedpokla´da´, zˇe voda spojiteˇ vyplnˇuje celou oblast, kterou prote´ka´. Proudeˇnı´ vody po´rovity´m prostrˇedı´m aproximujeme spojity´m fiktivnı´m proudeˇnı´m vody v cele´m tomto prostrˇedı´. Zanedba´me prostorove´ rozlozˇenı´ po´ru˚ a puklin v hornineˇ. Darcyho za´kon byl zformulovań Henrym Darcym v letech 1855-1856 na za´kladeˇ jeho pokusu˚ s vodou a pı´skem [WikFE]. Cı´lem modelu je urcˇit v modelovane´ oblasti hodnoty vektoru filtracˇnı´ rychlosti q ˇ ´ıha]. Tyto (vyjadrˇuje smeˇr a rychlost proudeˇnı´) a piezometrickou vy´sˇku hladiny H [R dveˇ velicˇiny na´m poskytujı´ u´plny´ popis stavu a proudeˇnı´ podzemnı´ vody. 10

1. VYTVORˇENI´ MATEMATICKE´HO MODELU 1.3.1 Rovnice kontinuity Vyćhozı´ rovnicı´ proudeˇnı´ podzemnı´ vody je rovnice kontinuity, ktera´ je odvozena z principu spojitosti a za´kona zachovańı´ hmotnosti: celkova´ zmeˇna hmotnosti vody obsazˇene´ v oblasti (∆m ) se musı´ rovnat soucˇtu hmotnosti vody vte´kajıćı´ (mi ) a hmotnosti zdroju˚ (ms ) v oblasti (naprˇ´ıklad sra´zˇky), mıńus hmotnost vody odte´kajıćı´ z oblasti (mo ) a hmotnost propadu˚ (md ) (naprˇ´ıklad odparˇovańı´ nebo cˇerpańı´ vody). ∆m = (mi − mo ) + (ms − md ).

(1.1)

Necht’ nada´le Ω znacˇ´ı modelovanou zvodenˇ a Λ je libovolna´ podoblast zvodneˇ Ω a ∂Λ jejı´ hranice. Pak rovnici (1.1) lze na podoblasti Λ zapsat takto [HoNo]: Z Z Z ∂ n ρ dV = − q ρ dS + P ρ dV ; (1.2) ∂t Λ

∂Λ

Λ

n = n(x, y, z, t) znacˇ´ı po´rovitost ([x, y, z] jsou prostorove´ sourˇadnice a t je cˇas), ρ = ρ(x, y, z, t) hustotu vody, q = (qx , qy , qz )T = q(x, y, x, t) je vektor filtracˇnı´ rychlosti a P = P (t) znacRˇ´ı celkovy´ objem vody ze zdroju˚ nebo propadu˚ Rna jednotkovy´ objem horniny. ∂ (Cˇlen ∂t nρ dV vyjadrˇujeRzmeˇnu hmotnosti v cˇase; − ∂Λ qρ dS je celkova´ hmotnost Λ vody, ktera´ protecˇe oblastı´; Λ P ρ dV vyjadrˇuje celkovou bilanci zdroju˚ a propadu˚.) Uzˇitı´m Gauss-Ostrogradske´ho veˇty dosta´va´me: Z Z q ρ dS = ∇(q ρ) dV. (1.3) Λ

∂Λ

Dosazenı´m (1.3) do (1.2), zameˇneˇnı´m porˇadı´ derivace a integra´lu u prvnı´ho cˇlenu (derivujeme podle jine´ promeˇnne´) a u´pravou dostaneme: Z ∂ (n ρ) + ∇ (q ρ) − P ρ dV = 0. (1.4) ∂t Λ

Rovnice (1.4) platı´ pro libovolnou oblast Λ ⊆ Ω, cozˇ je mozˇne´ jen tak, zˇe integrovana´ funkce je na Ω rovna nule. Dosta´va´me tak rovnici kontinuity: ∇(q ρ) +

∂(n ρ) = P ρ, ∂t

(1.5)

v rozepsane´m tvaru: ∂(qx ρ) ∂(qy ρ) ∂(qz ρ) ∂(n ρ) + + + = P ρ. ∂x ∂y ∂z ∂t

(1.6)

∂(n ρ) ∂H =Sρ , ∂t ∂t

(1.7)

ˇ ´ıha]: Da´le platı´ [R

11

1. VYTVORˇENI´ MATEMATICKE´HO MODELU kde S = S(x, y, z, t) je meˇrna´ objemova´ za´sobnost a H = (x, y, z, t) je piezometricka´ vy´sˇka hladiny podzemnı´ vody. Po dosazenı´ (1.7) do (1.5) dostaneme rovnici ve tvaru: ∇(q ρ) = −S ρ

∂H + P ρ, ∂t

(1.8)

1.3.2 Darcyho za´kon Darcyho za´kon je vztah, ktery´ popisuje proudeˇnı´ kapaliny pore´znı´m prostrˇedı´m. Uda´va´ linea´rnı´ za´vislost rychlosti proudeˇnı´ na koeficientu filtrace prostrˇedı´ a hydraulicke´m ˇ ´ıha]: gradientu [R h1 − h2 , L kde Q je objemovy´ pru˚tok, k je koeficient filtrace, P je plocha, kterou kapalina prote´ka´, 2 a h1 −h je hydraulicky´ gradient, viz obra´zek 1.4. L Q = −kP

L

h1 k

h2

Obra´zek 1.4: Darcyho pokus. 2 Vydeˇlenı´m rovnice plochou P dostaneme q = −k h1 −h , q je vektor filtracˇnı´ rychlosti. L V trojrozmeˇrne´m prostrˇedı´ platı´ zobecneˇny´ Darcyho vztah:

q = κ∇H; 

(1.9)



kxx kxy kxz  κ = κ(x, y, z) = kxy kyy kyz  je tenzor filtrace a ∇H je hydraulicky´ gradient. kxz kyz kzz V rozepsane´m tvaru: ∂H ∂H ∂H + kxy + kxz , ∂x ∂y ∂z ∂H ∂H ∂H qy = kxy + kyy + kyz , ∂x ∂y ∂z ∂H ∂H ∂H qz = kxz + kyz + kzz . ∂x ∂y ∂z Prˇi vysˇsˇ´ıch rychlostech proudeˇnı´ (zvy´sˇ´ı se objemovy´ pru˚tok) prˇesta´va´ v Darcyho vztahu platit linea´rnı´ za´vislost. Odklon od linearity je vsˇak velmi pozvolny´, a pro veˇtsˇinu ˇ ´ıha]. modelu˚ se proto zanedba´va´ [R qx = kxx

12

1. VYTVORˇENI´ MATEMATICKE´HO MODELU 1.3.3 Rovnice proudeˇnı´ podzemnı´ vody Spojenı´m rovnice kontinuity a Darcyho za´kona dostaneme rovnici neusta´lene´ho proudeˇnı´ stlacˇitelne´ kapaliny anizotropnı´m stlacˇitelny´m prostrˇedı´m. (V rovnici (1.8) nahradı´me q vyja´drˇenı´m (1.9) z Darcyho zobecneˇne´ rovnice.) ∇ (κ∇Hρ) = −Sρ

∂H + Pρ ∂t

(1.10)

V rozepsane´m tvaru: ∂ ∂x

∂H ∂H ∂ ∂H ∂H ∂H kxx ∂H + k + k ρ + k + k + k ρ + xy xz ∂z xy yy ∂y yz ∂z ∂x ∂y ∂y ∂x ∂ + ∂z kxz ∂H + kyz ∂H + kzz ∂H ρ = −Sρ ∂H + P ρ. ∂x ∂y ∂z ∂t

(1.11)

Rovnice (1.10), respektive (1.11), nejle´pe popisuje rea´lnou situaci proudeˇnı´ podzemnı´ vody. Je zrˇejme´, zˇe rˇesˇenı´ rovnice (1.10) je dost na´rocˇne´ a v praxi navıć cˇasto chybı´ dostatecˇne´ mnozˇstvı´ nameˇrˇenyćh dat. Podle pozˇadavku˚ na model a rozsahu dat tedy rovnici da´le zjednodusˇujeme. Prˇi neusta´lene´m proudeˇnı´ jsou velicˇiny v rovnici (1.10) (κ, H, ρ, . . . ) za´visle´ na cˇase. Pro velkou cˇa´st modelu˚ lze bez proble´mu˚ pouzˇ´ıt usta´lene´ proudeˇnı´. Rovnice se tak zjednodusˇ´ı o jednu promeˇnnou a navıć vypadne cˇlen Sρ ∂H . ∂t Zejmeńa u plosˇneˇ rozsa´hlyćh modelovanyćh oblastı´ je cˇasto obtı´zˇne´ stanovit hodnotu koeficientu filtrace pro jeden smeˇr proudeˇnı´ vody na cele´m u´zemı´, natozˇ pro vıće smeˇru˚. Tenzor filtrace se proto veˇtsˇinou uda´va´ pouze jako diagona´lnı´ matice s koeficienty filtrace pouze pro smeˇry sourˇadnyćh os. V rovnici tak vypadnou cˇleny s kij , i 6= j. Toto zjednodusˇenı´ je vhodne´ te´meˇrˇ pro vsˇechna prostrˇedı´ s vy´jimkou proudeˇnı´ v prasklinaćh (jako naprˇ´ıklad ska´ly), a proto budeme nada´le uva´deˇt rovnici pouze v tomto tvaru. U mnoha zvodnı´ navıć prˇevla´da´ proudeˇnı´ pouze jednı´m smeˇrem, a lze tedy pouzˇ´ıt koeficient filtrace pro tento smeˇr. V rovnici tak nahradı´me tenzor κ jednorozmeˇrny´m k. Dalsˇ´ı obvyklou idealizacı´ je neuvazˇovańı´ stlacˇitelnosti (deformace) prostrˇedı´ a kapaliny. V takove´m prˇ´ıpadeˇ je hustota vody ρ konstantnı´ a lze ji z rovnice vypustit. Nejcˇasteˇjsˇ´ım zjednodusˇenı´m mnoha modelu˚ je tzv. Dupuitu˚v teore´m (Dupuitu˚v prˇedpoklad). Tam, kde je vertika´lnı´ proudeˇnı´ podzemnı´ vody (tedy proudeˇnı´ ve smeˇru osy z) male´ nebo zanedbatelne´, platı´ ∂H = 0. ∂z Je tedy vhodne´ zetove´ sourˇadnice z rovnice zcela vypustit. Aby prˇitom byla dodrzˇena hmotnostnı´ rovnova´ha rovnice, je trˇeba prˇidat faktor vy´sˇky zvodneˇ b = b(x, y) (tedy sı´ly nasycene´ vrstvy). Vy´sˇka zvodneˇ za´visı´ na tom, zda se jedna´ o zvodenˇ s volnou, cˇi napjatou hladinou. V prˇ´ıpadeˇ napjate´ hladiny je voda stlacˇena mezi dveˇma ma´lo propustny´mi vrstvami (prˇ´ıpadneˇ relativneˇ ma´lo propustny´mi vu˚cˇi okolı´) a vy´sˇka zvodneˇ je tedy vzda´lenost nepropustnyćh vrstev a nenı´ za´visla´ na vy´sˇce hladiny podzemnı´ vody. 13

1. VYTVORˇENI´ MATEMATICKE´HO MODELU Naopak v prˇ´ıpadeˇ zvodneˇ s volnou hladinou je hornı´ hranicı´ zvodneˇ u´rovenˇ hladiny podzemnı´ vody. Vy´sˇku zvodneˇ lze tedy vyja´drˇit jako b = H − b0 := h, kde b0 je spodnı´ hranice zvodneˇ (tedy u´rovenˇ spodnı´ nepropustne´ vrstvy) a H je vy´sˇka hladiny podzemnı´ vody. Pak platı´, zˇe filtracˇnı´ rychlost je za´visla´ na b, tedy qx = kxx b ∂H , ∂x qy = kyy b ∂H . ∂y Dosazenı´m do rovnice (1.10) dostaneme: ∂ ∂H ∂ ∂H ∂H kxx b ρ + kyy b ρ = −Sρ + P ρ. ∂x ∂x ∂y ∂y ∂t Kombinacı´ jednotlivyćh aproximacı´ lze zı´skat na´sledujıćı´ typy rovnic: Neusta´lene´ proudeˇnı´ Stlacˇitelnost

Nestlacˇitelnost

Anizotropnı´ prostrˇedı´

Izotropnı´ Anizotropnı´ prostrˇedı´ prostrˇedı´

∇(κ∇Hρ) = −Sρ ∂H + Pρ ∂t

x

Izotropnı´ prostrˇedı´ ∇(∇H)k = −S ∂H +P ∂t

∇(κ∇H) = −S ∂H +P ∂t Usta´lene´ proudeˇnı´

Stlacˇitelnost

Nestlacˇitelnost

Anizotropnı´ prostrˇedı´

Izotropnı´ Anizotropnı´ prostrˇedı´ prostrˇedı´

∇(κ∇Hρ) = P ρ

x

∇(κ∇H) = P

Izotropnı´ prostrˇedı´ ∇(∇H)k = P

V rozepsane´m tvaru: Cˇa´st (a) je pro trojrozmeˇrne´ proudeˇnı´, cˇa´st (b) pro dvojrozmeˇrne´ proudeˇnı´ s napjatou hladinou a cˇa´st (c) pro dvojrozmeˇrne´ proudeˇnı´ s volnou hladinou. 1.

Usta´lene´ proudeˇnı´ stlacˇitelne´ kapaliny anizotropnı´m stlacˇitelny´m prostrˇedı´m: (a)

(b)

14

∂ ∂x

∂H ∂ ∂H ∂ ∂H kxx ρ + kyy ρ + kzz ρ = P ρ. ∂x ∂y ∂y ∂z ∂z ∂ ∂x

∂H ∂ ∂H kxx b ρ + kyy b ρ = P ρ. ∂x ∂y ∂y

(1.12)

1. VYTVORˇENI´ MATEMATICKE´HO MODELU (c) ∂ ∂x 2.

Usta´lene´ proudeˇnı´ nestlacˇitelne´ kapaliny anizotropnı´m nestlacˇitelny´m prostrˇedı´m: (a) ∂ ∂x

∂H kxx ∂x

(b)

∂ + ∂y

∂H kxx b ∂x

∂H kxx h ∂x

∂ ∂x

∂ ∂x

(c)

3.

∂H ∂ ∂H kxx h ρ + kyy h ρ = P ρ. ∂x ∂y ∂y

∂H kyy ∂y

∂ + ∂y

∂ + ∂y

∂ + ∂z

∂H kzz ∂z

∂H kyy b ∂y

∂H kyy h ∂y

= P.

= P.

(1.13)

= P.

Usta´lene´ proudeˇnı´ nestlacˇitelne´ kapaliny izotropnı´m nedeformujıćı´m se prostrˇedı´m: (a)

∂2H ∂2H ∂2H + + ∂x2 ∂y 2 ∂z 2

k=P

Je-li v te´to rovnici P = 0 (tj. oblast bez zdroju˚ nebo propadu˚), jedna´ se o Laplaceovu rovnici. (b) ∂ ∂x

∂ ∂x

(c)

4.

∂H b ∂x

∂H h ∂x

∂ + ∂x

∂ + ∂x

∂H b ∂y

∂H h ∂y

=

P k

=

P k

(1.14)

Neusta´lene´ proudeˇnı´ stlacˇitelne´ kapaliny anizotropnı´m stlacˇitelny´m prostrˇedı´m: (a) ∂ ∂x

∂H ∂ ∂H ∂ ∂H ∂H kxx ρ + kyy ρ + kzz ρ = −Sρ + P ρ. ∂x ∂y ∂y ∂z ∂z ∂t

(b) ∂ ∂x

∂ ∂x

(c)

∂H kxx bρ ∂x

∂H kxx hρ ∂x

∂ + ∂y

∂ + ∂y

∂H kyy bρ ∂y

∂H kyy hρ ∂y

= −Sρ

∂H + P ρ. ∂t

= −Sρ

(1.15)

∂H + P ρ. ∂t 15

1. VYTVORˇENI´ MATEMATICKE´HO MODELU 5.

Neusta´lene´ proudeˇnı´ nestlacˇitelne´ kapaliny anizotropnı´m nestlacˇitelny´m prostrˇedı´m: (a)

(b)

(c)

6.

∂ ∂x

∂ + ∂y

∂H kxx b ∂x

∂H kxx h ∂x

∂ ∂x

∂ ∂x

∂H kyy ∂y

∂ + ∂y

∂ + ∂y

∂ + ∂z

∂H kzz ∂z

∂H kyy b ∂y

∂H kyy h ∂y

= −S

∂H + P. ∂t

= −S

∂H + P. ∂t

= −S

∂H + P. ∂t

(1.16)

Neusta´lene´ proudeˇnı´ nestlacˇitelne´ kapaliny izotropnı´m nestlacˇitelny´m prostrˇedı´m: (a)

(b)

(c)

7.

∂H kxx ∂x

∂2H ∂2H ∂2H + + ∂x2 ∂y 2 ∂z 2

∂ ∂x

∂H b ∂x

∂ ∂x

∂H h ∂x

∂ + ∂y

∂ + ∂y

k = −S

∂H b ∂y

∂H h ∂y

∂H +P ∂t

k = −S

∂H +P ∂t

k = −S

∂H +P ∂t

(1.17)

Proudeˇnı´ stlacˇitelne´ kapaliny izotropnı´m stlacˇitelny´m prostrˇedı´m se obvykle nepouzˇ´ıva´.

1.3.4 Okrajove´ a pocˇa´tecˇnı´ podmıńky Necht’ ∂Ω znacˇ´ı hranici za´jmove´ oblasti Ω a ∂Ω1 ∪ ∂Ω2 ∪ ∂Ω3 = ∂Ω. Pak okrajove´ podmıńky mohou by´t na´sledujıćı´ [HoNo]: 1.

Jsou dańy hodnoty funkce H = H0 pro kazˇde´ (x, y, z) ∈ ∂Ω1 . (Dirichletova podmıńka, okrajova´ podmıńka prvnı´ho typu.)

2.

Je dań objemovy´ pru˚tok prˇes cˇa´st hranice −κ∇Hν = −qν = q0 pro kazˇde´ (x, y, z) ∈ ∂Ω2 , kde ν je jednotkovy´ vektor norma´ly k hranici. (Neumannova podmıńka, okrajova´ podmıńka druhe´ho typu.)

3.

Je dańa za´vislost objemove´ho pru˚toku na vy´sˇce hladiny pode´l hranice −κ∇Hν + H = −qν + H = q0 , pro kazˇde´ (x, y, z) ∈ ∂Ω3 , kde ν je jednotkovy´ vektor norma´ly k hranici a H je vy´sˇka hladiny pode´l hranice. (Smı´sˇena´ okrajova´ podmıńka, okrajova´ podmıńka trˇetı´ho typu.)

16

1. VYTVORˇENI´ MATEMATICKE´HO MODELU Podmıńka druhe´ho typu se pouzˇ´ıva´ zejmeńa pro hranice bez proudeˇnı´ (qν = 0). Podmıńka trˇetı´ho typu se pouzˇ´ıva´ pro polopropustne´ hranice, naprˇ´ıklad tvorˇ´ı-li hranici rˇeka, nebo jezero. Prˇi rˇesˇenı´ u´lohy neusta´lene´ho proudeˇnı´ je trˇeba kromeˇ okrajovyćh podmıńek zadat take´ pocˇa´tecˇnı´ podmıńku. Tou je piezometricka´ vy´sˇka hladiny H(x, y, z, 0) na za´jmove´ oblasti Ω v cˇase t = 0. (Pro usta´lene´ proudeˇnı´ nema´ pocˇa´tecˇnı´ podmıńka smysl.) Zvolena´ rovnice proudeˇnı´ podzemnı´ vody spolu s okrajovy´mi podmıńkami a prˇ´ıpadneˇ pocˇa´tecˇnı´ podmıńkou tvorˇ´ı matematicky´ model proudeˇnı´ podzemnı´ vody.

17

Kapitola 2

Metody rˇesˇenı´ modelu Je zrˇejme´, zˇe vzhledem ke slozˇitosti rovnic modelu je zejmeńa pro prakticke´ ućˇely nevyhnutelne´ numericke´ rˇesˇenı´. Nejvıće se pouzˇ´ıva´ metoda konecˇnyćh prvku˚ a metoda konecˇnyćh diferencı´. V dnesˇnı´ dobeˇ existuje neˇkolik programu˚, ktere´ nabı´zı´ kompletnı´ rˇesˇenı´ modelu proudeˇnı´ podzemnı´ vody pomocı´ teˇchto metod (MODFLOW, FLOWPATH, FEFLOW, GMS a dalsˇ´ı). Aplikace numerickyćh metod prˇi rˇesˇenı´ modelu proudeˇnı´ podzemnı´ vody ma´ oproti klasicky´m proble´mu˚m znacˇne´ odlisˇnosti. Prvnı´m du˚lezˇity´m faktorem je velikost zkoumane´ oblasti Ω, ktera´ se pohybuje v rˇa´du stovek metru˚ azˇ desı´tek kilometru˚. Je tedy zrˇejme´, zˇe jake´koli deˇlenı´ takove´to oblasti musı´ mı´t odpovı´dajıćı´ rozmeˇry (tedy obvykle se pohybuje v rˇa´du metru˚). S tı´m souvisı´ take´ prˇesnost u´daju˚, at’uzˇ se jedna´ o zemeˇpisne´ sourˇadnice cˇi nadmorˇske´ vy´sˇky. Druhy´m nemeńeˇ du˚lezˇity´m faktorem prˇi tvorbeˇ modelu jsou dostupne´ finance. Pro rˇesˇenı´ je nezbytne´ zı´skat geologicke´ a hydrologicke´ informace o za´jmove´ oblasti. Veˇtsˇina u´daju˚ se zı´ska´va´ z pru˚zkumnyćh vrtu˚ a hydrodynamicky´mi zkousˇkami. Spra´vneˇ provedeny´ vrt musı´ zasahovat azˇ do nepropustne´ho podlozˇ´ı (jinak by informace z neˇj zı´skane´ byly zkreslene´), cozˇ je cˇasto neˇkolik metru˚ pod zemı´ (obvykle kolem 10 azˇ 20 metru˚). Rozbor zeminy a provedenı´ hydrodynamickyćh zkousˇek (k nim patrˇ´ı zejmeńa ru˚zne´ cˇerpacı´ zkousˇky) pak prˇedstavujı´ dalsˇ´ı nezanedbatelne´ na´klady. Cena jednoho kompletnı´ho vrtu se pak pohybuje kolem 60 azˇ 100 000Kcˇ. V praxi jsou tedy informace o za´jmove´ oblasti pomeˇrneˇ spore´ a za´lezˇ´ı na odbornosti geologu˚ a zkusˇenosti modela´rˇe, zda zı´skana´ data spra´vneˇ interpretujı´, aby vytvorˇili funkcˇnı´ model. V neposlednı´ rˇadeˇ je trˇeba prˇihle´dnout take´ k faktoru˚m jako jsou inzˇeny´rske´ sı´teˇ v oblasti, budovy a soukrome´ pozemky, ktere´ mohou zı´ska´vańı´ dat znacˇneˇ komplikovat. V na´sledujıćı´m popı´sˇeme, jak vsˇechny tyto faktory ovlivnı´ metody rˇesˇenı´.

2.1 Metoda konecˇnyćh diferencı´ Princip metody konecˇnyćh diferencı´, zna´me´ take´ jako metoda sı´tı´, je velmi jednoduchy´. V oblasti, ve ktere´ hleda´me rˇesˇenı´, zvolı´me neˇjakou konecˇnou mnozˇinu bodu˚, kterou nazy´va´me sı´tı´ a prˇ´ıslusˇne´ body nazy´va´me uzly sı´teˇ. Pro rˇesˇenou diferencia´lnı´ rovnici a okrajove´ podmıńky nahradı´me derivace v jednotlivyćh uzlech sı´teˇ diferencˇnı´mi podı´ly (tj. linea´rnı´ kombinacı´ funkcˇnıćh hodnot hledane´ funkce v okolnıćh bodech). Dostaneme tak soustavu konecˇneˇ mnoha rovnic pro hodnoty hledane´ funkce v uzlech sı´teˇ. Vznikla´ 19

2. METODY RˇESˇENI´ MODELU soustava rovnic zrˇejmeˇ nenı´ linea´rnı´, pokud ani vyćhozı´ diferencia´lnı´ rovnice nebyla linea´rnı´. Volbu sı´teˇ prova´dı´me podle geometrie dane´ oblasti. Pro dvojrozmeˇrne´ proble´my nejcˇasteˇji volı´me obde´lnı´kovou nebo cˇtvercovou sı´t’, prˇ´ıpadneˇ lze zvolit take´ troju´helnı´kovou sı´t’. Pro trˇ´ırozmeˇrne´ proble´my ma´me sˇiroke´ pole mozˇnostı´ volby sı´teˇ, nejcˇasteˇji ovsˇem uzˇ´ıva´me kva´drove´ sı´teˇ. Uka´zˇeme si uzˇitı´ metody konecˇnyćh diferencı´ na dvoudimenziona´lnı´m proudeˇnı´ na rovnicıćh (1.17)(b) a (1.13)(c). Meˇjme na oblasti Ω dańu obde´lnı´kovou ekvidistantnı´ sı´t’ s krokem l ve smeˇru x a krokem k ve smeˇru y. Obde´lnı´k sı´teˇ o rozmeˇrech l × k budeme nazy´vat bunˇkou sı´teˇ. Uzly sı´teˇ nebudeme rozumeˇt vrcholy teˇchto obde´lnı´ku˚, jak je obvykle´, ale strˇedy(teˇzˇisˇteˇ) jednotlivyćh buneˇk sı´teˇ a budeme je znacˇit xi,j . Hodnoty hledane´ funkce H v jednotlivyćh uzlech budeme znacˇit hi,j . V prˇ´ıpadeˇ neusta´lene´ho proudeˇnı´ uvazˇujeme navıć cˇasovou promeˇnnou t a v rovnici vystupuje take´ derivace funkce H podle t. Zvolme tedy cˇasovy´ krok u > 0 a oznacˇme hi,j (t) hodnotu funkce H v uzlu xi,j v cˇase t. Dva uzly sı´teˇ nazveme sousednı´mi, je-li jejich vzda´lenost ve smeˇru osy x resp. y rovna l resp. k. Kazˇdy´ uzel sı´teˇ ma´ pra´veˇ cˇtyrˇi sousednı´ uzly. Lezˇ´ı-li pro dany´ uzel vsˇechny jeho sousednı´ uzly v oblasti Ω, pak se tento uzel nazy´va´ vnitrˇnı´ (viz obra´zek 2.1). Hranicı´ oblasti Ω obvykle rozumı´me uzavrˇenou krˇivku. V nasˇem prˇ´ıpadeˇ budeme jako hranici cha´pat sjednocenı´ buneˇk sı´teˇ pode´l te´to krˇivky. Prˇ´ıslusˇne´ uzly budeme nazy´vat hranicˇnı´. (Tato aproximace je zpu˚sobena zejmeńa plosˇnou rozsa´hlostı´ oblasti Ω a neprˇesnostı´ informacı´ o hranici, jak jizˇ bylo uvedeno v u´vodu te´to kapitoly.) Pro vnitrˇnı´ uzly lze derivace nahradit diferencˇnı´mi podı´ly obsahujıćı´mi hodnoty v sousednıćh uzlech. Pro hranicˇnı´ uzly vyuzˇijeme okrajovyćh podmıńek, abychom zı´skali stejny´ pocˇet rovnic jako uzlu˚ sı´teˇ. hraniční uzel vnitřní uzel oblast

Obra´zek 2.1: Prˇ´ıklad sı´teˇ s vyznacˇenı´m sousednıćh uzlu˚. ˇ esˇme rovnici (1.17)(b). Pro jednoduchost budeme prˇedpokla´dat, zˇe se jedna´ o zvoR denˇ s napjatou hladinou s konstantnı´ vy´sˇkou b. Dostaneme tak rovnici: 2 ∂ H ∂2H ∂H + kb = −S + P. ∂x2 ∂y 2 ∂t Pro vnitrˇnı´ uzel xi,j nahradı´me derivace vyskytujıćı´ se v rovnici diferencˇnı´mi podı´ly. Pro prostorove´ derivace lze uzˇ´ıt centra´lnı´ diference, pro cˇasovou derivaci uzˇijeme zpeˇtne´ 20

2. METODY RˇESˇENI´ MODELU diference:

∂2H (xi,j , t) ∂x2

≈

hi−1,j (t)−2hi,j (t)+hi+1,j (t) , l2

∂2H (xi,j , t) ∂y 2

≈

hi,j−1 (t)−2hi,j (t)+hi,j+1 (t) , k2

∂H (xi,j , t) ∂t

≈

hi,j (t)−hi,j (t−u) . u

Necht’ da´le Pi,j (t) znacˇ´ı objem vody ze zdroju˚ nebo propadu˚ a Si,j znacˇ´ı meˇrnou objemovou za´sobnost bunˇky sı´teˇ se strˇedem xi,j v cˇase t. Pro kazˇdy´ vnitrˇnı´ uzel xi,j v cˇase t pak dosta´va´me rovnici: hi−1,j (t)−2hi,j (t)+hi+1,j (t) hi,j−1 (t)−2hi,j (t)+hi,j+1 (t) + kb = l2 k2 −Si,j

hi,j (t)−hi,j (t−u) u

+ Pi,j (t).

Prˇepis okrajovyćh podmıńek provedeme na´sledovneˇ. Okrajova´ podmıńka prvnı´ho typu uda´va´ hodnotu hledane´ funkce H pode´l prˇ´ıslusˇne´ hranice. Zna´me tedy prˇ´ımo hodnoty hi,j (t) pro dane´ uzly. Prˇepis okrajove´ podmıńky druhe´ho typu −k

∂H = q0 ∂ν

je slozˇiteˇjsˇ´ı, nebot’podmıńka obsahuje derivaci ve smeˇru vneˇjsˇ´ı norma´ly k hranici. Pro aproximaci derivace v hranicˇnı´m uzlu xi,j pouzˇijeme jednosmeˇrne´ diference ∂H hi,j − ha,b ≈ , ∂ν d kde ha,b je hodnota v uzlu xa,b a d je vzda´lenost uzlu˚ xi,j a xa,b . Uzel xa,b vybereme podle prˇesnosti u´daju˚ z hranice takto: •

Urcˇenı´ hranice je dosti neprˇesne´ – vybereme sousednı´ uzel ve smeˇru x nebo y podle prˇiblizˇne´ho smeˇru norma´ly.

•

Urcˇenı´ hranice je pomeˇrneˇ prˇesne´ – vybereme nejblizˇsˇ´ı uzel smeˇru norma´ly (viz obra´zek 2.2). Podmıńku pak mu˚zˇeme nahradit rovnicı´ −k

hi,j (t) − ha,b (t) = qi,j (t), d

kde qi,j (t) je prˇ´ıslusˇna´ hodnota okrajove´ podmıńky. Poznamenejme jesˇteˇ, zˇe veˇtsˇina programu˚ postupuje prˇi vy´beˇru uzlu xa,b podle prvnı´ho bodu. 21

2. METODY RˇESˇENI´ MODELU n

n hraniční uzel

n

vybranný uzel

n

oblast

Obra´zek 2.2: Prˇ´ıklad vy´beˇru uzlu˚ ve smeˇru norma´ly pro neprˇesneˇ i prˇesneˇ zadanou hranici. Okrajovou podmıńku trˇetı´ho typu −k

∂H + H = q0 ∂ν

prˇepı´sˇeme podobneˇ na rovnici −k

hi,j (t) − ha,b (t) + hi,j (t) = qi,j (t), d

kde qi,j (t) znacˇ´ı intenzitu objemove´ho pru˚toku prˇes hranici bunˇky se strˇedem v xi,j v cˇase t a ha,b (t) znacˇ´ı hodnotu funkce H v uzlu xa,b , ktery´ vybereme analogicky jako v prˇedchozı´m prˇ´ıpadeˇ. Protozˇe se jedna´ o rovnici neusta´lene´ho proudeˇnı´, doplnı´me uvedene´ rovnice pocˇa´tecˇnı´ podmıńkou – vy´sˇkou hladiny H v cˇase t = 0. Dosta´va´me tak hodnoty hi,j (0) pro vsˇechny uzly sı´teˇ. Je-li n pocˇet uzlu˚, pak celkem ma´me soustavu n linea´rnıćh rovnic pro hranicˇnı´ a vnitrˇnı´ uzly v cˇase t a n pocˇa´tecˇnıćh hodnot v cˇase t = 0. Soustavu tedy rˇesˇ´ıme postupneˇ pro t = u, 2u, 3u, . . . . ˇ esˇme nynı´ rovnici (1.13)(c), prˇicˇemzˇ budeme prˇedpokla´dat, zˇe dno zvodneˇ ma´ R konstantnı´ vy´sˇku b0 = 0, od ktere´ budeme meˇrˇit hladinu podzemnı´ vody H. Pro vy´sˇku zvodneˇ pak platı´ h = H − 0 = H a rovnice je tedy tvaru ∂ ∂H ∂ ∂H kxx H + kyy H = P. ∂x ∂x ∂y ∂y a w(x, y) := kyy H ∂H . Pak derivace funkcı´ v a w lze aproxiOznacˇme v(x, y) := kxx H ∂H ∂x ∂y movat centra´lnı´mi diferencemi: v(x + 2l , y) − v(x − 2l , y) ∂v(x, y) ≈ , ∂x l w(x, y + k2 ) − w(x, y − k2 ) ∂w(x, y) ≈ . ∂y k 22

2. METODY RˇESˇENI´ MODELU Meˇjme nynı´ vnitrˇnı´ uzel xi,j a oznacˇme ki− 1 ,j koeficient filtrace mezi bunˇkami se strˇedy 2 xi,j a xi−1,j , podobneˇ ki+ 1 ,j koeficient filtrace mezi bunˇkami se strˇedy xi,j a xi+1,j . Ana2 logicky oznacˇme koeficienty filtrace ki,j− 1 a ki,j+ 1 . Da´le oznacˇme hi− 1 ,j , hi+ 1 ,j , hi,j− 1 a 2 2 2 2 2 hi,j+ 1 hodnoty hledane´ funkce H mezi prˇ´ıslusˇny´mi bunˇkami (viz obra´zek 2.3). 2

xi j−1 h ij−_12 xi−1j

k ij−_12

ki−_12 j

ki+ _12 j xi j

hi−_21 j hij+ _12

xi+1 j

hi+ _12 j

kij+ _12

xij+1

Obra´zek 2.3: Oznacˇenı´ indexu˚. Diference v a w v uzlu xi,j lze pak zapsat: ki+ 1 ,j hi+ 1 ,j v(x + 2l , y) − v(x − 2l , y) 2 2 → l

∂H(x+ 2l ,y) ∂x

− ki− 1 ,j hi− 1 ,j 2

2

∂H(x− 2l ,y) ∂x

l ∂H(x,y+ k )

2 ki,j+ 1 hi,j+ 1 − ki,j− 1 hi,j− 1 w(x, y + k2 ) − w(x, y − k2 ) ∂y 2 2 2 2 → k k Pro derivace funkce H opeˇt uzˇijeme centra´lnı´ diference:

∂H(x + 2l , y) H(x + l, y) − H(x, y) ≈ , ∂x l

,

∂H(x,y− k2 ) ∂y

.

∂H(x − 2l , y) H(x, y) − H(x − l, y) ≈ , ∂x l

∂H(x, y + k2 ) ∂H(x, y − k2 ) H(x, y + k) − H(x, y) H(x, y) − H(x, y − k) ≈ , ≈ . ∂y k ∂y k Dosadı´me-li vsˇechny diference do rˇesˇene´ rovnice dostaneme rovnici pro uzel xi,j : ki+ 1 ,j hi+ 1 ,j (hi+1,j − hi,j ) − ki− 1 ,j hi− 1 ,j (hi,j − hi−1,j ) 2

2

2

2

l2

+

ki,j+ 1 hi,j+ 1 (hi,j+1 − hi,j ) − ki,j− 1 hi,j− 1 (hi,j − hi,j−1 )

2 2 = Pi,j , k2 kde Pi,j znacˇ´ı zdroje nebo propady vody v bunˇce se strˇedem v xi,j . Hodnoty hi− 1 ,j , hi+ 1 ,j , 2 2 hi,j− 1 a hi,j+ 1 obvykle pocˇ´ıta´me jako aritmeticke´ pru˚meˇry sousednıćh hodnot, tj.:

+

2

hi− 1 ,j = 2

2

2

2

hi,j + hi−1,j hi+1,j + hi,j hi,j + hi,j−1 hi,j+1 + hi,j , hi+ 1 ,j = , hi,j− 1 = , hi,j+ 1 = . 2 2 2 2 2 2 2 23

2. METODY RˇESˇENI´ MODELU Rovnici lze tedy zapsat pouze s pouzˇitı´m hodnot v uzlech sı´teˇ: ki+ 1 ,j (h2i+1,j − h2i,j ) − ki− 1 ,j (h2i,j − h2i−1,j ) 2

2

2l2

+

ki,j+ 1 (h2i,j+1 − h2i,j ) − ki,j− 1 (h2i,j − h2i,j−1 ) 2

2

2k 2

=

= Pi,j . Rovnice pro vnitrˇnı´ uzly sı´teˇ je tedy nelinea´rnı´. Abychom zı´skali dostatecˇny´ pocˇet rovnic pro vsˇechny uzly sı´teˇ, prˇepı´sˇeme take´ okrajove´ podmıńky stejneˇ jako u prˇedchozı´ rovnice. Jedna´ se o usta´lene´ proudeˇnı´, cˇasova´ promeˇnna´ v rovnici nevystupuje, a tedy nenı´ trˇeba zada´vat pocˇa´tecˇnı´ podmıńku. Uka´zali jsme, jak lze uzˇitı´m metody konecˇnyćh diferencı´ prˇeve´st matematicky´ model na soustavu rovnic at’ jizˇ linea´rnıćh cˇi ne. Zu˚sta´va´ ale ota´zka, jak tuto soustavu rˇesˇit. U linea´rnıćh rovnic lze uzˇ´ıt klasicke´ metody jako je naprˇ´ıklad Gaussova eliminacˇnı´ metoda. Tuto metodu ale nelze uzˇ´ıt v prˇ´ıpadeˇ nelinea´rnıćh rovnic a i u linea´rnıćh rovnic klade tato metoda znacˇne´ na´roky na pameˇt’ pocˇ´ıtacˇe. Proto prˇi rˇesˇenı´ vznikle´ soustavy obvykle opeˇt volı´me vhodne´ numericke´ metody. Prˇecˇ´ıslujeme-li uzly sı´teˇ jednı´m indexem tak, zˇe postupujeme po rˇa´dcıćh sı´teˇ zleva doprava a odshora dolu˚, pak je matice soustavy blokoveˇ trˇ´ıdiagona´lnı´. V linea´rnı´m prˇ´ıpadeˇ lze uzˇ´ıt k jejı´mu rˇesˇenı´ metody pro pa´sove´ matice, ktere´ vy´pocˇet znacˇneˇ usnadnı´. Obecneˇ nejcˇasteˇji uzˇ´ıvany´mi metodami jsou ovsˇem iteracˇnı´ metody. Jako pocˇa´tecˇnı´ iteraci h0 volı´me bud’ hodnotu pocˇa´tecˇnı´ podmıńky, je-li u´loha nestaciona´rnı´, nebo prˇiblizˇne´ hodnoty funkce H zı´skane´ prˇ´ımy´m meˇrˇenı´m. Prˇi vy´pocˇtu pak postupujeme sı´tı´ (naprˇ. zleva doprava a shora dolu˚) a z rovnice soustavy pro kazˇdy´ uzel xi,j vypocˇ´ıta´me novou hodnotu hki,j pomocı´ hodnot v sousednıćh uzlech. Pro vy´pocˇet jedne´ nove´ hodnoty tak uzˇ´ıva´me jak nove´ sousednı´ hodnoty hka,b tak stare´ hodnoty hk−1 ´ riem pro zastavenı´ vy´pocˇtu je obvykle vztah c,d . Krite khk − hk+1 k < ε, kde k.k znacˇ´ı zvolenou normu a ε je pozˇadovana´ prˇesnost. Rychlost konvergence lze zlepsˇit uzˇitı´m va´zˇene´ iterace. Novy´m hodnota´m, ktere´ vystupujı´ ve vy´pocˇtu jednoho uzlu, da´me veˇtsˇ´ı va´hu nezˇ hodnota´m stary´m.

2.2 Metoda konecˇnyćh prvku˚ Metoda konecˇnyćh prvku˚ se v soucˇasnosti cˇasto pouzˇ´ıva´ k rˇesˇenı´ mnoha fyzika´lnıćh proble´mu˚. Objasnı´me strucˇneˇ za´kladnı´ postup metody, podrobneˇjsˇ´ı popis mu˚zˇe cˇtena´rˇ najı´t naprˇ´ıklad v [Rektorys]. Meˇjme Hilbertu˚v prostor L2 (Ω) se skala´rnı´m soucˇinem Z (f, g) = f (x)g(x) dx Ω

24

2. METODY RˇESˇENI´ MODELU a okrajovou u´lohu Au = f u = 0

na Ω na ∂Ω,

kde A je symetricky´ pozitivneˇ definitnı´ diferencia´lnı´ opera´tor definovany´ na podprostoru D funkcı´ z L2 (Ω), ktere´ jsou na hranici Ω rovny nule, f je funkce z L2 (Ω). Veˇta o minimu kvadraticke´ho funkciona´lu: Rovnice Au = f ma´ rˇesˇenı´ v ∈ D, pra´veˇ kdyzˇ funkciona´l F u = (Au, u) − 2(f, u) naby´va´ na D minima´lnı´ hodnoty ve v. Veˇta na´m da´va´ na´vod, jak prˇeve´st diferencia´lnı´ proble´m na proble´m variacˇnı´ho pocˇtu. Podmıńka pozitivnı´ definitnosti opera´toru A zajisˇt’uje jednoznacˇnost rˇesˇenı´ v. Je-li diferencia´lnı´ opera´tor elipticky´, pak je pozitivneˇ definitnı´. Rovnice usta´lene´ho proudeˇnı´ je elipticka´, a tato podmıńka je tedy splneˇna. Rovnice neusta´lene´ho proudeˇnı´ je ovsˇem parabolicka´. Metodu lze vyuzˇ´ıt i v tomto prˇ´ıpadeˇ, nicmeńeˇ rˇesˇenı´ parabolickyćh rovnic vede obecneˇ ke komplikacı´m, a proto se pro tento typ rovnice metoda konecˇnyćh prvku˚ prˇ´ılisˇ nepouzˇ´ıva´. Proble´m nasta´va´ take´ v prˇ´ıpadeˇ rovnice dvoudimenziona´lnı´ho proudeˇnı´ s volnou hladinou, kdy v rovnici vystupuje vy´sˇka zvodneˇ, ktera´ v prˇ´ıpadeˇ volne´ hladiny za´visı´ na hledane´ funkci H a rovnice nenı´ linea´rnı´. (Da´le se teˇmito prˇ´ıpady nebudeme zaby´vat.) Poznamenejme jesˇteˇ, zˇe veˇta nerˇ´ıka´ nic o tom, zda neˇjake´ rˇesˇenı´ vu˚bec existuje. Zvolme konecˇneˇdimenziona´lnı´ podprostor Dn prostoru D. Budeme hledat funkci vn ∈ Dn , ktera´ minimalizuje funkciona´l F na podprostoru Dn . Funkce v nenı´ obecneˇ prvkem Dn , a proto bude funkce vn pouze aproximacı´ rˇesˇenı´ v. Zvolme ba´zi ϕ1 , ϕ2 , . . . , ϕn podprostoru Dn . Pak kazˇdou funkci vn ∈ Dn lze zapsat vn = c1 ϕ 1 + c2 ϕ2 + · · · + cn ϕn . Hleda´me koeficienty c1 , . . . , cn ∈ R takove´, aby funkciona´l F vn = (Avn , vn ) − 2(f, vn ) = = (A (c1 ϕ1 + c2 ϕ2 + · · · + cn ϕn ) , c1 ϕ1 + c2 ϕ2 + · · · + cn ϕn )−2(f, c1 ϕ1 +c2 ϕ2 +· · ·+cn ϕn ) = = c21 (Aϕ1 , ϕ1 ) + 2c1 c2 (Aϕ1 , ϕ2 ) + · · · + 2c1 cn (Aϕ1 , ϕn ) + c22 (Aϕ2 , ϕ2 ) + 2c2 c3 (Aϕ2 , ϕ3 ) + . . . · · · + c2n (Aϕn , ϕn ) − 2c1 (f, ϕ1 ) − 2c2 (f, ϕ2 ) − · · · − 2cn (f, ϕn ) naby´val minima´lnı´ hodnoty. Musı´ tedy platit rovnice ∂F vn =0 ∂c1

∧

∂F vn =0 ∂c2

∧

...

∧

∂F vn = 0, ∂cn 25

2. METODY RˇESˇENI´ MODELU tj. dosta´va´me Ritzovu soustavu rovnic ([Rektorys]): 2c1 (Aϕ1 , ϕ1 ) + 2c2 (Aϕ1 , ϕ2 ) + · · · + 2cn (Aϕ1 , ϕn ) = 2(f, ϕ1 ), 2c2 (Aϕ1 , ϕ2 ) + 2c2 (Aϕ2 , ϕ2 ) + · · · + 2cn (Aϕ2 , ϕn ) = 2(f, ϕ2 ), .. .

(2.1)

2cn (Aϕ1 , ϕn ) + 2c2 (Aϕ1 , ϕn ) + · · · + 2cn (Aϕn , ϕn ) = 2(f, ϕn ). Uvazˇme naprˇ´ıklad rovnici (1.13)(a) resp. (b): ∂ ∂H ∂ ∂H ∂ ∂H kxx + kyy + kzz =P ∂x ∂x ∂y ∂y ∂z ∂z resp. ∂ ∂x

∂ Au = ∂x

∂H kxx b ∂x

∂ + ∂y

∂H kyy b ∂y

=P

Pak je ∂u kxx ∂x

∂ + ∂y

∂u kyy ∂y

∂u kxx b ∂x

∂ + ∂y

∂ + ∂z

∂u kzz ∂z

resp. ∂ Au = ∂x

∂u kyy b ∂y

a f = P . Pro koeficienty (Aϕi , ϕj ) bude platit: Z ∂ ∂ϕi ∂ ∂ϕi ∂ ∂ϕi kxx + kyy + kzz ϕj dx (Aϕi , ϕj ) = ∂x ∂x ∂y ∂y ∂z ∂z Ω

resp. Z (Aϕi , ϕj ) =

∂ ∂x

∂ϕi kxx b ∂x

∂ + ∂y

∂ϕi kyy b ∂y

ϕj dx.

Ω

Uzˇitı´m integrace per partes dostaneme Z Z Z ∂ ∂ϕi ∂ϕi ∂ϕi ∂ϕj kxx ϕj dx = kxx nx ϕj dS − kxx dx, ∂x ∂x ∂x ∂x ∂x Ω

∂Ω

Ω

kde n = (nx , ny , nz ) je jednotkovy´ vektor vneˇjsˇ´ı norma´ly k hranici ∂Ω. Funkce ϕi lezˇ´ı v Dn , cozˇ je podprostor D. D je podprostor funkcı´ z L2 (Ω), ktere´ jsou na hranici Ω rovny nule. Tedy funkce ϕi jsou na hranici Ω rovny nule, a tudı´zˇ platı´ Z ∂ϕi kxx nx ϕj dS = 0. ∂x ∂Ω

26

2. METODY RˇESˇENI´ MODELU Celkem tedy platı´ Z (Aϕi , ϕj ) = −

kxx

∂ϕi ∂ϕj ∂ϕi ∂ϕj ∂ϕi ∂ϕj + kyy + kzz dx. ∂x ∂x ∂y ∂y ∂z ∂z

Ω

Respektive analogicky pro prˇ´ıpad (b) Z ∂ϕi ∂ϕj ∂ϕi ∂ϕj + kyy dx. (Aϕi , ϕj ) = − kxx ∂x ∂x ∂y ∂y Ω

Princip metody konecˇnyćh prvku˚ spocˇ´ıva´ v tom, zˇe se cela´ oblast Ω rozdeˇlı´ na navza´jem disjunktnı´ cˇa´sti – prvky. V trˇ´ırozmeˇrnyćh oblastech nejcˇasteˇji volı´me cˇtyrˇsteˇny, prˇ´ıpadneˇ ru˚zne´ hranoly. Pro dvourozmeˇrne´ proble´my deˇlı´me oblast nejcˇasteˇji na troju´helnı´ky, lze ale take´ uzˇ´ıt cˇtyrˇu´helnı´ky a sˇestiu´helnı´ky. Funkce ϕ1 , . . . , ϕn se volı´ tak, aby co nejvıće koeficientu˚ (Aϕi , ϕj ) bylo rovno nule, tedy aby matice soustavy (2.1) byla pa´sova´ nebo velmi rˇ´ıdka´. Nejcˇasteˇji volı´me ru˚zne´ splinove´ funkce, jejichzˇ nosicˇ je omezeny´ na maly´ pocˇet prvku˚. Prˇ´ıklad, jak lze tyto funkce volit, uka´zˇeme pro dvojrozmeˇrny´ proble´m.

Obra´zek 2.4: Prˇ´ıklad triangulace. Meˇjme oblast Ω rozdeˇlenou na konecˇny´ pocˇet troju´helnı´ku˚ (viz obra´zek 2.4). Vrcholy troju´helnı´ku˚ budeme nazy´vat uzly. V dobrˇe sestrojene´ sı´ti by se nemeˇl prˇ´ılisˇ meˇnit pocˇet troju´helnı´ku˚ sty´kajıćıćh se v jednom uzlu (v nasˇem prˇ´ıpadeˇ sˇest, tedy kazˇdy´ vnitrˇnı´ uzel ma´ sˇest sousednıćh uzlu˚). Da´le by v sı´ti nemeˇlo by´t prˇ´ılisˇ mnoho prota´hlyćh troju´helnı´ku˚. (Pomeˇr nejveˇtsˇ´ı a nejmensˇ´ı vy´sˇky troju´helnı´ku by nemeˇl by´t veˇtsˇ´ı nezˇ cˇtyrˇi [Kazda].) Takovouto sı´t’lze sestrojit naprˇ´ıklad pomocı´ Delaunayovy triangulace. Ta maximalizuje nejmensˇ´ı u´hel sestrojenyćh troju´helnı´ku˚ a prˇitom ji lze algoritmicky pomeˇrneˇ rychle sestrojit (v cˇase O(n log n) [WikFE], kde n je pocˇet uzlu˚). Uzly v sestrojene´ sı´ti postupneˇ ocˇ´ıslujeme (postupujeme naprˇ. zleva doprava a shora dolu˚), prˇicˇemzˇ si pro kazˇdy´ uzel musı´me pamatovat vsˇechny sousednı´ uzly (pro na´sledujıćı´ vy´pocˇty). Pro dany´ uzel A definujeme ba´zovou funkci ϕA , jejı´zˇ funkcˇnı´ hodnoty si mu˚zˇeme prˇedstavit jako pla´sˇt’sˇestiboke´ho jehlanu nad troju´helnı´ky obsahujıćı´mi vrchol 27

2. METODY RˇESˇENI´ MODELU A, ve ktere´m je hodnota ϕA rovna jedne´. (Mimo uvedene´ troju´helnı´ky je ϕA nulova´ – viz obra´zek 2.5.) Funkce ϕi , i = 1, . . . , n jsou tedy po cˇa´stech linea´rnı´ a jejich derivace jsou konstantnı´. Tato volba umozˇnˇuje snadny´ vy´pocˇet koeficientu˚ (Aϕi , ϕj ). Soustavu (2.1) pak rˇesˇ´ıme vhodnou iteracˇnı´ metodou. Hledana´ aproximace vn funkce v pak bude po cˇa´stech linea´rnı´ a bude urcˇena hodnotami v jednotlivyćh uzlech sı´teˇ. A

A

Obra´zek 2.5: Ba´zova´ funkce ϕA .

B

hranice

B

Obra´zek 2.6: Hranicˇnı´ ba´zova´ funkce ϕB .

Ma´me-li k dispozici dostatecˇne´ mnozˇstvı´ dat, mu˚zˇeme kromeˇ vrcholu˚ troju´helnı´ku˚ sı´teˇ pouzˇ´ıt pro vy´pocˇet i dalsˇ´ı body, naprˇ´ıklad strˇedy stran troju´helnı´ku˚, jejich teˇzˇisˇteˇ a podobneˇ. Funkce ϕi konstruujeme analogicky, pouze zvy´sˇ´ıme stupenˇ polynomu˚, ze kteryćh se funkce skla´dajı´, tedy ϕi budou po cˇa´stech kvadraticke´, kubicke´ atd., podle pocˇtu pouzˇityćh bodu˚. Vy´sledna´ funkce vn pak bude hladsˇ´ı a prˇesneˇjsˇ´ı. V prˇ´ıpadeˇ nehomogennı´ okrajove´ podmıńky volı´me navıć ba´zove´ funkce ϕi take´ pro krajnı´ body, podobneˇ jako pro vnitrˇnı´ uzly. Ty budou nenulove´ pouze nad prˇ´ıslusˇny´m pocˇtem hranicˇnıćh troju´helnı´ku˚ (viz obra´zek 2.6). Integra´ly teˇchto funkcı´ prˇes hranici oblasti Ω pak budou rovny dany´m okrajovy´m podmıńka´m.

2.3 Srovnańı´ metod Prˇestozˇe princip obou metod je pomeˇrneˇ odlisˇny´, majı´ metoda konecˇnyćh diferencı´ a metoda konecˇnyćh prvku˚ mnoho spolecˇne´ho. (Metodu konecˇnyćh diferencı´ lze dokonce prˇi vhodne´ volbeˇ ba´zovyćh funkcı´ povazˇovat za specia´lnı´ prˇ´ıpad metody konecˇnyćh prvku˚.) Obeˇ metody vedou k rˇesˇenı´ soustavy s pa´sovou maticı´, prˇ´ıpadneˇ s velmi rˇ´ıdkou maticı´. Metoda konecˇnyćh diferencı´ aproximuje zadanou diferencia´lnı´ rovnici, zatı´mco metoda konecˇnyćh prvku˚ aproximuje jejı´ prˇesne´ rˇesˇenı´. Vy´hodou metody konecˇnyćh diferencı´ je jejı´ velmi jednoducha´ mysˇlenka a snadna´ implementace. Jejı´ nevy´hodou oproti metodeˇ konecˇnyćh prvku˚ je ovsˇem znacˇneˇ mala´ geometricka´ prˇizpu˚sobivost. Sı´t’ musı´ by´t zcela pravidelna´ a jejı´ vystizˇenı´ hranice oblasti je velmi hrube´. Oproti tomu metoda konecˇnyćh prvku˚ umozˇnˇuje nejen lepsˇ´ı vystizˇenı´ hranice, ale zejmeńa prˇipousˇtı´ nepravidelnost sı´teˇ. Prˇi modelovańı´ proudeˇnı´ podzemnı´ vody je cˇasto trˇeba do modelu zahrnout velke´ oblasti, abychom zı´skali uceleny´ obraz modelovane´ situace. Za´rovenˇ je prˇitom vhodne´ neˇktere´ cˇa´sti, naprˇ´ıklad kde se vyskytuje kontaminace, modelovat 28

2. METODY RˇESˇENI´ MODELU prˇesneˇji. V neposlednı´ rˇadeˇ je metoda konecˇnyćh prvku˚ znacˇneˇ prˇizpu˚sobiva´, co se tyćˇe konkre´tnı´ polohy uzlu˚. Sı´t’lze tedy snadno prˇizpu˚sobit poloze vrtu˚, studnı´ a podobneˇ. Zatı´mco drˇ´ıve nebyla metoda konecˇnyćh prvku˚ pro proble´m modelovańı´ pohybu kapalin vu˚bec vyuzˇ´ıvańa, vy´sˇe uvedene´ du˚vody zpu˚sobily, zˇe se v poslednı´ dobeˇ stala rovnocennou metodou k metodeˇ konecˇnyćh diferencı´. Neexistuje ovsˇem zˇa´dny´ algoritmus, ktery´ by na´m umozˇnil rozhodnout, ktera´ z metod bude vhodneˇjsˇ´ı pro rˇesˇenı´ dane´ho proble´mu.

29

Kapitola 3

Model oblasti Slavkov u Brna V te´to kapitole uka´zˇeme pouzˇitı´ matematicke´ho modelu na oblasti Slavkov u Brna. Data pouzˇita´ v te´to kapitole pocha´zı´ ze zpra´vy prˇedsanacˇnı´ho dopru˚zkumu spolecˇnosti GEOtest Brno [GEOtest]. Za´jmova´ oblast modelu o rozloze asi 0,3 km2 se nacha´zı´ v jizˇnı´ cˇa´sti meˇsta Slavkov u Brna v okolı´ zˇeleznicˇnı´ stanice (viz obra´zek 3.3). V te´to oblasti se nacha´zı´ area´l spolecˇnosti EMP, s.r.o., ktery´ drˇ´ıve patrˇil podniku MEZ Brno. Cˇinnostı´ spolecˇnosti MEZ byla vy´roba, monta´zˇ a u´drzˇba elektrotechnickyćh zarˇ´ızenı´. Prˇitom docha´zelo k nezodpoveˇdne´mu zacha´zenı´ s nebezpecˇny´mi la´tkami a to zejmeńa s nepola´rnı´mi extrahovatelny´mi la´tkami (chlorovane´ uhlovodı´ky, ropne´ uhlovodı´ky, atd.), ktere´ se vyuzˇ´ıvaly jako rozpousˇteˇdla a odmasˇt’ovadla. Du˚sledkem toho vznikla v okolı´ stara´ ekologicka´ za´teˇzˇ. Spolecˇnost GEOtest provedla v te´to oblasti prˇedsanacˇnı´ pru˚zkum, jehozˇ soucˇa´stı´ byl odhad u´rovneˇ znecˇisˇteˇnı´ podzemnı´ vody pomocı´ modelu proudeˇnı´. V ra´mci tohoto pru˚zkumu byly provedeny vrtne´ a sonda´zˇnı´ praće, odbeˇr vzorku˚ zemin, stavebnıćh konstrukcı´, podzemnı´ a povrchove´ vody, hydrogeologicky´ pru˚zkum a geodeticke´ zameˇrˇenı´. V prvnı´ fa´zi pru˚zkumu byla lokalizovańa mı´sta, kde se zacha´zelo s nebezpecˇny´mi la´tkami, jako prima´rnı´ zdroje kontaminace. Na´sledneˇ se prˇedbeˇzˇneˇ odhadl postup kontaminace. Nejdu˚lezˇiteˇjsˇ´ım faktorem prˇitom bylo stanovenı´ privilegovanyćh cest proudeˇnı´ podzemnı´ vody (viz obra´zek 3.4). Privilegovane´ cesty proudeˇnı´ se nacha´zı´ naprˇ´ıklad pode´l tektonickyćh zlomu˚, pode´l nichzˇ voda proudı´ rychleji, a mohou slouzˇit jako koridor pro cestu kontaminace. Urcˇujı´ se morfohydrogeometrickou analy´zou ze spa´du tereńu, geologicke´ stavby a tektoniky oblasti. Nalezenı´ privilegovanyćh cest je pomeˇrneˇ obtı´zˇne´ a vyzˇaduje odborne´ geologicke´ posouzenı´. Jejich pru˚beˇh se prˇitom v modelu proudeˇnı´ a v mapaćh izohydrohyps nemusı´ vu˚bec projevit. S prˇihle´dnutı´m ke stanovenı´ teˇchto cest a odhadovany´m rozsˇ´ırˇenı´m kontaminace se na´sledneˇ urcˇil pocˇet a umı´steˇnı´ vrtu˚ a sond pru˚zkumu. Prˇi jejich plańovańı´ bylo trˇeba bra´t zrˇetel take´ na inzˇeny´rske´ sı´teˇ v oblasti a zejmeńa na trat’ zˇeleznice. Celkem bylo uskutecˇneˇno sˇest novyćh hydrogeologickyćh vrtu˚ a 43 sond. Byl proveden odbeˇr vzorku˚ a hydrodynamicke´ zkousˇky pro urcˇenı´ parametru˚ modelu. V pru˚beˇhu vrtańı´ a hydrodynamickyćh zkousˇek byly meˇrˇeny na´sledujıćı´ u´daje: •

poloha nepropustne´ho podlozˇ´ı a stropu

•

vy´sˇka hladiny podzemnı´ vody

•

u´daje pro vy´pocˇet koeficientu˚ filtrace. 31

3. MODEL OBLASTI SLAVKOV U BRNA Na´sledneˇ byl vytvorˇen popis zkoumane´ oblasti, na jehozˇ za´kladeˇ byl sestaven koncepcˇnı´ model.

Obra´zek 3.1: Odebrane´ vzorky zemin. Za´jmova´ oblast se deˇlı´ na dveˇ zcela ru˚znorode´ cˇa´sti – jizˇnı´ a severnı´ (viz obra´zek 3.3). Jizˇnı´ cˇa´st je z veˇtsˇiny zastaveˇna´, a proto zde te´meˇrˇ nedocha´zı´ ke vsakovańı´ sra´zˇek a dotaci podzemnı´ vody. Nacha´zı´ se zde zvodenˇ s volnou hladinou, jejı´zˇ tlousˇt’ka se pohybuje kolem peˇti metru˚. Celkova´ propustnost v te´to oblasti je velmi slaba´ (rˇa´doveˇ 10−7 azˇ 10−8 ms−1 ) s obcˇasny´mi pı´scˇity´mi vlozˇkami, ktere´ ji zvysˇujı´. Dı´ky volne´ hladineˇ a pomeˇrneˇ male´mu mnozˇstvı´ vody zde nenı´ trˇeba bra´t v u´vahu nepropustny´ strop zvodneˇ. Naopak severnı´ cˇa´st oblasti je tvorˇena prˇeva´zˇneˇ pı´sky a sˇteˇrky s velmi silnou propustnostı´ (rˇa´doveˇ 10−4 ms−1 ), ktere´ jsou shora ohranicˇene´ povodnˇovy´mi hlıńami. Ty tvorˇ´ı hornı´ nepropustnou vrstvu zvodneˇ, dı´ky cˇemuzˇ se zde vyvinula napjata´ zvodenˇ. Sra´zˇky zde rovneˇzˇ nezapocˇ´ıta´va´me, pra´veˇ dı´ky hornı´ nepropustne´ vrstveˇ, ktera´ zabranˇuje jejich pru˚saku do zvodneˇ. Prˇiblı´zˇenı´ situace lze videˇt v geologicke´m profilu u´zemı´ na obra´zku 3.6. Z nameˇrˇenyćh vy´sˇek podzemnı´ vody a nepropustnyćh podlozˇ´ı v jednotlivyćh vrtech byla programem SURFER interpolovańa hladina podzemnı´ vody (viz obra´zek 3.5) a vy´sˇka nepropustne´ho dna a stropu zvodneˇ. Program SURFER umozˇnˇuje modelovat trojrozmeˇrne´ povrchy vyuzˇitı´m ru˚znyćh interpolacˇnıćh metod. V tomto prˇ´ıpadeˇ bylo vyuzˇito metody kriging, ktera´ vycha´zı´ z metody nejmensˇ´ıch cˇtvercu˚ [WikFE]. Tato metoda umozˇnˇuje interpolovat hodnoty nezna´me´ funkce na pravidelne´ sı´ti bodu˚, ktera´ odpovı´da´ modelu, zatı´mco zna´me´ hodnoty jsou obvykle rozmı´steˇny zcela nesystematicky.

32

3. MODEL OBLASTI SLAVKOV U BRNA

Obra´zek 3.2: Vrtna´ souprava.

33


3.1 Koncepcˇnı´ model Cı´lem modelu bylo simulovat proudeˇnı´ podzemnı´ vody za beˇzˇnyćh podmıńek (tedy bez cˇerpańı´ cˇi ru˚znyćh sanacˇnıćh za´sahu˚), aby bylo mozˇno urcˇit prˇiblizˇny´ rozsah a postup kontaminace. Vzhledem k pomeˇrneˇ male´ vy´sˇce zvodneˇ je vertika´lnı´ proudeˇnı´ sice velmi slabe´, ale z hlediska kontaminace je nezanedbatelne´. Hodnoty koeficientu˚ filtrace byly stanoveny na za´kladeˇ cˇerpacıćh zkousˇek v neˇkteryćh vrtech. Jejich hodnoty byly voleny konstantnı´ na jednotlivyćh oblastech a stejne´ pro vsˇechny vodorovne´ smeˇry. Pro vertika´lnı´ smeˇr byly z du˚vodu male´ komunikace voleny o rˇa´d nizˇsˇ´ı. Celkem byl zvolen staciona´rnı´ trˇ´ırozmeˇrny´ model proudeˇnı´ nestlacˇitelny´m anizotropnı´m prostrˇedı´m. Jako popisujıćı´ rovnice modelu byla tedy vybrańa rovnice (1.13)(a). Hladina podzemnı´ vody v za´jmove´ oblasti je smı´sˇena´ (prˇechod mezi volnou a napjatou hladinou). S ohledem na tuto skutecˇnost byla pro rˇesˇenı´ modelu zvolena metoda konecˇnyćh diferencı´. Byla volena pravidelna´ cˇtvercova´ sı´t’modelu. Vzhledem k prˇesnosti geodeticke´ho zameˇrˇenı´ a rozloze zkoumane´ oblasti byl zvolen krok sı´teˇ deset metru˚. Ve zkoumane´ oblasti nedocha´zı´ k odbeˇru˚m podzemnı´ vody. V jizˇnı´ cˇa´sti je vsakovańı´ sra´zˇek velmi slabe´ nebo zˇa´dne´. V severnı´ cˇa´sti znemozˇnˇuje vsakovańı´ sra´zˇek hornı´ nepropustna´ vrstva zvodneˇ. Prava´ strana rovnice (ktera´ vyjadrˇuje zdroje a propady vody v oblasti) je tedy rovna nule. Zdrojem vody v oblasti je pouze prˇ´ıtok podzemnı´ vody prˇes jizˇnı´ a vyćhodnı´ hranici modelu. Okrajove´ podmıńky modelu byly voleny umeˇle na za´kladeˇ mapy nameˇrˇenyćh izohydrohyps (viz obra´zek 3.5). Byla volena podmıńka prvnı´ho typu, tedy vy´sˇka hladiny na hranici modelu. Hranice modelu byly voleny podle tvaru za´jmove´ oblasti a polohy izohydrohyps. Vy´sledna´ sı´t’modelu je zobrazena na obra´zku 3.7.

ˇ esˇenı´ a kalibrace modelu 3.2 R Pro rˇesˇenı´ modelu byl zvolen program Visual MODFLOW Pro. Program MODFLOW je zdarma dostupny´ software, ktery´ vyvinula americka´ neza´visla´ veˇdeckovy´zkumna´ vla´dnı´ agentura U.S. Geological Survey. MODFLOW je jednı´m z nejpouzˇ´ıvaneˇjsˇ´ıch programu˚ pro vy´pocˇet rovnice proudeˇnı´ podzemnı´ vody, ktery´ prˇeva´dı´ parcia´lnı´ diferencia´lnı´ rovnici metodou konecˇnyćh diferencı´ na soustavu rovnic. Komercˇnı´ nadstavba Visual MODFLOW Pro da´le obsahuje numericke´ metody pro rˇesˇenı´ soustav z programu MODFLOW a umozˇnˇuje vizualizaci spocˇ´ıtane´ hladiny podzemnı´ vody (vykreslenı´ izohydrohyps, vektoru˚ proudeˇnı´ apod.). V prvnı´ fa´zi byly do modelu zadańy okrajove´ podmıńky a hodnoty nepropustnyćh podlozˇ´ı v jednotlivyćh bodech sı´teˇ z programu SURFER. Koeficienty filtrace byly zadańy konstantnı´mi hodnotami podle nameˇrˇenyćh u´daju˚ (viz obra´zek 3.8). Pote´ byla vypocˇtena hladina podzemnı´ vody s teˇmito hodnotami a vygenerovańo jejı´ zna´zorneˇnı´ (viz obra´zek 3.9). Srovnańı´m nameˇrˇenyćh izohydrohyps (viz obra´zek 3.5) a namodelovane´ho stavu vidı´me, zˇe namodelovany´ stav neodpovı´da´ skutecˇnosti. To je dańo neprˇesny´mi hodnotami koeficientu˚ propustnosti. Bylo tedy nutne´ model opravit, aby modelovany´ stav 34

3. MODEL OBLASTI SLAVKOV U BRNA le´pe vystihoval nameˇrˇene´ hodnoty. Tento proces nazy´va´me kalibrace modelu. Postupneˇ byly opravovańy jednotlive´ hodnoty koeficientu˚ filtrace a model byl prˇepocˇ´ıta´vań. Po neˇkolika desı´tkaćh prˇiblı´zˇenı´ byl model nastaven do stavu, ktery´ se pomeˇrneˇ dobrˇe shoduje s nameˇrˇeny´mi izohydrohypsami (viz obra´zek 3.10, konecˇne´ hodnoty koeficientu˚ filtrace vidı´me na obra´zku 3.11). Vizualizaci namodelovane´ hladiny podzemnı´ vody ukazujı´ obra´zky 3.13, 3.14 a 3.15. Pro rˇesˇenı´ syste´mu rovnic vygenerovanyćh programem MODFLOW byl vyuzˇit solver WHS. Ten vyuzˇ´ıva´ Stoneovy metody neu´plne´ho rozkladu matice soustavy na soucˇin dolnı´ a hornı´ troju´helnı´kove´ matice, ktera´ je velmi oblı´bena´ pro rˇesˇenı´ syste´mu˚, ktere´ vznikly z parcia´lnıćh diferencia´lnıćh rovnic metodou konecˇnyćh diferencı´. Pro soustavu Ax = b hleda´me rozklad matice A ve tvaru A = LU − M, kde L je dolnı´ a U hornı´ troju´helnı´kova´ matice a M je matice, pro kterou platı´ kLU k >> kM k v neˇjake´ maticove´ normeˇ [WikFE]. Rozklad matice A prˇitom hleda´me tak, aby matice LU byla le´pe podmıńeˇna´ nezˇ matice A, cozˇ zvysˇuje stabilitu rˇesˇenı´. Program WHS rˇesˇ´ı soustavu ve dvou krocıćh. Prvnı´ krok hleda´ optima´lnı´ rozklad matice A (vneˇjsˇ´ı iterace). Druhy´ krok rˇesˇ´ı danou soustavu specia´lnı´ iteracˇnı´ metodou s prˇ´ıslusˇny´m rozkladem matice soustavy (vnitrˇnı´ iterace) [VM Pro]. Jako krite´rium pro zastavenı´ vy´pocˇtu byly pouzˇity na´sledujıćı´ rozdı´ly: •

rozdı´l vy´sˇek hladin podzemnı´ vody v kazˇde´m bodeˇ sı´teˇ mezi jednotlivy´mi iteracemi je mensˇ´ı nezˇ 10−4 a za´rovenˇ

•

rozdı´l reziduı´ pro jednotlive´ iterace je mensˇ´ı nezˇ 10−4 . (Oznacˇ´ıme-li xk hodnotu k-te´ iterace rˇesˇenı´, pak rk = b − Axk je k-te´ reziduum.)

Pro dobrou konvergenci metody je take´ vhodne´ zvolit dobrou pocˇa´tecˇnı´ aproximaci. (Naprˇ. pouze konstantnı´ hodnota pocˇa´tecˇnı´ aproximace mu˚zˇe zpu˚sobit divergenci metody.) V nasˇem prˇ´ıpadeˇ byla jako pocˇa´tecˇnı´ aproximace zvolena hodnota nameˇrˇene´ (resp. interpolovane´) hladiny podzemnı´ vody. Zvolena´ metoda velmi dobrˇe a rychle konvergovala. (Pru˚beˇh vy´pocˇtu a dalsˇ´ı parametry ukazuje obra´zek 3.12.) Za´veˇrem lze rˇ´ıci, zˇe vybrany´ model a numericke´ metody splnily ocˇeka´vańı´. Namodelovany´ stav proudeˇnı´ podzemnı´ vody dobrˇe odpovı´dal rea´lne´ situaci a mohl by´t vyuzˇit pro simulaci postupu kontaminace (viz obra´zek 3.16). Na za´kladeˇ tohoto modelu bude v budoucnu naplańovań sanacˇnı´ za´sah, jehozˇ cı´lem bude odstraneˇnı´ kontaminace a uvedenı´ prostrˇedı´ do pu˚vodnı´ho stavu.

35


Obra´zek 3.3: Mapa za´jmove´ oblasti s rozdeˇlenı´m jizˇnı´ a severnı´ cˇa´sti. 36


Obra´zek 3.4: Privilegovane´ cesty proudeˇnı´ podzemnı´ vody. 37


Obra´zek 3.5: Mapa nameˇrˇenyćh izohydrohyps. 38


Obra´zek 3.6: Geologicky´ rˇez za´jmovou oblastı´.

39


Obra´zek 3.7: Sı´t’modelu s vyznacˇenı´m okrajovyćh podmıńek. 40


Obra´zek 3.8: Prvnı´ odhad koeficientu˚ propustnosti.

41


Obra´zek 3.9: Prvnı´ modelovane´ izohydrohypsy.

42


Obra´zek 3.10: Konecˇny´ stav modelovanyćh izohydrohyps. 43


Obra´zek 3.11: Konecˇne´ nastavenı´ koeficientu˚ propustnosti.

44


Obra´zek 3.12: Uka´zka pru˚beˇhu vy´pocˇtu.

45


Obra´zek 3.13: Vizualizace modelovane´ hladiny podzemnı´ vody.

Obra´zek 3.14: Vizualizace modelovane´ hladiny podzemnı´ vody, nepropustne´ho podlozˇ´ı a pozice vrtu˚. 46


Obra´zek 3.15: Vizualizace modelovane´ hladiny podzemnı´ vody a zemske´ho povrchu. Tam, kde voda prˇesahuje povrch, je hladina napjata´.

47


Obra´zek 3.16: Modelovany´ rozsah kontaminace. 48

Literatura ˇ ´ıha] [R

ˇ ´ıha a kolektiv. Matematicke´ modelovańı´ hydrodynamickyćh a disperzIng. Jaromı´r R nıćh jevu˚. Vysoke´ ucˇenı´ technicke´ v Brneˇ, 1997.

[Vita´sek] RNDr. Emil Vita´sek, CSc. Numericke´ metody. Nakladatelstvı´ technicke´ literatury, Praha, 1987. [Rektorys] prof. RNDr. Karel Rektorys, DrSc. Variacˇnı´ metody v inzˇeny´rskyćh proble´mech a v proble´mech matematicke´ fyziky. Academia, Praha, 1999. [Kazda] Ing. Ivo Kazda, DrSc. Podzemnı´ hydraulika v ekologickyćh a inzˇeny´rskyćh aplikacıćh. Academia, Praha, 1997. [GEOtest] GEOtest Brno, a.s. EMP s.r.o. – prˇedsanacˇnı´ dopru˚zkum a projektova´ dokumentace sanacˇnı´ho za´sahu. Za´veˇrecˇna´ zpra´va prˇedsanacˇnı´ho dopru˚zkumu, Brno, 2009. [WikOE] Wikipedie, otevrˇena´ encyklopedie. http://cs.wikipedia.org/wiki/Hlavn%C3%AD strana.

[WikFE] Wikipedia, the free encyclopedia. http://en.wikipedia.org/wiki/Main Page.

[VM Pro] Waterloo Hydrogeologic Inc. Visual MODFLOW Pro user’s manual. Waterloo Hydrogeologic Inc., Waterloo, Ontario, Kanada, 2003. [HoNo] Milan Hokr a Jan Nova´k. Transportnı´ procesy v pore´znı´m prostrˇedı´ - fyzika´lnı´ popis. www.fm.vslib.cz/kmo/czech/vyuka/trp/ \ hokr novak transportni procesy fyzika 03-01-15.pdf.

[Matousˇek] Dr. ing. Vaćlav Matousˇek. Proudeˇnı´ podzemnı´ vody. http://hydraulika.fsv.cvut.cz/users/matousek/downloads/ \ 08 Podzemni voda vm.pdf.

[La´gner] Antonıń La´gner. Kolobeˇh vody v prˇ´ırodeˇ - minera´lnı´ a podzemnı´ vody. http://www.priroda.cz.

[PoVo]

Podzemnı´ voda. http://geotech.fce.vutbr.cz/wwwroot/scripta/ \ geologie/podvoda.htm.

49

3. MODEL OBLASTI SLAVKOV U BRNA [CHU]

Cˇesky´ hydrometeorologicky´ u´stav. Obeˇh vody. http://ga.water.usgs.gov/edu/watercycleczech.html

[EPA]

Jacob Bear, Milovan S. Beljin, Randall R. Ross. Fundamentals of Ground-Water Modeling. http://www.epa.gov/tio/tsp/download/issue13.pdf.

50

MASARYKOVA UNIVERZITA V BRNĚ. vody. Mgr. Eva Cvachovcová

Recommend Documents