Kombinovaný model proud ní podzemní vody založený na primární formulaci MKP

TECHNICKÁ UNIVERZITA V LIBERCI

Fakulta mechatroniky a mezioborových inženýrských studií

Studijní program: N2612 – Elektrotechnika a informatika Studijní obor: 3901T025 – P írodov dné inženýrství

Kombinovaný model proud ní podzemní vody založený na primární formulaci MKP

Combined model of the groundwater flow based on primar formulation of FEM

Diplomová práce

Autor: Bc. Jan Palek

Vedoucí práce: Ing. Otto Severýn, Ph.D.

Konzultant: Ing. Ji ina Královcová, Ph.D.

V Liberci 16. 5. 2008

Prohlášení

Byl jsem seznámen s tím, že na mou diplomovou 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 TUL má právo na uzav ení licen ní smlouvy o užití mé DP a prohlašuji, že s o u h l a s í m s p ípadným užitím mé diplomové práce (prodej, zap j ení apod.).

Jsem si v dom toho, že užít své diplomové práce i poskytnout licenci k jejímu využití mohu jen se souhlasem TUL, která má právo ode mne požadovat p im ený p ísp vek na úhradu náklad , vynaložených univerzitou na vytvo ení díla (až do jejich skute né výše).

Diplomovou práci jsem vypracoval samostatn s použitím uvedené literatury a na základ konzultací s vedoucím diplomové práce a konzultantem.

Datum:

Podpis:

Pod kování

Rád bych pod koval vedoucímu mé diplomové práce panu Ing. Ottu Severýnovi, Ph.D. za jeho odbornou pomoc, cenné rady a poskytnuté informace.

Dále bych rád pod koval panu Prof. Dr. Ing. Ji ímu Maryškovi, CSc. za jeho odbornou pomoc týkající se matematické ásti této práce a paní Ing. Ji in Královcové, Ph.D. za poskytnuté sít na testování.

Anotace

Tato diplomová práce se zabývá vývojem modelu proud ní podzemní vody, který je realizován programem Flow123D a je založen na primární formulaci MKP.

Bylo nutné navrhnout a implementovat algoritmus, který zajistí správnou funk nost modelu pro multidimenzionální sít .

Po seznámení s použitým fyzikálním a matematickým modelem uvedeným v první kapitole se v nujeme popisu návrhu a implementace zmín ného modelu do programu Flow123D. Implementovaný algoritmus je v záv ru práce otestován na benchmarkových úlohách.

Program Flow123D, testovací úlohy i software pro vykreslení výsledk modelu je možné nalézt na p iložením CD.

Klí ová slova: primární formulace, MKP, proud ní podzemní vody, multidimenzionální p ístup

Abstract

The aim of this Diploma Thesis is to develop a model of the groundwater flow.

The model of this physical process is realized by program Flow123D and it is based on primar formulation of the FEM It was necessary to create and to implement an algorithm for correct calculation in case of multidimensional meshes. Used physical and mathematical model is briefly explained in the first chapter. The second chapter gives a concept and description of implementation of the new algorithm into the Flow123D program. There is tested the correct functionality of the implemented algorithm on several benchmark problems in the third chapter.

Flow123D program, input files, results of testing problems and visualization software is possible to find on enclosed CD.

Keywords: primar formulation, FEM, groundwater flow, multidimensional approach

Obsah

OBSAH ... - 6 -

SEZNAM OBRÁZK ... - 7 -

ÚVOD ... - 8 -

1 MODEL PROUD NÍ PODZEMNÍCH VOD ... - 11 -

1.1 Fyzikální popis úlohy ... - 11 -

1.2 Slabá formulace ešení, aplikace MKP ... - 14 -

1.3 Po íta ová implementace modelu ... - 19 -

2 NÁVRH A IMPLEMENTACE ALGORITMU ... - 21 -

2.1 Návrh algoritmu pro kombinaci 2D a 1D element ... - 23 -

2.2 Návrh algoritmu pro kombinaci 3D a 2D element ... - 26 -

3 IMPLEMENTACE ALGORITMU ... - 29 -

3.1 Použití primárního modelu ... - 29 -

3.2 Algoritmus z hlediska programátora ... - 30 -

3.2.1 Multinode_test ... - 32 -

3.2.2 Add_nodes ... - 32 -

3.2.3 Create_elements ... - 33 -

3.2.4 Assign_multinodes ... - 33 -

3.2.5 Alocate_new_nodes ... - 33 -

3.2.6 Fill_new_nodes ... - 33 -

3.2.7 Area2D ... - 34 -

3.2.8 Area3D ... - 35 -

3.2.9 Create_msh_file ... - 36 -

3.2.10 Create_bcd_file ... - 36 -

3.2.11 Topology ... - 36 -

3.2.12 Multi_n_row ... - 37 -

3.2.13 Multi_fill_col_num ... - 37 -

3.2.14 Multi_values ... - 38 -

4 TESTOVACÍ ÚLOHY ... - 39 -

4.1 Test z úvodu ... - 39 -

4.2 Test 2D-1D na rozsáhlejší oblasti ... - 40 -

4.3 Test 2D-3D na rozsáhlejší oblasti ... - 42 -

5 ZÁV R ... - 51 -

Seznam obrázk

Obr. 1-1, sí ... - 9 -

Obr. 1-2, výsledek kombinace 1D a 2D ... - 9 -

Obr. 1-3 - kompatibilní spojení 2D a 1D ... - 20 -

Obr. 1-4 - nekompatibilní spojení 2D a 1D ... - 20 -

Obr. 1-5 - kompatibilní spojení 3D a 2D ... - 20 -

Obr. 1-6 - nekompatibilní spojení 3D a 2D ... - 20 -

Obr. 2-1, sestavení globální matice ... - 21 -

Obr. 2-2, spole né uzly pro 1D a 2D elementy ... - 22 -

Obr. 2-3, princip odd lení element 1D a 2D ... - 23 -

Obr. 2-4, znázorn ní p idávání krajních uzl 1D elementu ... - 25 -

Obr. 2-5, požadovaná sí ... - 25 -

Obr. 2-6, princip odd lení element 2D a 3D ... - 26 -

Obr. 3-1, vývojový diagram ... - 29 -

Obr. 3-2, 1. fáze, vývojový diagram ... - 31 -

Obr. 3-3, Area2D ... - 34 -

Obr. 3-4, Area3D ... - 35 -

Obr. 3-5, 2.fáze, vývojový diagram ... - 37 -

Obr. 4-2, sí ... - 39 -

Obr. 4-3, sí ... - 40 -

Obr. 4-4, = 0.01 ... - 41 -

Obr. 4-5, = 1.0 ... - 41 -

Obr. 4-6, = 10.0 ... - 42 -

Obr. 4-7, sí ... - 43 -

Obr. 4-8, puklina v síti ... - 43 -

Obr. 4-9, = 0.001,Primar ... - 44 -

Obr. 4-10, = 0.001,MH ... - 44 -

Obr. 4-11, = 1.0,Primar ... - 45 -

Obr. 4-12, = 1.0,MH ... - 45 -

Obr. 4-13, = 10.0,Primar ... - 46 -

Obr. 4-14, = 10.0,MH ... - 46 -

Obr. 4-15, = 0.001,Primar ... - 47 -

Obr. 4-16, = 0.001,MH ... - 47 -

Obr. 4-17, = 1.0,Primar ... - 48 -

Obr. 4-18, = 1.0,MH ... - 48 -

Obr. 4-19, = 10.0,Primar ... - 49 -

Obr. 4-20, = 10.0,MH ... - 49 -

ÚVOD

Tématem této diplomové práce je vývoj modelu proud ní podzemních vod.

Model je založen na primární formulaci metody kone ných prvk (MKP). Modelování tohoto p írodního jevu je základem pro ešení hydrogeologických, geochemických a hlavn ekologických problém .

Pro tyto ú ely byl na kated e modelování proces fakulty Mechatroniky a mezioborových inženýrských studií a dodnes stále je vyvíjen program Flow123D, více p edstaven v dalším textu, nebo v [2]. P vodn byl výpo et aproximovaného ešení ustáleného filtra ního proud ní založen pouze na smíšené hybridní (MH) formulaci MKP, více v [5]. Do Flow123D byly v rámci této formulace implementovány možnosti použití 1D element (linie), 2D element (trojúhelníky) a 3D element ( ty st ny) i jejich kombinace, tzn. umí pracovat s oblastmi, ve kterých se pro 3D prost edí nacházejí krom porézního média aproximovaného 3D elementy i pukliny aproximované 2D elementy a pro 2D prost edí je porézní médium aproximováno 2D elementem a puklina 1D elementem. Hovo íme o tzv. multidimenzionálním p ístupu a multidimenzionálních sítích.

Pozd ji byla do programu z d vod nevýhod (MH) formulace (mj. vede na velké indefinitní matice) implementována i primární formulace MKP, podrobn ji v [6].

I v rámci této formulace byly do programu implementovány stejné t i prvky zastupující každou dimenzi. Funk nost této implementace byla ov ena na n kolika základních testovacích úlohách. Testy byly provád ny nejprve na t ch úlohách, jejichž sít obsahovaly pouze elementy jedné dimenze, výsledky t chto aproximací byly ve všech t ech p ípadech shodné s analytickými výsledky. Druhým typem test byly testy na multidimenzionálních sítích. V n kterých situacích byly výsledky korektní, v jiných nikoliv. Konkrétn v p ípad p í ného p etoku p es pukliny – a v kombinaci 3D a 2D nebo v kombinaci 2D a 1D element – již docházelo k rozpor m numerických výsledk s výsledky p edpokládanými. Tok p es puklinu byl nezávislý na propustnosti pukliny, což neodpovídalo reálnému chování.

Pro názornost zde uvedeme elementární p íklad, na kterém je z eteln vid t výše zmín ný nedostatek.

M jme sí na obr.1-1, Dirichletovy okrajové podmínky v uzlech P6 = 10 m, P2 = 1 m – t mi docílíme toku z levého horního rohu do pravého dolního rohu sít ,

propustnost na 2D elementech KT = 1 0

0

1 m/den, na 1D elementech (na obr. 1-1 jsou

znázorn ny erven ) postupn KL = 10m/den, KL = 1 m/den a KL = 0.01 m/den.

Sí (elementy jsou pro názornost zmenšeny a uzly posunuty):

Obr. 1-1, sí

Výpo et tlakových polí a vektor rychlostí:

K = 10 K = 1 K = 0.01

Obr. 1-2, výsledek kombinace 1D a 2D

Z obr.1-2 je patrné, že výsledek je stále stejný, není zde žádná závislost na propustnosti pukliny. Intuitivn se dá o ekávat, že málo propustná puklina bude v tomto p ípad fungovat jako hydraulický izolátor obou ástí.

Tímto jsme se dostali k úkol m, které má tato práce vy ešit:

o Navrhnout algoritmus pro aproximaci úlohy filtra ního proud ní pomocí primární MKP a pro p idávání uzl sít na styku element r zných dimenzí a pro stanovení p etoku mezi takovýmito uzly

o Implementovat navržené algoritmy do programu Flow123D o Otestovat tyto algoritmy na modelových úlohách

Práce je roz len na do t ech kapitol. V první kapitole je uveden fyzikální popis problému a jeho aproximace pomocí primární formulace MKP a seznámení s programem Flow123D, do kterého byla provedena implementace našeho algoritmu.

Ve druhé kapitole je nejprve popsán zp sob, na kterém byla navržena komunikace element r zných dimenzí v primární formulaci p vodn [6] a dále navržen postup, kterým je dosaženo korektního chování modelu. V této kapitole je následn popsán algoritmus jak z hlediska uživatele programu Flow123D, tak i z hlediska programátora algoritmu.

Ve t etí kapitole jsou uvedeny testovací úlohy, na kterých je ukázána funk nost implementovaného algoritmu.

1 MODEL PROUD NÍ PODZEMNÍCH VOD

V této kapitole je zpracován základní matematicko-fyzikální popis filtra ního proud ní. Kapitola je rozd lena do t í ástí:

o Fyzikální popis úlohy

o Slabá formulace ešení, aplikace MKP o Po íta ová implementace modelu

1.1 Fyzikální popis úlohy

Filtra ní proud ní kapaliny budeme studovat v oblasti

Ω

, která je omezena hranicí

∂ Ω

. Fyzikální popis úlohy je dán dv ma zákony:

1. Darcyho zákonem:

∇ p

= K −

u

na

Ω

, (1)

kde

u [ ] ^ms

⁻¹ … vektor filtra ní rychlosti proud ní,

K [ ]

^{ms−1 … tenzor}

propustnosti a

∇ p

… gradient tlakové výšky. Uvažujeme lineární Darcyho zákon (

K

není funkcí p, ale pouze polohy ).

2. rovnicí kontinuity:

= q

∇ u ⋅

^na

Ω

, (2)

kde q

[ ]

^s−1 … hustota zdroj kapaliny

Rovnice kontinuity je vyjád ením zákona zachování hmoty respektive zákona zachování objemu pro nestla itelnou tekutinu.

K t mto rovnicím je dále nutné uvažovat okrajové podmínky. V našem modelu využíváme dva druhy okrajových podmínek:

o Dirichletova okrajová podmínka:

p

D

p =

na

∂ Ω

_D , (3) která p edepisuje hodnotu tlakové výšky na zadané ásti hranice.

o Neumannova okrajová podmínka:

u

N

⋅ n =

u

^na

∂ Ω

N , (4)

která p edepisuje hodnotu filtra ního toku na zbývající ásti hranice. Zde pD

[ ]

^m …hodnoty tlakové výšky na hranici oblasti,

u

_N

[ ]

^ms−1 …hodnoty toku p es hranici oblasti,

n

… jednotkový vektor vn jší normály k ivky tvo ící ást hranice

Ω

N

∂

.

Dirichletova okrajová podmínka modelované reálné oblasti je zadávána na liniích tok ek , na hladinách jezer, v m ených vrtech a podobn .

Neumannovu okrajovou podmínku obvykle zadáváme zejména pro vyjád ení nepropustné hranice odpovídající nulovému toku (homogenní Neumannova okrajová podmínka).

Pro úplnost zde ješt zmíníme nejobecn jší okrajovou podmínku, kterou však náš model nevyužívá. Jedná se o Newtonovu okrajovou podmínku, která vyjad uje stav, kdy vn jší tok hranicí závisí na spádu tlakové výšky a její tvar je:

) ( p − p

_D

=

⋅ n σ

u

na

∂ Ω

W , (5)

kde

σ [ ]

^s−1 je p estupový koeficient.

Filtra ní proud ní kapaliny p edstavuje pr tok kapaliny oblastí

Ω

, kde tato oblast je vypln na porézním médiem – strukturou složenou ze zrn nebo vláken pevné látky, mezi nimiž je volný prostor (póry), který m že být obsahovat vzduch nebo kapalinu. Dále m že oblast

Ω

obsahovat pukliny. Ty mohou být vypln ny málo propustným materiálem a tím tvo it hydraulický izolátor nebo naopak velmi dob e propoušt t kapalinu a tím být hydraulickým kolektorem. V našem p ípad je vždy celý pórový i puklinový prostor zapln n kapalinou a íkáme mu saturované prost edí. V této práci se zabýváme pouze stacionárním filtra ním proud ním.

To, jak je porézní prost edí pro kapalinu propustné i nepropustné, vyjad uje vlastnost – hydraulická vodivost zna ená

K

. Obecn se jedná o pozitivn definitní tenzor druhého ádu. Pokud je prost edí homogenní a navíc izotropní, lze vlastnost daného prost edí vyjád it jen jedním nezáporným reálným íslem a. Pro propustnost pak platí:

K

=

aE

, (6)

kde

E

je jednotková matice op t ádu shodného s dimenzí oblasti.

Cílem našeho modelu je získat hodnoty tlakových výšek

p [ ]

m v oblasti

Ω

. Tlaková výška je výška volné hladiny nad referen ním místem. Pro up esn ní uvedeme, že tlaková výška vychází z obecn jší piezometrické výšky

Φ [ ]

m , pro kterou byl Darcyho zákon p vodn odvozen. Její tvar je:

g h p

ρ +

∗

=

Φ

⁽⁷⁾

kde

h [ ]

m je referen ní výška,

p

^∗

[ ]

Pa je tlak, ρ

[

^{kg⋅ m−3}

]

je hustota kapaliny a g

[

^{m⋅ s−2}

]

je tíhové zrychlení. Výraz

g

p ρ

∗

[ ]

^m vyjad uje tlakovou výšku

p

. Darcyho zákon vyjád en pomocí piezometrické výšky je:

)

~ (

g h p ρ ∇ + ρ

^∗

−

= K g

u

⁽⁸⁾

kde

K ~

je tenzor propustnosti, hustota kapaliny je konstantní a stejn tak gravita ní zrychlení považujeme v celé oblasti

Ω

za konstantní. Pokud dále za referen ní výšku

h

dosazujeme nulovou výšku, lze Darcyho zákon vyjád it zp sobem jako je ve vztahu (1), kde je tenzor propustnosti

K

pouze korekcí pomocí hustoty a gravita ního zrychlení:

g K

K = ^~ ρ

⁽⁹⁾

Z vypo tených tlakových výšek je dále t eba vypo ítat i pole rychlosti filtra ního proud ní

u [ ] ^ms

⁻¹ v naší oblasti

Ω

. I když zde uvádíme jednotky rychlosti v

[ ] ^ms

⁻¹ , ve skute nosti se tato veli ina pohybuje ádov v

[

^{cm⋅ den−1}

]

^.

Výpo tem filtra ní rychlosti a tlakové výšky získáme pot ebné informace k posouzení hydrogeologické situace na oblasti

Ω

.

1.2 Slabá formulace ešení, aplikace MKP

V tomto odstavci si ukážeme, jak v p ípad naší konkrétní úlohy filtra ního proud ní dosp jeme od matematicko fyzikálního modelu uvedeného v p edchozí kapitole k modelu numerickému. Cílem této práce není p esné matematické odvození modelu, proto jsou uvedeny pouze hlavní kroky tohoto odvození.

Nejprve definujeme v dalším postupu použité prostory funkcí:

( ) ^Ω

L

2 ⁽¹⁰⁾

je prostor funkcí integrovatelných s druhou mocninou na oblasti

Ω

^.

( ) ^{Ω = ∈} ( ) ^Ω ( ^{+ ∇} ) ^{< ∞}

Ω

dx w w

L w

H

¹ ₂

;

^{2 2}

(11)

je Sobolev v prostor funkcí

w ^∈ L

2

( ) ^Ω

kvadraticky integrabilní v etn svého gradientu. Pro formulaci slabého ešení pot ebujeme dále prostor:

( ) ^{Ω =} { ^∈

¹

( ) ^{Ω ; = 0 ∈ ∂ Ω} }

1

0

w H w

H γ

⁽¹²⁾

kde

γ

je operátor stop vytvá ející restrikci na hranici oblasti. Dále pro tu ást hranice oblasti

∂ Ω

_D, kde je p edepsána Dirichletova okrajová podmínka, definujeme prostor:

( ) { ( )

^{D D}

}

D

w H w p

H

¹_,∗

Ω = ∈

¹

Ω ; γ = ∈ ∂ Ω

⁽¹³⁾

v tomto prostoru jsou spln ny Dirichletovy okrajové podmínky, zatímco:

( ) { ( )

^D

}

D

w H w

H

¹_,0

Ω = ∈

¹

Ω ; γ = 0 ∈ ∂ Ω

⁽¹⁴⁾

je prostor funkcí, které mají na hranici nulové stopy.

Dále si pro zjednodušení zápisu zavádíme skalární sou in funkcí v prostoru

( ) Ω L

2 :

( )

Ω

= u x v x dx v

u , ( ) ( )

⁽¹⁵⁾

skalární sou in funkcí v prostoru

L

2

( ) ^{∂ Ω}

:

Ω

∂

= u x v x dx v

u , ( ) ( )

⁽¹⁶⁾

Nyní se budeme v novat již samotnému postupu odvození modelu.

Pokud Darcyho zákon dosadíme do rovnice kontinuity, získáme tím rovnici v divergentním tvaru:

( ∇ ^p ) = ^q

⋅

∇

− K

⁽¹⁷⁾

Ob strany rovnice vynásobíme testovací funkcí

w ^∈ H

¹0

( ) ^Ω

a integrujeme p es oblast

Ω

, získáme rovnici ve tvaru:

( )

( ^{− ∇ ⋅ K ∇ p , w} ) ( )

_Ω

^{= q , w}

_Ω ⁽¹⁸⁾

Na vztah (18) aplikujeme Greenovu v tu a získáme:

( ∇ ∇ ) ( )

Ω

=

Ω

+

⋅

∇

− K p n , w K p , w q , w

_. ⁽¹⁹⁾

Kde len

K ∇ p ⋅ n

je zadaný tok

u

N na hranici oblasti, který p evedeme na pravou stranu:

( ^K ∇ ^p , ∇ ^w ) ( )

Ω

= ^q , ^w

Ω

+ ^u

N

, ^w

∂Ω_{N .} ⁽²⁰⁾

Definice:

Slabým ešením nazveme takovou funkci

^∈

1 ∗

( ) ^Ω

,

H

D

p

, která spl uje rovnici (20)

pro všechny testovací funkce

w ^∈ H

¹D,0

( ) ^Ω

_.

Nyní z d vodu ešitelnosti úlohy p evedeme nekone n dimenzionální problém na kone n dimenzionální použitím numerické metody – metody kone ných prvk . Aproximace ešení spo ívá v rozd lení spojité oblasti na kone ný po et prvk (element ), tyto elementy jsou jednozna n ur eny uzly (vrcholy element ). V primární formulaci hledáme hodnoty tlakových výšek na uzlech element .

Z matematického hlediska spo ívá aproximace ešení rovnice (20) v nahrazení hledané funkce

p

lineární kombinací bázových funkcí. Ze stejného d vodu je vhodné mít zadány v lineární kombinaci bázových funkcí i funkce zdroj

q

a funkce zadaných tok

u

:

∈

⋅

=

I i

i

P

i

p ϕ

_,

∈

⋅

=

I i

i

Q

i

q ϕ

_a

∈

⋅

=

IN i

i

G

i

u ϕ

⁽²¹⁾

kde

ϕ

i

^∈ H

h¹

( ) ^Ω

je bázová funkce a

H

¹h

( ) ^Ω

je podprostor

^H

¹

( ) ^Ω

s diskretiza ním parametrem h a

P

ⁱ,

Q

ⁱ a

G

ⁱ jsou hodnoty tlakových výšek, zdroj a tok , p i emž hodnoty tlakových výšek, které nejsou zadány okrajovou podmínkou, je nutné dopo ítat, ostatní hodnoty jsou z okrajových podmínek známé nebo v našem p ípad hodnoty

Q

ⁱ považovány za nulové (jako nulové zdroje v oblasti), ty je tedy možné z výpo tu vynechat.

Po nahrazení testovací funkce bázovou

w = ϕ

dosadíme vztahy (21) do rovnice (20) a dostaneme:

( ) ( ) ( ) (

i j

)

I i

i j

i I

i i j

i I

i i j

i I

I i

i

N D

D

G Q

K P K

P ϕ , ϕ ϕ , ϕ ϕ , ϕ ϕ , ϕ

∈

∈

∈

−

∈

+

=

∇

∇ +

∇

∇

⁽²²⁾

kde I je množina všech index odpovídající uzl m element diskretizované oblasti

Ω

. IN je množina index na hranici oblasti, kde je zadána Neumannova okrajová

podmínka a I_D je množina index na hranici oblasti, kde je zadána Dirichletova okrajová podmínka.

Dále lze leny

( ϕ

i

, ϕ

j

)

nahradit maticí

Φ

, koeficienty

Q

ⁱ vektorem

Q

_,

koeficienty

G

ⁱ vektorem

G

. Získáme tvar:

( ^{∇ ∇} ) ⁺ ( ^{∇ ∇} ) ^{= Φ + Φ}

∈

−

∈

T T

j i I

i i j

i I

I i

i

K P K Q G

P

D D

ϕ ϕ ϕ

ϕ , ,

⁽²³⁾

A sou et len na pravé stran nahradíme jedním vektorem

b

.

b G

Q

^T

Φ +

^T

Φ =

⁽²⁴⁾

Dále se budeme zabývat levou stranou rovnice (22). V ní vystupují dva leny, len

(

i j

)

I I i

i

K

P

D

ϕ ϕ ∇

∇

−

∈

,

obsahuje neznámé hodnoty

P

ⁱ a len

(

i j

)

I i

i

K

P

D

ϕ ϕ ∇

∇

∈

,

obsahuje známé hodnoty tlakových výšek, zadané v Dirichletov okrajové podmínce. Prozatím len se zadanými okrajovými podmínkami vynecháme, len rovnice (22)

( K ∇ , ϕ

i

∇ ϕ

j

)

nahradíme maticí

A ~

, známé i neznámé koeficienty

P

ⁱ

lze vyjád it jako vektor

P

. Získáme tak systém stavových rovnic:

b P

A ~ ⋅ = ~

⁽²⁵⁾

Do systému (11) ješt chybí dodat Dirichletovy okrajové podmínky. Pro p ípad zadané tlakové výšky v uzlu je pot eba vzít všechny nenulové prvky matice

A ~

daného uzlu (každý uzel má jeden ádek a jeden sloupec v matici) a jednotliv je vynásobit se zadanou hodnotou tlakové výšky a ode íst od p íslušného ádku vektoru pravé strany.

Následn je nutné vynulovat všechny prvky v ádku a ve sloupci náležící danému uzlu a pouze na diagonálu vložit hodnotu 1. Na ádek vektoru pravé strany, který pat í

danému uzlu, se vloží zadaná hodnota tlakové výšky. Po provedení tohoto algoritmu získáme kone ný tvar lineárního systému rovnic:

b P

A ⋅ =

⁽²⁶⁾

Vektor neznámých hodnot tlakových výšek

P

vypo ítán eši em lineárních systém , v našem p ípad je využit Matlab.

Matici

A

nazýváme stavovou nebo také globální maticí, je složena z tzv.

lokálních matic, což jsou jednotlivé p ísp vky k výrazu

( K ∇ , ϕ

i

∇ ϕ

j

)

po ítané pro jednotlivé elementy celé diskretizované oblasti (sít ).

Bázové funkce používáme lineární ve tvaru závislém na dimenzi elementu. Pro úse ku:

( ) ^x _{i i} ^x

i α β

ϕ = +

⁽²⁷⁾

pro trojúhelník:

( ) ^{x y} _{i i} ^x _i ^y

i α β γ

ϕ , = + +

⁽²⁸⁾

a pro ty st n:

( ^{x y z} ) _{i i} ^x _i ^y _i ^z

i α β γ δ

ϕ , , = + + +

⁽²⁹⁾

Její koeficienty

α

,

β

,

γ

a p ípadn

δ

dopo ítáváme ze známých hodnot, které musí bázová funkce spl ovat nad danými uzly jednotlivých element . Funkce musí spl ovat tuto podmínku:

j ij

i x δ

ϕ =

⁽³⁰⁾

kde

ϕ i

je bázová funkce i-tého uzlu,

x

j je j-tý uzel,

δ ij

je Kroneckerovo delta ( = 1 pro i = j; = 0 pro i j).

Detailn jší rozepsání tvar bázových funkcí, výpo t jejich koeficient i struktury lokálních matic jednotlivých element je možné najít v [6].

Všeobecn lze uvést, že výhodou primární formulace je její relativn snadná realizace – jednoduché odvození formulace úlohy a nep íliš složitý algoritmus sestavení globální matice. Diskretizací rovnice ešíme soustavu pro primární veli inu p s menším

po tem prom nných, než tomu je ve smíšené nebo smíšené-hybridní formulaci MKP a navíc je stavová matice primární formulace symetricky pozitivn definitní a tím urychluje možnost výpo tu na rozdíl od MH formulace, jejíž stavová matice je indefinitní.

Nevýhodou však je skute nost, že vypo teme primární veli inu p a filtra ní rychlost u musíme následn vypo ítat zp tným dosazením do Darcyho zákona. Tím je snížena hladkost, filtra ní rychlost se stává po elementech konstantní. Další nevýhodou je, že bilan ní rovnice je v této formulaci spln na v obecném p ípad na celé oblasti a v každé oblasti, která obklopuje uzel, ale není spln na na jednotlivých elementech v síti, a proto není možné pomocí této formulace po ítat transport látky, který tento p edpoklad vyžaduje.

P esto poskytuje dostate né výsledky k posouzení hydrogeologické situace.

1.3 Po íta ová implementace modelu

Po seznámení s matematicko fyzikálním modelem uvedeném v p edchozí kapitole si p edstavíme program Flow123D, který je po íta ovou implementací tohoto modelu. Je psán v programovacím jazyku C a v sou asné dob umož uje krom výpo tu ustáleného proud ní také simulaci neustáleného proud ní nebo transportu látek v podzemí. V naší práci vyvíjíme tu ást programu, která se zabývá výpo tem ustáleného proud ní.

Pro náš model je navíc tento program výhodný z toho d vodu, že je v n m již implementován algoritmus pro vytvo ení pot ebného *.ngh souboru, seznamu všech sousedností element z celé sít . Postup spušt ní tohoto algoritmu, stejn tak postup spušt ní výpo tu ustáleného proud ní je uveden v p íloze této práce.

Charakteristickou vlastností programu Flow123D je použití kombinovaného (multidimenzionálního) p ístupu. Používá prvky t í dimenzí, kde 3D odpovídá poréznímu médiu, 2D puklinám a 1D kanál m. Program umož uje dva zp soby propojení a to bu kompatibilní, nebo nekompatibilní. Prvky r zných dimenzí jsou propojeny kompatibiln tehdy, má-li prvek s vyšší dimenzí celou jednu st nu spole nou práv s celým jedním prvkem nižší dimenze. Pro ilustraci jsou uvedeny p íklady spojení. V rozsáhlých úlohách je mnohdy nevyhnutelné prvky spojovat nekompatibiln a to z d vodu vysoké algoritmické náro nosti kompatibilního spojení prvk . V této práci je uvažováno pouze kompatibilní spojování.

Obr. 1-3 - kompatibilní spojení 2D a 1D Obr. 1-4 - nekompatibilní spojení 2D a 1D

Obr. 1-5 - kompatibilní spojení 3D a 2D Obr. 1-6 - nekompatibilní spojení 3D a 2D

2 NÁVRH A IMPLEMENTACE ALGORITMU

V této ásti práce si rozd líme problém na dv ásti:

o Návrh algoritmu pro kombinaci 2D a 1D element o Návrh algoritmu pro kombinaci 3D a 2D element

P ed popisem zmín ných ástí této práce je však t eba vysv tlit zp sob, jakým byl realizován kombinovaný model doposavad.

P vodní implementace primární formulace má komunikaci mezi elementy r zných dimenzí uskute n nou díky spole ným uzl m sousedících element . Do globální matice se zapisují prvky lokálních matic všech element a to i r zných dimenzí, jak je nazna eno na obr. 2-1. Prvky lokální matice 1D elementu (zelen ) jsou v globální matici se teny s p íslušnými prvky lokální matice 2D elementu ( erven ).

Obr. 2-1, sestavení globální matice

Zp sob, jakým je dosaženo komunikace mezi 2D a 3D elementem, je analogický k výše uvedenému, rozdíl je pouze v rozm ru lokálních matic element , pro ty st n je to matice tvrtého ádu a pro trojúhelník z stává matice t etího ádu, jako tomu je na obr. 2-1. Tím v globální matici s ítáme 9 pozic (p ísp vek 2D elementu) namísto výše uvedených 4 pozic (p ísp vek 1D elementu).

Podle výsledk testu pro p í ný p etok p es element nižší dimenze (puklinu) z úvodu této práce je tento zp sob komunikace nevhodný. P í ina spo ívá v tomto: na obr. 2-2 jsou pro elementy 0, 1 a 2 spole né uzly 0 a 3 a dochází v toku sítí nap í

liniovým elementem (nap . od uzlu 2 k uzlu 1) a jelikož je primární formulace založena na výpo tu tlakových výšek na jednotlivých uzlech v síti, tak vypo tená tlaková výška na uzlu 0 a 3 je stejná jak pro liniový element, tak i pro oba trojúhelníkové elementy, nezávisle na propustnosti linie. Jednoduše vysv tleno, tok z elementu 1 „p esko í“

rovnou na element 0 a nezohlední se tím již zmín ná propustnost elementu 2.

Obr. 2-2, spole né uzly pro 1D a 2D elementy

ešením tohoto problému je následující postup. Na obr. 2-3 je v levé ásti nazna eno p vodní spojení 2D a 1D element . Na pravé stran nazna ujeme, že 1D element od 2D elementu fyzicky odd líme. Tím je pro ádný výpo et t eba dodefinovat velikost p etoku mezi t mito elementy, který je úm rný rozdílu tlakových výšek.

Použijeme tento vztah:

− + +

= +

2 3

2 1 3 2

1II II II I I

c

p p p p σ p

u ⁽³¹⁾

Kde

u [ ]

^ms−1 je p etok mezi 2D a 1D elementem, σ_c

[ ]

^s−1 je p estupový koeficient mezi elementy, p₁^II

[ ]

^m je tlaková výška v uzlu íslo 1 ve 2D elementu a p₁I

[ ]

^m je tlaková výška v uzlu 1 v 1D elementu. První len v závorce ve vztahu (31) je celková tlaková výška 2D elementu a druhý len je celková tlaková výška 1D elementu.

V podstat se jedná o analogii Newtonovy okrajové podmínky uvedené ve vztahu (5) první kapitoly s tím rozdílem, že ji používáme pro p etok mezi elementy r zných dimenzí uvnit sít .

Obr. 2-3, princip odd lení element 1D a 2D

ím menší hodnota p estupového koeficientu σ_c, tím dochází u dané sousednosti element r zných dimenzí k menšímu p etoku.

Tímto jsme uvedli zp sob zpracování nové komunikace mezi elementy r zných dimenzí obecn a v následujícím textu toto rozvedeme pro kombinaci 1D s 2D elementy a 2D s 3D elementy.

2.1 Návrh algoritmu pro kombinaci 2D a 1D element

ešení spo ívá ve vytvo ení nových uzl v místech napojení 2D element na 1D elementy, tím se docílí jednozna ného odd lení 1D a 2D element . Elementy spolu v modelu nebudou již propojeny p es spole né uzly, ale p es p estupový koeficient σ_c charakterizující množství vody p eteklé mezi elementy (uzly).

P ipome me, že primární model MKP vede na soustavu n rovnic o n neznámých v maticovém zápisu:

b P A ⋅ =

Kde A je globální matice, P jsou hledané tlakové výšky v uzlech sít a b je vektor pravé strany. Takovéto propojení element funguje jako zdroj

q

, který je na pravé stran rovnice.

Ze vztahu (31) podle obr.2-3 získáme na pravé stran rovnice (26):

− + +

− +

=

=

= 3 3 2

1 ₁^II ₂^II ₃^II ₁^I ₂^I

c

p p p p

σ

p

II3 II2

II1 u u

u ⁽³²⁾

− + + + +

=

= 2 3 2

1 ₁^II ₂^II ₃^II ₁^I ₂^I

c

p p p p

σ

p

2I I1 u

u ⁽³³⁾

po rozepsání:

6 6

9 9

9

2 1

3 2

1 I

c I c II c II c II

cp

σ

p

σ

p

σ

p

σ

p

σ

_{− − + +}

−

=

II1

u ⁽³⁴⁾

6 6

9 9

9

2 1

3 2

1 I

I c II c

II c II c

cp

σ

p

σ

p

σ

p

σ

p

σ

_{− − + +}

−

=

II2

u

6 6

9 9

9

2 1

3 2

1 I

I c II c

II c II c

cp

σ

p

σ

p

σ

p

σ

p

σ

_{− − + +}

−

=

II3

u

4 4

6 6

6

2 1

3 2

1 I

c I c II c II c II

cp

σ

p

σ

p

σ

p

σ

p

σ

_{+ + − −}

+

=

I1

u

4 4

6 6

6

2 1

3 2

1

I c I c II c II c II

cp

σ

p

σ

p

σ

p

σ

p

σ

_{+ + − −}

+

=

2I

u

Výše uvedené vztahy p evedeme na levou stranu rovnice (26) a v globální matici se projeví takto:

=

+ +

−

−

−

+ +

−

−

−

−

− +

+ +

−

− +

+ +

−

− +

+ +

I I II II II

I I II II II

c c

c c

c

c c

c c

c

c c

c c

c

c c

c c

c

c c

c c

c

r r r r r

p p p p p

2 1 3 2 1

2 1 3 2 1

4 4

6 6

6

4 4

6 6

6

6 6

9 9

9

6 6

9 9

9

6 6

9 9

9

σ σ

σ σ

σ

σ σ

σ σ

σ

σ σ

σ σ

σ

σ σ

σ σ

σ

σ σ

σ σ

σ

(35)

Všechny hodnoty v globální matici budou tedy p i teny (ode teny) na dané pozice ke stávajícím hodnotám získaným z lokálních matic element . Hodnoty

r

jsou

práv p íslušným p etokem

u

. Po et uzl , které je t eba kolem originálního uzlu vytvo it pro jednozna né odd lení element r zných dimenzí, je roven po tu liniových element p ipojených k danému uzlu. Výjimkou jsou krajní uzly 1D element , kde místo jednoho nového uzlu p idáváme uzly dva obr. 2-4. Nový uzel znamená uzel s identicky stejnými geometrickými vlastnostmi jako má originální uzel, tj. na obrázcích teoreticky nejsou nové uzly vid t, proto je pro p ehlednost zobrazujeme s malým posunutím sm rem k element m, se kterými mají být spojeny.

Obr. 2-4, znázorn ní p idávání krajních uzl 1D elementu

Pro názornost použijme originální sí z typového p íkladu v úvodu. Vidíme, že uzly liniových element jsou zárove uzly trojúhelníkových element . Cílem nového algoritmu je zobrazení na obr.2-5. Jak je z obrázku patrné, originální uzly ponecháváme liniovým element m a nové uzly p i azujeme trojúhelníkovým element m. Každý uzel má své identifika ní íslo, o íslování nových uzl je provedeno tak, že se p e ísluje celá sí od za átku i s násobnými uzly jako na obr. 2-5.

Obr. 2-5, požadovaná sí

Více je celý postup b hu algoritmu pospán níže, v kapitole 2.3.

2.2 Návrh algoritmu pro kombinaci 3D a 2D element

U kombinace 2D a 3D element je situace analogická výše popsané kombinaci 1D a 2D element . Op t dostane každý element v místech styku s elementem jiné dimenze nový uzel. Rovnice (31) se p em ní na:

+

− + + +

= +

3 4

3 2 1 4 3 2

1III III III III II II II

c

p p p p p p σ p

u ⁽³⁶⁾

Kde

u [ ]

^ms−1 je p etok mezi 3D a 2D elementem, σ_c

[ ]

^s−1 je p estupový koeficient mezi elementy, p₁^III

[ ]

m je tlaková výška v uzlu íslo 1 ve 3D elementu a p₁II

[ ]

^m je tlaková výška v uzlu 1 v 2D elementu. První len v závorce ve vztahu (36) je celková tlaková výška 3D elementu a druhý len je celková tlaková výška 2D elementu.

Kombinace t chto dvou dimenzí je nazna ena na obr.2-5.

Obr. 2-6, princip odd lení element 2D a 3D

Stejn jako u p edešlé kombinace si odvodíme koeficienty, které se budou p i ítat (ode ítat) od stávajících hodnot na daných pozicích v globální matici.

Vztah (36) si op t p epíšeme do tvaru:

+

− + + +

− +

=

=

=

= 4 4 3

1 ₁^III ₂^III ₃^III ₄^III ₁^II ₂^II ₃^II

c

p p p p

p p

σ

p

III4 III3 III2

III1 u u u

u ⁽³⁷⁾

+ +

− + + + +

=

=

= 3 3 4

1 ₁^II ₂^II ₃^II ₁^III ₂^III ₃^III ₄^III

c

p p p p p p

σ

p

II3 II2

II1 u u

u ⁽³⁸⁾

Po rozepsání:

12 12

12 16

16 16

16

3 2

1 4

3 2

1 II

II c II c

III c III c

III c III c

cp

σ

p

σ

p

σ

p

σ

p

σ

p

σ

p

σ

_{− − − + + +}

−

=

III1

u ⁽³⁹⁾

12 12

12 16

16 16

16

3 2

1 4

3 2

1 II

II c II c

III c III c

III c III c

cp

σ

p

σ

p

σ

p

σ

p

σ

p

σ

p

σ

_{− − − + + +}

−

=

III2

u

12 12

12 16

16 16

16

3 2

1 4

3 2

1 II

c II c II c III c III c III c III

cp

σ

p

σ

p

σ

p

σ

p

σ

p

σ

p

σ

_{− − − + + +}

−

=

III3

u

12 12

12 16

16 16

16

3 2

1 4

3 2

1 II

c II c II c III c III c III c III

cp

σ

p

σ

p

σ

p

σ

p

σ

p

σ

p

σ

_{− − − + + +}

−

=

III4

u

9 9

9 12

12 12

12

3 2

1 4

3 2

1 II

c II c II c III c III c III c III

cp

σ

p

σ

p

σ

p

σ

p

σ

p

σ

p

σ

_{+ + + − − −}

+

=

II1

u ⁽⁴⁰⁾

9 9

9 12

12 12

12

3 2

1 4

3 2

1 II

c II c II c III c III c III c III

cp

σ

p

σ

p

σ

p

σ

p

σ

p

σ

p

σ

_{+ + + − − −}

+

=

II2

u

9 9

9 12

12 12

12

3 2

1 4

3 2

1 II

c II c II c III c III c III c III

cp

σ

p

σ

p

σ

p

σ

p

σ

p

σ

p

σ

_{+ + + − − −}

+

=

II3

u

A protože tlakové výšky

p

jsou neznámé, výše uvedené vztahy p evedeme na levou stranu rovnice (26) a v globální matici se projeví takto:

=

+ +

+

−

−

−

−

+ +

+

−

−

−

−

+ +

+

−

−

−

−

−

−

− +

+ +

+

−

−

− +

+ +

+

−

−

− +

+ +

+

−

−

− +

+ +

+

II II II III III III III

II II II III III III III

c c

c c

c c

c

c c

c c

c c

c

c c

c c

c c

c

c c

c c

c c

c

c c

c c

c c

c

c c

c c

c c

c

c c

c c

c c

c

r r r r r r r

p p p p p p p

3 2 1 4 3 2 1

3 2 1 4 3 2 1

9 9

9 12

12 12

12

9 9

9 12

12 12

12

9 9

9 12

12 12

12

12 12

12 16

16 16

16

12 12

12 16

16 16

16

12 12

12 16

16 16

16

12 12

12 16

16 16

16

σ σ

σ σ

σ σ

σ

σ σ

σ σ

σ σ

σ

σ σ

σ σ

σ σ

σ

σ σ

σ σ

σ σ

σ

σ σ

σ σ

σ σ

σ

σ σ

σ σ

σ σ

σ

σ σ

σ σ

σ σ

σ ⁽⁴¹⁾

Hodnoty

r

jsou hodnoty pravé strany p íslušného uzlu, které jsou sou tem p vodní pravé strany a práv p etokem

u

. Tím je hotov návrh komunikace 1D s 2D elementy stejn tak 2D s 3D elementy a m žeme p istoupit k popisu implementace t chto myšlenek do struktur programu Flow123D.

3 IMPLEMENTACE ALGORITMU 3.1 Použití primárního modelu

Podívejme se, jakým zp sobem je náš algoritmus patrný z hlediska uživatele programu Flow123D.

Podle obr. 3-1 je p ed samotným spušt ním programu pot eba zadat do INI souboru názvy soubor sít (*.msh), okrajových podmínek (*.bcd), sousedností (*.ngh) a materiál (*.mtr). Poznamenejme ješt , že soubor sousedností dokáže Flow123D vytvo it také, ovšem za p edpokladu, že zbývající t i soubory již máme vytvo ené.

Obr. 3-1, vývojový diagram

Pokud máme tyto soubory zadány do *.ini souboru, je možné spustit program.

Podrobný popis spušt ní programu stejn jako popis všech pot ebných soubor lze nalézt v p íloze.

Pokud jsou v síti, na které se po ítá, p ítomny elementy r zných dimenzí spojené spole nými uzly, dojde ke spušt ní algoritmu p idávání uzl . Jeho struktura je vysv tlena níže. Výsledkem našeho kódu je vytvo ení nové sít s p idanými uzly uložené do nového souboru sít „m.msh“. Zárove je vytvo en i nový soubor okrajových podmínek „m.bcd“. Následn program automaticky p epíše v *.ini soubor na „m.ini“ (z d vodu zálohy p vodního *.ini souboru) a v n m název p vodního souboru sít na „m.msh“ a stejn tak p vodní soubor okrajových podmínek na vytvo ený „m.bcd“ a spustí program Flow123D znovu. Tentokrát dojde v výpo tu úlohy proud ní, jejíž výsledek zapíše do souboru „m.pos“.

3.2 Algoritmus z hlediska programátora

Na obr. 3-2 je nazna en chod programu tak, jak byl do programu Flow123D implementován. Následn si popíšeme každý krok tohoto algoritmu. Je koncipován do dvou fází. V první fázi je nejprve otestována sí , zda je nebo není t eba odd lovat elementy r zných dimenzí, tuto fázi zobrazuje obr. 3-2. Pro p ípad, že je t eba odd lovat, vytvo í již zmín né dva soubory „m.msh“ a „m.bcd“. Struktura obou soubor je uvedena v p íloze této práce.

Ve druhé fázi na obr. 3-5 je již sí vyhodnocená jako zp sobilá k výpo tu, ale je ješt pot eba zapsat správné hodnoty p estupových koeficient na správné pozice v globální matici, jak je uvedeno v (35) resp. (41). Nejprve za n me s popisem první fáze.

Obr. 3-2, 1. fáze, vývojový diagram

3.2.1 Multinode_test

Pro v tšinu celého algoritmu byl pro p ehlednost vytvo en nový *.c soubor nazvaný multinode.c. Obsah tohoto souboru nyní popíšeme. Po na tení vstupních soubor úlohy a provedení všech výpo t týkajících se topologie sít a element dochází v procedu e multinode_test k projití všech uzl sít a pokud je minimáln na jednom z nich p ipojen zárove 2D a 1D nebo 2D a 3D element, pak je pot eba provést náš algoritmus.

V tomto p ípad je stále v této procedu e zjišt n po et uzl , který je t eba p idat ke každému originálnímu uzlu pomocí procedury Add_nodes_1D_2D resp.

Add_nodes_2D_3D. Tyto procedury využívají výstupu funkce How_many_planes. Ta vrací po et rovin napojených na daný uzel, po et element p ipojených k danému uzlu a vý et identifika ních ísel t chto element v každé rovin . Rovina je tvo ena vždy 2D elementy. V pojetí 3D prost edí tedy rovina vyjad uje puklinu a v pojetí 2D prost edí vyjad uje libovoln nato enou rovinu porézního media.

3.2.2 Add_nodes

Tyto procedury p idávají nové uzly k originálním uzl m tím zp sobem, že pouze zvyšují hodnotu p íznaku added_nodes ve struktu e každého uzlu, tím si každý uzel s sebou nese informaci o po tu pot ebných uzl . Po et nových pot ebných uzl z hlediska struktury sít je pro kombinaci 1D a 2D element roven po tu liniových element ležících v jedné rovin a p ipojených k danému uzlu, výjimku tvo í pouze uzly, na které je p ipojena pouze jedna linie jako na obr. 2.4, pak je z ejmé, že se jedná o krajní linii (kraj pukliny) v oblasti a k tomuto uzlu p idáváme dva nové uzly z toho d vodu, aby puklina ani na jejím okraji nemohla být obtékána.

U kombinace 2D a 3D element je po et nových pot ebných uzl roven pro p ípad jedné roviny procházející daným uzlem dv ma p idaným uzl m, pro dv roviny procházejícím tímto uzlem je roven ty em novým uzl m a pro t i roviny procházející tímto uzlem je roven osmi novým p idaným uzl m, více rovinám procházející jedním uzlem se tato práce nezabývá.

Dále jsou v rámci multinode_test volány i procedury Create_elements, Create_msh_file a Create_bcd_file. V p ípad , že není t eba vytvá et nové uzly (na žádném z uzl sít není ani jeden uzel spole ný pro 1D a 2D elementy resp. 2D

3.2.3 Create_elements

Úkolem této procedury je zavést nové íslování uzl assign_multinodes, alokovat v pam ti místo pro nové uzly element alocate_new_nodes a p i adit tyto nové uzly element m fill_new_nodes.

Po p idání nových uzl každému originálnímu uzlu vyvstala otázka, jak p istoupit k o íslování nových uzl v novém m.msh souboru. Nabízely se dv varianty.

První byla p vodní o íslované uzly ponechat a o íslovat pouze nové uzly a druhá byla íslovat všechny postupn od nuly, což by vedlo k tomu, že by uzly se stejnou polohou byly hodnotami identifika ních ísel u sebe. Zvolili jsme z d vodu p ehlednosti a snazšího dohledávání v rozsáhlejších sítích druhou variantu, p e íslovat vše od nuly.

3.2.4 Assign_multinodes

V této procedu e se alokují místa v pam ti pro novou celo íselnou prom nnou node_number ve struktu e uzlu, do které se zapíše jeho nové identifika ní íslo. Po et alokovaných míst je roven po tu pot ebných uzl k p idání (added_nodes) + originální uzel, který je znásobován. Takže uzel, ke kterému není t eba p idávat žádný násobný uzel, bude mít pouze jednu prom nnou node_number, ve které bude uloženo jeho nové id. íslo.

3.2.5 Alocate_new_nodes

V této procedu e se již zabýváme elementy, alokujeme zde u každého elementu sít celo íselné prom nné element_node. Jejich po et je roven po tu uzl daného elementu (linie dva, trojúhelník t i, ty st n ty i). Do této sady prom nných je t eba zapsat identifika ní ísla nových uzl , které tvo í daný element. Tím lze pak jednoduše vytvo it pot ebný m.msh soubor.

3.2.6 Fill_new_nodes

Poslední procedura je nejrozsáhlejší a nejsložit jší z celého algoritmu, protože je zde nutné p i adit ty správné uzly k jednotlivým element m (neboli naplnit v p edešlé procedu e vytvo ené sady prom nných element_node jednotlivých element ). Prochází se uzel po uzlu v celé síti a testuje se, zda je daný uzel znásobený. Nastávají dv možnosti výsledku testu. Bu je daný uzel samotný, pak je rovnou p i azen element m,