KURS IZRADE MODELA LIZZA – WODA
NA PRIMERU TESTA BUNARA PORED REKE
Cilj: izrada HD modela podzemnih voda, na kome se simulira testiranje bunara
Pretpostavka: Slušaoci su završili izradu konceptualnog – hidrogeološkog modela.
Pre kursa svaki polaznik mora da ima instaliranu Lizu, MingW, Mintty, Wodu i Paraview.

1. OSNOVNI POJMOVI O SOFTVERU

1.1. Procesi i osnovne jednačine

1.2. Generalizovani problemi koje je moguće rešavati

2. KONTURE MODELA

2.1 Karta terena
Izbor (satelitski snimak, skenirana karta: topografska, geološka, hidrogeološka, geomorfološka, ...), veličina, razmera (detaljnost).
*Georeferenciranje
*Postavljanje

2.2 Postavljanje kontura modela
Pojam kontura modela, kao graničnog uslova.
Veličina terena obuhvaćenog modelom. Posledice.
Napomena: na izrađenom modelom uporediti rezultate proračuna sa različitim konturama modela.
CSV format.
*Postavljanje više varijanti kontura modela, sa definisanom ulogom delova kontura.

2.3 Planare (granične površi) slojeva
Šta su granične površi slojeva? Šta ih čini takvim (hidrogeološki, hidrodinamički)?
Kako unosimo granične površi? -> Bušotine (Gaus-Krigerove koordinate, nadmorska visina, ili drugi pravougli koordinatni sistem), količina, distribucija u planu, dubina (obuhvaćeni slojevi).
Lokalno izostajanje sloja, sočivo (međusloj). Interpolacija i konflikt površi.
Kako se formira mreža modela u planu i zašto je to važno. OVO BIH STAVIO KASNIJE
* Unošenje podataka za granične površi slojeva

2.4 Slojevi modela
Pojam sloja u modelu i veza sa prirodnim slojem. Parametri modela: Kx, Ky, Kz, S, Syld.
Vertikalna diskretizacija. OVO BIH OSTAVIO ZA PRICU O NUMERICI
Početni raspored parametara, interpolacija.
Promena parametara u prostoru (vremenski nisu promenljivi): u planu, po profilu (po visini linearno).
*Unošenje podataka za slojeve.

2.5 Jednačina strujanja u zasićenoj sredini
Darsijev zakon originalno je formulisan za jednodimenzionalno strujanje:
u = − K h izlaz − h ulaz d ,	(1)
gde je d pređeno rastojanje, u je brzina strujanja, h = z + p − p a ρ g je hidraulički potencijal, z je visina merena u odnosu na proizvoljnu ali jedinstvenu kotu, p je pritisak u vodi, p a je atmosferski pritisak, ρ je gustina vode, g je intenzitet gravitacionog ubrzanja, K je koeficijent filtracije
K = k ρ g μ ,	(2)
k je koeficijent propusnosti (karakteristika sredine), a μ je koeficijent viskoznosti. Prelaskom na graničnu vrednost kada d → 0 dobija se Darsijev zakon u diferencijalnom obliku:
u = − K dh dx .	(3)
Darsijev zakon je primenljiv i u višedimenzionalnom slučaju i tada glasi
u = − K ∇ h.	(4)
Ovde je u vektor brzine, a ∇ h je gradijent hidrauličkog potencijala.
U izotropnoj sredini koeficijent filtracije je skalar, međutim ako je sredina neizotropna onda su K i k simetrični pozitivno-definitni tenzori drugog reda tj. simetrične pozitivno-definitne matrice
K = K xx K xy K xz K xy K yy K yz K xz K yz K zz .	(5)
Pozitivna definitnost je osobina tenzora K koja garantuje da će vektor K ∇ h imati sve pozitivne koordinate ako je takav i vektor ∇ h . Ako ovo ne bi važilo voda bi mogla da teče u smeru u kome raste pritisak, što nije fizičko ponašanje.
Pozitivna definitnost i simetrija matrice K proističu iz njenog prikaza kao proizvoda matrica rotacije i dijagonalne matrice koeficijenta filtracije:
K = R − 1 K xx ¯ 0 0 0 K yy ¯ 0 0 0 K zz ¯ R , K xx ¯ > 0, K yy ¯ > 0, K zz ¯ > 0,	(6)
gde je R matrica rotacije
R = R xy R xz R yz = cos α − sin α 0 sin α cos α 0 0 0 1 cos β 0 − sin β 0 1 0 sin β 0 cos β 1 0 0 0 cos γ − sin γ 0 sin γ cos γ ,	(7)
a α, β i γ su usmereni uglovi za koje je potrebno rotirati koordinatne ose u xy, xz i yz ravni respektivno da bi se one dovele do poklapanja sa osama neizotropije.
Ukoliko se ose neizotropije poklapaju sa koordinatnim osama, rotacija nije potrebna i matrica K je dijagonalna. Jedino ovaj slučaj je implementiran u Lizzi.
Koeficijent filtracije se zadaje nezavisno u svakom sloju u ključnim tačkama između kojih se vrši linearna interpolacija. Moguće je zadati i linearnu varijaciju po visini u okviru sloja. Povrh toga, mogu se zadati zone u okviru sloja koje imaju drugi koeficijent filtracije. U ovom slučaju keficijent filtracije ima prekid u okviru sloja.
Jednačina održanja mase u poroznoj sredini glasi
∂ θ ρ ∂ t + ∇ ⋅ { ρ u } = ρ q ,	(8)
gde q predstavlja izvore i ponore, t je vreme, a θ je sadržaj vode tj. deo jedinične zapremine koji sadrži vodu. U zasićenom slučaju ovo je isto što i poroznost ε .
Ako pretpostavimo da je gustina konstantna i da nema izvora i ponora, jednačina održanja mase postaje
∂ θ ∂ t + ∇ ⋅ u = 0.	(9)
Da bi se razumeo matematički zapis zakona održanja neophodno je shvatiti suštinu pojma divergencije, koja se sastoji u teoremi o divergenciji:
∫ Ω ∇ ⋅ u d Ω = ∮ ∂ Ω u ⋅ n d Γ ,	(10)
gde je n jedinični vektor normalan na granicu ∂ Ω usmeren ka spolja (tzv. spoljašnja normala), a u ⋅ n predstavlja brzinu isticaja iz Ω . Ova jednakost govori da je zapreminski integral divergencije brzine u jednak zbiru proticaja kroz granice te zapremine. Ako se jednakost podeli sa zapreminom Ω i ako se pusti da ova zapremina teži nuli, dobija se sa leve strane ∇ ⋅ u . Dakle
∇ ⋅ u = lim ∣ Ω ∣ → 0 1 ∣ Ω ∣ ∮ ∂ Ω u ⋅ n d Γ .	(11)
Stoga o divergenciji treba razmišljati kao o sumi proticaja u beskonačno malu zapreminu i iz nje, a ne kao o sumi parcijalnih izvoda
∇ ⋅ u = ∂ u x ∂ x + ∂ u y ∂ y + ∂ u z ∂ z .	(12)
Zamenom Darsijevog zakona u jednačinu (9) dobija se jednačina podzemnog strujanja
S ∂ h ∂ t = ∇ ⋅ ( K ∇ h ) ,	(13)
gde je S koeficijent specifične izdašnosti izdani pod pritiskom, tj. S = ∂ ε ∂ h .
Ako se strujna slika ne menja u vremenu, radi se o stacionarnom rešenju koje odgovara slučaju kada je vremenski izvod jednak nuli. Nasuprot tome, rešenje koje se menja u vremenu se naziva tranzijentno rešenje.

2.6 Nezasićena sredina
Zasićenost vodom s se definiše kao količnik zapremine vode i zapremine koju zauzimaju svi fluidi zajedno (tj. vazduh i voda). Prema tome
θ = s ε .	(14)
U zasićenoj sredini s = 1 , a u nezasićenoj sredini 0 < s < 1 .
Koeficijent filtracije zavisi od zasićenosti i manji je u nezasićenoj sredini. Ukoliko je zasićenost veoma mala, vode nema dovoljno da bi uopšte tekla. Darsijev zakon za nezasićenu sredinu tako postaje
u = − k r ( s ) K ∇ h ,	(15)
gde je K je koeficijent filtracije u zasićenom slučaju, k r ( s ) je relativna propustljivost, ψ = h − z je dubina, a z je visina u odnosu na proizvoljno izabranu kotu, jedinstvenu za ceo model. Predstavljanje tenzora koeficijenta filtracije kao proizvoda skalarne funkcije k r ( s ) i konstantnog tenzora K je opravdano u izotropnom slučaju, ali postoje dokazi da u neizotropnom slučaju komponente tenzora različito menjaju u zavisnosti od zasićenosti, i dakle ovakav prikaz nije adekvatan. Međutim, ovakav pristup je uobičajen, alternativni pristupi su u povoju i samo ovo je implementirano.
Privlačne sile koje deluju između molekula vode su u ravnoteži pod vodom, gde je svaki molekul vode okružen istim takvim molekulima vode. Međutim, budući da su privlačne sile između molekula vazduha i molekula vode mnogo slabije nego između dva molekula vode, rezultanta privlačnih sila na površini vode je različita od nule i orijentisana suprotno od vazduha. Pritisak koji nastaje delovanjem ove sile na površinu vode naziva se porni pritisak. Budući da deluju u istom smeru na površinu vode, porni pritisak p c i pritisak vazduha p a su u ravnoteži sa pritiskom u vodi p w :
p c = p a − p w .	(16)
Odavde se vidi da je pritisak u vodi tik do površine manji od pritiska vazduha za porni pritisak.
Porni pritisak se može izraziti kao
p c = 2 σ cos θ r ,	(17)
gde je r poluprečnik pore, θ je specifični kontaktni ugao koji slobodna površina zaklapa sa čvrstim telom, a σ je površinski napon koji ne zavisi od veličine pore. Iz ove formule se vidi da površinski napon neograničeno raste kako se smanjuje poluprečnik pore. Kako se tlo isušuje, voda se povlači u sve manje pore i porni pritisak raste. U prirodi bi voda isparila kada pritisak u njoj postane dovoljno mali, ali kako isparavanje nije deo ovog fizičkog modela, porni pritisak može da raste neograničeno a pritisak u ostacima vode može da postane proizvoljno mali, čak i negativan. U praksi ovo ne predstavlja problem zato što u toj situaciji relativna propustljivost teži nuli što zaustavlja tok, sprečavajući dalje isušivanje.
Tok kroz nezasićenu sredinu u suštini je tok dve faze – vode i vazduha. Za modeliranje dvofaznog toka potrebne su dve jednačine podzemnog strujanja – jedna za vazduh i jedna za vodu – međusobno povezane jednačinom (16). Međutim, uobičajeni pristup proračunu strujanja podzemne vode podrazumeva da se pretpostavi da je pritisak vazduha konstantan i jednak atmosferskom pritisku. Ovo eliminiše jednačinu strujanja vazduha.
Da bi se zatvorio sistem, ostaje da se odredi zavisnost zasićenosti od dubine i funkcija k r ( s ) . Postoji više modela kojima se izražavaju ove zavisnosti. Jedan od njih je Van Genuchten – Mualem model:
s e = { 1 [ 1 + ∣ α ψ ∣ n ] m ako je ψ < 0 1 ako je ψ > 0 ,	(18)
k r = s e 1 2 { 1 − [ 1 − s e 1 m ] m } 2 .	(19)
Ovde je s e = s − s r s s − s r efektivna zasićenost, s r je rezidualna zasićenost koja odgovara beskonačnom pornom pritisku, s s je maksimalna moguća zasićenost (praktično uvek 1), n i α su parametri modela, a m = 1 − 1 n . Iako su parametri n i α zapravo karakteristike sredine i trebalo bi da variraju od sloja do sloja, trenutno korisniku nije omogućeno da izabere ove vrednosti, već su one fiksirane u programu.
Konačno, jednačina podzemnog strujanja u nezasićenoj sredini glasi
∂ θ ∂ t = ∇ ⋅ ( k r ( ψ ) K ∇ h ) .	(20)
Ovo je nelinearna jednačina i predstavlja veliki izazov za numeričko rešavanje, najviše zbog toga što u slabo zasićenoj sredini pritisak može neograničeno da opada a da to ne dovede do osetnog smanjenja zasićenosti. Uz neophodnu pažnju koristeći solver W.O.D.A. moguće je uspešno rešiti ovu jednačinu.
Jedan od nedostataka ovakvog fizičkog modela jeste da nema održanja mase vazduha, odnosno vazduh može da nestaje i nastaje niotkud kada pritisak dovoljno poraste odnosno padne. Ovo je direktna posledica relacije (18). U najvećem broju slučajeva ovo ne predstavlja problem, međutim treba imati u vidu da se ovako ne mogu modelirati recimo vazdušni džepovi.

3. GRANIČNI USLOVI

3.1 O graničnim uslovima
Rešenje jednačine (20) nije jedinstveno određeno dok se ne zadaju granični uslovi i početni uslov. Pretpostavlja se da hidraulički potencijal h zadovoljava jednačinu (20) unutar nekog ograničenog povezanog skupa koji se naziva domen. U svakoj tački granice domena neophodno je definisati nekakav granični uslov. U našem kodu postoje tri tipa graničnih uslova:
Dirihleov granični uslov, kada se zadaje vrednost hidrauličkog potencijala h;
Nojmanov granični uslov, kada se zadaje proticaj;
Robinov granični uslov, ili granični uslov reke, gde se zadaje nivo u reci (ili jezeru isl) h 0 i koeficijent transfera Ψ = K / d , gde je K nepoznati koeficijent filtracije kolmiranog sloja korita reke, a d je njegova nepoznata debljina. U ovom slučaju ni reka ni kolmirani sloj ne pripadaju modelu, a granica modela je dno kolmiranog sloja. Ni potencijal ni proticaj kroz granicu modela nisu određeni ovim graničnim uslovom, već samo njihova međusobna zavisnost u = Ψ ( h − h 0 ) .
Granični uslov se može zadati samo na granici modela. U Lizzi postoji mogućnost da se granični uslovi zadaju u unutrašnjosti modela što praktično uvodi unutrašnju granicu u model, ali ova mogućnost nije podržano u solveru W.O.D.A. i ne treba je koristiti.

3.2. Uslov na granici modela - potencijal
Prirodna obeležja (šta može da bude u prirodi).
Karakteristike u modelu.
Zadavanje vremenske funkcije.
* Unošenje podataka za potencijal na konturama modela.

3.3. Uslov na granici modela – podzemni doticaj/oticaj
U ovom slučaju zadaje se zapreminski proticaj po jedinici površine u jedinici vremena [m/s].
Prirodna obeležja (šta može da bude u prirodi).
Karakteristike u modelu.
* Unošenje podataka za potencijal na konturama modela.

3.4. Reka (slobodna vodena površina)
Objašenjenje uloge u modelu. Kako radi. Ovo je Robinov granični uslov.
Parametri: podužno prostiranje, kontura korita (profili), kolmiranost, obale, inundacija, nasip (granica razlivanja).
* Unošenje podataka: centralna linija, profili, obalske linije

3.5. Bunar
Pojam bunara u modelu. Kako radi.
Parametri: filter, kolmiranost, isprana zona.
Mogućnosti zadavanja: kao potencijal, kao proticaj u jedinici vremena [ m 3 / s ] . Ograničenja.
Granični uslov u bunaru zadaje se na dnu najnižeg filtra u slučaju cevastog bunara, odnosno na nivou najnižeg otvora za dren u šahtu Reni bunara.
Bunar kao pijezometar.
Unošenje podataka za bunar.

3.6. Početni uslov i eventualne završne napomene o modelu
Da bi razvoj strujne slike kroz vreme bio određen neophodno je zadati rešenje u početnom trenutku t = 0 . Početno rešenje može da bude nivo zadat jedinstveno za ceo model, ili strujna slika koja je prethodno dobijena nekakvim proračunom.
Stacionarno rešenje ne zavisi od početnog rešenja. Međutim, budući da se za rešavanje jednačine podzemnog strujanja koristi iterativni postupak što znači da se nalaze sve tačnija i tačnija približna rešenja dok se ne dostigne zadata tačnost, početno rešenje se koristi i u stacionarnom slučaju kao početna aproksimacija stacionarnog rešenja. Tada je potrebno da početno rešenje bude što je moguće bliže traženom stacionarnom rešenju kako bi se iterativni postupak brže završio ili kako bi se uopšte završio. Kasnije će biti više reči o tome kako doći do dovoljno dobrog početnog rešenja.

4. PRORAČUN

4.1 Mreža
Mreža predstavlja podelu domena na poliedarske ćelije. U slučaju Lizze ove ćelije su uspravne trostrane prizme.
Kako se formira mreža modela u planu i zašto je to važno.
Vertikalna diskretizacija.

4.2 Diskretizacija metodom konačnih zapremina
Zakoni održanja spadaju u tzv. egzaktne fizične zakone. Oni su sasvim tačni i ne uključuju bilo kakvu aproksimaciju. Na primer, u odsustvu izvora i ponora mase, promena mase u određenoj zapremini jednaka je zbiru proticaja kroz njene granice sasvim tačno, do na poslednji atom. Sa druge strane, Darsijev zakon spada u tzv. konstitutivne relacije, zajedno sa recimo Omovim ili Fikovim zakonom. Ovo su približni zakoni.
Kada se konstruiše numerički metod, poželjno je da se diskretizacija izvrši na način koji ne unosi nikakvu grešku u egzaktne fizičke zakone, već se greška aproksimacije pravi samo u konstitutivnim relacijama, koje su ionako približne.
Diskretizacija jednačine podzemnog strujanja (20) metodom konačnih zapremina vrši se tako što se svakoj ćeliji mreže Ω i dodeli jedan stepen slobode hidrauličkog potencijala h i . Zatim se na zapreminski integral jednačine održanja (9) primeni teorema o divergenciji
∫ Ω i ∂ θ ∂ t d Ω = ∫ Ω i ∇ ⋅ u d Ω = ∮ ∂ Ω i n ⋅ u d Γ .	(21)
Budući da je ova ćelija poliedar (u našem slučaju trostrana prizma), poslednji integral se može razbiti na sumu proticaja kroz strane ovog poliedra:
∫ Ω i ∂ θ ∂ t d Ω = ∑ f ∫ f n ⋅ u dS = ∑ f u f ,	(22)
gde je u f isticaj iz Ω i kroz stranu f.
Primenom Njutn-Lajbnicove formule na levu stranu integrala ove jednakosti po vremenu od t n do t n + 1 = t n + Δ t dobija se
θ i n + 1 − θ i n = Δ t ∑ f u f ,	(23)
gde je θ i = ∫ Ω i θ d Ω , a gornji indeks predstavlja vremenski korak.
Primetimo da u ovoj diskretizaciji do sada nije korištena nikakva aproksimacija, tj. ona je za sada potpuno tačna. Budući da je zakon održanja mase već diskretizovan tj. zamenjen sumom proticaja, ova egzaktnost nam omogućava da dobijemo rešenje koje ne sadrži nikakvu grešku u balansu, što nije moguće u slučaju metoda konačnih elemenata. Međutim, moramo napraviti aproksimaciju da bismo predstavili proticaje u f kao kombinacije vrednosti hidrauličkog potencijala u susednim ćelijama.
Najjednostavniji način da se aproksimira proticaj u f = − ∫ f k r ( s ) n T K ∇ h dS je
u f ≈ ∣ f ∣ k r ¯ ∥ K n ∥ h i − h out ∥ x i − x out ∥ ,	(24)
gde je ∣ f ∣ površina strane f, h out je hidraulički potencijal u ćeliji koja se nalazi sa spoljašnje strane strane f u odnosu na ćeliju Ω i , x i i x out su tačke kolokacije u ove dve ćelije tj. tačke za koje se vezuje vrednost hidrauličkog potencijala (u našem slučaju težišta ćelija), a k r ¯ je nekakva srednja vrednost relativne propustljivosti između ove dve ćelije.
Dva su problema vezana za aproksimaciju (24). Prvo, kako izračunati k r ¯ ? Relativna propustljivost se može izračunati pomoću Van Genuchtenovog modela u tačkama gde je poznat hidraulički potencijal, a to je u ćelijama, međutim u strani gde se traži relativna propustljivost tog podatka nema. Standardni pristupi su da se k r ¯ dobije kao neka od standardnih sredina (aritmetička, geometrijska, harmonijska itd) iz vrednosti k r u centrima ćelija, ili da se prvo hidraulički potencijal usrednji nekako i da se onda izračuna k r iz Van Genuchtenove relacije. Svi ovi pristupi daju nezadovoljavajuće i veoma različite rezultate. Krivac za ove razlike je to što se vrednosti k r veoma brzo menjaju u blizini tačke zasićenja. Dosta tačni rezultati se dobijaju kada se uzme integralna sredina
k r ¯ = 1 h out − h i ∫ h i h out k r ( h ) dh .	(25)
Ovaj integral je neophodno izračunati numerički. Solver W.O.D.A. koristi jedan savršeniji metod zasnovan na integralnoj sredini.
Drugi problem vezan za aproksimaciju (24) je pitanje tačnosti aproksimacije izvoda n T K ∇ h konačnom razlikom − ∥ K n ∥ h i − h out ∥ x i − x out ∥ . Naime, n T K ∇ h predstavlja izvod hidrauličkog potencijala u pravcu vektora K n = n T K pomnožen sa intenzitetom ovog vektora. U izotropnom slučaju ovo je prosto izvod u pravcu normale n pomnožen sa koeficijentom filtracije. Sa druge strane, konačna razlika h i − h out ∥ x i − x out ∥ je dobra aproksimacija izvoda u pravcu vektora x i − x out , ali se ovaj vektor ne poklapa uvek sa normalom, ili sa vektorom K n u izotropnom slučaju. Štaviše, kada se mreža profini (svaka ćelija se podeli na više sitnijih ćelija), ovaj problem se ne smanjuje, i dakle ni greška ovakve aproksimacije se neće smanjiti, što znači da je ovakva aproksimacija nekonzistentna.
Ipak, ovakva aproksimacija je upotrebljiva za potrebe proračuna podzemnog strujanja. Osnovni kvalitet ove aproksimacije je to što ona zadovoljava princip maksimuma. To je matimatička teorema koja garantuje da hidraulički potencijal u unutrašnjosti modela ne može da ima lokalni minimum ili maksimum, osim ako ovakav minimum ili maksimum postoji u početnom rešenju. Tačno rešenje jednačine podzemnog strujanja zadovoljava princip maksimuma, ali numerička rešenja dobijena drugim numeričkim metodama to obično ne zadovoljavaju. Činjenica da aproksimacija (24) zadovoljava princip maksimuma znači da će dobijeno rešenje biti u skladu sa fizikom, i da je ono zapravo tačno rešenje koje odgovara koeficijentu filtracije koji nije onaj zadati već varira od ćelije do ćelije onoliko koliko pravac vektora x out − x i odstupa od pravca vektora K n . Sa druge strane, koeficijent filtracije koji je homogen unutar sloja ne modelira prirodu ništa tačnije. Vrednost ovakvog homogenog koeficijenta filtracije se može dobiti merenjem veoma približno, ili se dobija tariranjem modela. Greške koje se na taj tačin prave teško da mogu biti manje od greške ove aproksimacije. I konačno, cilj modeliranja nije tačan koeficijent filtracije, već model čije ponašanje što je moguće vernije oslikava realnost, a diskutabilno je da li se ovaj cilj može bolje dostići ako se koristi tačnija aproksimacija proticaja.
Ipak, tačnija aproksimacija je implementirana u solveru W.O.D.A. To je aproksimacija oblika
u f ≈ M 1 h i − M 2 h out + r ,	(26)
gde čak i u zasićenom slučaju koeficijenti M 1 , M 2 i r zavise ne samo od geometrije i materijalnih karakteristika kao u slučaju aproksimacije (24), već i od hidrauličkog potencijala. Zbog toga je ova aproksimacija nelinearna i znatno složenija, i nećemo ulaziti u njene detalje.
Sadržaj vode θ i koji se pojavljuje u jednačini (23) računa se pomoću Van Genuchten - Mualemove relacije (18) iz hidrauličkog potencijala h i . O ovome će biti više reči u poglavlju 4.3.
* Eksportovanje modela iz Lizze i proračun sa Wodom.

4.3 Nelinearnost
Proračun kreće od početnog rešenja h 0 i zatim se računa h 1 , h 2 i tako dalje, dok se ne dostigne unapred zadato krajnje vreme t=T. Veličina vremenskog koraka Δ t se može zadati ručno ili se može odrediti automatski, o čemu će biti reči u poglavlju 4.4.
Najjednostavnije bi bilo izračunati u f iz poznatih vrednosti potencijala, a zatim iz jednačine (23) direktno izračunati θ i n + 1 i odatle dobiti h n + 1 numeričkim rešavanje Van Genuchtenovog modela. Međutim, ovakav pristup nameće izuzetno strogo ograničenje na veličinu vremenskog koraka (tzv. CFL uslov) i praktično je neprimenljiv na velike modele jer bi proračun predugo trajao. Zbog toga se diskretizacija proticaja vrši implicitno. Jednačina (23) postaje
	(27)
gde je proticaj u f n + 1 izražen preko nepoznatih vrednosti potencijala u okolnim ćelijama na novom vremenskom koraku.
Jednačine oblika (27) čine sistem koji treba rešiti kako bi se dobio hidraulički potencijal na vremenskom koraku n+1. Osim u najjednostavnijem slučaju, ove jednačine su nelinearne iz više razloga:
Zasićenost koja je potrebna da se izračuna sadržaj vode θ i n + 1 zavisi od hidrauličkog potencijala h i n + 1 na nelinearan način određen Van Genuchtenovim modelom.
Relativna propustljivost koja se pojavljuje u izrazima za proticaj zavisi od zasićenosti na nelinearan način.
Ukoliko se koristi nelinearni metod za diskretizaciju proticaja (26) to unosi dodatnu nelinearnost.
Prema tome dobijeni sistem jednačina je linearan jedino u zasićenom slučaju ako se koristi linearna diskretizacija proticaja.
U solveru W.O.D.A. implementirana su tri metoda za rešavanje ovog sistema nelinearnih jednačina: Picarov metod, modifikovani Pikarov ili Celia metod i Njutnov metod. Preporučuje se korišćenje kombinacije drugog i trećeg metoda, ali kako bismo objasnili njihove razlike poći ćemo od Pikarovog metoda.
U Pikarovom metodu vrši se aproksimacija
θ i n + 1 − θ i n ≈ ∂ θ ∂ h ( h i n + 1 − h i n ) .	(28)
Izvod ∂ θ ∂ h se može dobiti iz Van Genuchtenovog modela i iz koeficijenta specifične izdašnosti izdani pod pritiskom S. Jednačina (27) postaje
∂ θ ∂ h ( h i n + 1 − h i n ) = Δ t ∑ f u f n + 1 .	(29)
Sistem ovakvih jednačina se može zapisati u matričnom obliku
A ( h n + 1 ) h n + 1 = r ( h n + 1 ) ,	(30)
gde je h n + 1 vektor koji se sastoji od po jedne vrednosti potencijala za svaku ćeliju mreže, A je matrica koja zavisi od
Nelinearni solveri: Picard, Celia i Newton. Tačnost i broj iteracija nelinearnog solvera. Zavisnost konvergencije od vremenskog koraka.
* Poređenje proračuna sa različitim nelinearnim solverima

4.4 Automatska kontrola vremenskog koraka
Algoritam i parametri automatske kontrole vremenskog koraka
saveStart, saveStep
Pseudo-stacionarni proračun
* Poređenje proračuna sa različitim parametrima automatske kontrole vremenskog koraka

4.5 Linearni solveri
Direktni i iterativni linearni solveri. BiCGSTAB i UMFPACK. Broj operacija i tačnost.
Definicija balansa i njegova uloga u proveri validnosti rešenja.
* Poređenje proračuna sa različitim parametrima linearnog solvera
* Učitavanje rezultata u Lizzu i Paraview

5. KALIBRACIJA

5.1. Uvodni pojmovi, izbor parametara

5.2. Analiza osetljivosti
Veličina modela. Dominantni parametri.
12