SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
SIMULAREA DINAMICII APELOR SUBTERANE
Simularea dinamicii apelor subterane reprezintă finalitatea procesului de
modelare a fenomenului de curgere prin mediul poros – permeabil.
Prin simulare se obţine răspunsul sistemului analizat la o serie de solicitari
Etapele modelării matematice sunt:
Model conceptual;
Model matematic;
Model numeric;
Validarea modelului;
Calibrarea (estimarea parametrilor);
Simularea proceselor analizate.
MODELUL CONCEPTUAL
Modelul conceptual este rezultatul schematizării structurilor acvifere
Schematizarea structurilor acvifere este o reprezentare simplificată a acestora,
referitoare la condiţiile geologice şi structurale ale depozitelor permeabile şi la cauzele
mişcării apei subterane.
Schematizarea este iterativa, finalizându-se în etapa de calibrare a modelului şi
asigură reprezentativitatea modelului
Modelul conceptual este reprezentativ dacă are o comportare similară cu
structura acviferă reală.
Etapele schematizării structurilor acvifere sunt:
Schematizarea spaţială are ca scop identificarea dezvoltării spaţiale a domeniului
modelului (model 1D, 2D, 3D) şi definirea geometriei frontierelor, respectiv a limitelor
modelului.
Schematizarea spaţială presupune:
localizarea şi dezvoltarea acviferelor pe baza secţiunilor geologice şi a
hărţilor structurale.
Dr. ing. Lacramioara Coarna Page 1
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Pe baza secţiunilor geologice se identifică structura acoperişului (în cazul
acviferelor sub presiune) şi a culcuşului impermeabil/semipermeabil pentru stratul
acvifer.
Principalele hărţi structurale utilizate în descrierea modelului conceptual sunt:
hărţile cu izohipse, izobate la culcuşul şi la coperişul stratului acvifer şi cu izopahite.
stabilirea relaţiilor hidraulice dintre acvifere şi a relaţiilor hidraulice ale acviferelor cu apele de suprafaţă.
Schematizarea parametrică urmăreşte estimarea proprietăţilor acvifere ale
formaţiunilor permeabile.
Prin schematizarea parametrică se urmăreşte:
distribuţia spaţială a proprietăţilor acvifere ale formaţiunilor litologice: porozitate efectivă (ne), conductivitate hidraulică (K), transmisivitate (T), coeficient de difuzivitate hidraulică (a), coeficient de inmagazinare (S), coeficient de realimentare (w)
variabiliatea spaţială a parametrilor (datorată condiţiilor de formare a acviferelor)
pe baza căreia se identifică tipul formaţiunilor:
- formaţiuni omogene/neomogene: prin analiza variabilităţii spaţiale a compoziţiei
litologice;
- formaţiuni izotrope / anizotrope / ortotrope: prin analiza variabilităţii spaţiale a
permeabilităţii (pe diferite direcţii).
În cadrul schematizării parametrice se analizează influenţa scării la care se
realizează modelul matematic.
Pentru modelul regional, prin schematizare se reduce variabilitatea
parametrilor, putând considera mediul omogen şi izotrop. Astfel, porozitatea şi
conductivitatea hidraulică au valori constante egale cu valoarea medie a valorilor
disponibile iar coeficientul de înmagazinare elastică influenţează câmpul piezometric al
acviferelor.
Pentru modelul local, zona modelată este divizată în subzone omogene si
izotrope/ortotrope, prin distorsiune parametrică şi distorsiune geometrică.
Dr. ing. Lacramioara Coarna Page 2
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Schematizarea hidrodinamică are ca scop precizarea contextului energetic al
curgerii prin stabilirea:
regimului de curgere, care poate fi laminar sau turbulent;
tipului de curgereÎn funcţie de evoluţia în timp a debitului, curgerea poate fi:
- staţionară (∂Q∂ t
=0);
- nestaţionară (∂Q∂ t
≠0).
În funcţie de condiţiile de alimentare / descărcare pe verticală din infiltrare de la
suprafaţa terenului sau prin drenanţă din / spre acviferele învecinate, curgerea poate fi
conservativă sau neconservativă.
Condiţii iniţiale care indică distribuţia spaţială a H pe domeniul de curgere la
un moment de început al procesului: H = H(x,y,z,t = 0)
o Condiţii de frontieră / de margine / pe limita domeniului. Acestea pot fi de
de trei tipuri: Dirichlet, Neumann şi Cauchy.
o Condiţii de tip Dirichlet (H impus):
Sarcina piezometrică (H) impusă pe frontieră este independentă de condiţiile de
curgere din acvifer şi poate fi constantă sau variabilă. Condiţia de margine de tip
Dirichlet se impune în următoarele situaţii:
- Sarcină piezometrică constantă: H = const. pentru:
- limita acvifer – apă de suprafaţă (lac, mare): H = Zapă
- zona de aflorare a unui strat acvifer: H = z+hp ()
- linia de izvoare prin care se descarcă un acvifer
- Sarcină piezometrică variabilă: H ǂ const. pentru:
- limita acvifer – râu: H = Hrâu
o Condiţii de tip Neumann (q = K(dH/dn) impus)
Debitul impus (q) pe frontieră, exprimat prin prima derivată a sarcinii
piezometrice pe normala la frontiera domeniului, poate fi nul sau nenul: Condiţia de margine de tip Neumann se impune în următoarele situaţii:
Dr. ing. Lacramioara Coarna Page 3
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
- Debit impus nul: q = qi = 0, pentru limite impermeabile sau paralele cu liniile de
curent (limită/frontieră tangentă la direcţia de curgere, paralelă cu liniile de curent).
Acest tip de condiţii se aplică în următoarele cazuri:
- limite impermeabile (culcus / coperis impermeabil);
- pereţi impermeabili ai rezervoarelor de apă îngropate/semiîngropate
- Debit impus nenul: q = qi ǂ 0, pentru limite traversate de un debit independent de
condiţiile de curgere
Acest tip de condiţii se aplică în următoarele cazuri:
- aporturi de apă în acvifer: infiltraţii, injectii
- pierderi de apă din acvifer: captări de apă: foraje, drenuri
- interfaţa a două formaţiuni de conductivitate net diferită: acvifere aluvionare
cu conductivitate mare plasate pe formaţiuni permeabile cu conductivitate redusă
(K2>>K1). Astfel, pentru acviferul aluvionar cu conductivitate mare K2: frontieră de tip
debit impus constant (q = qi = const.) iar pentru acviferul din formaţiuni cu
conductivitate redusă K1: frontieră de tip sarcină piezometrică impusă constantă (H =
Hi = const.).
o Condiţii de tip Cauchy (debit dependent de sarcina piezometrică: Q = Q (H)
impus) Acest tip de condiţii se aplică în următoarele situaţii:
- pentru un strat semipermeabil, de grosime M’ şi conductivitate K’, care separă
un acvifer cu nivel liber (ANL) de un acvifer inferior, sub presiune (ASP), caracterizat
prin sarcină piezometrică mai mare (HASP > HANL)
q=K ' H ASP−H ANL
M '
- la contactul acviferului cu apele de suprafaţă apare un strat subţire, de
grosime M* şi conductivitate redusă K*, datorat colmatării acviferului cu sedimente
purtate de râuri în timpul viiturilor
q=K ¿ H râu−H acv
M ¿
MODELUL MATEMATIC
Dr. ing. Lacramioara Coarna Page 4
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Modelul matematic constă în caracterizarea prin ecuaţii matematice a dinamicii
apelor subterane în domeniul modelului ale cărui proprietăţi sunt descrise prin modelul
conceptual. Soluţia unică a ecuaţiei de curgere este obţinută prin impunerea condiţiilor
iniţiale şi de margine.
Dinamica apei subterane prin mediul poros – permeabil (mediu continuu
echivalent) este descrisă de ecuaţia de difuzivitate.
¿ ( K⃗ gradh )=S ∂h∂ t
în care, S – coeficientul de înmagazinare (-)
Pentru acviferul cu nivel liber, coeficientul de înmagazinare este dat de
porozitatea efectivă:
S=ne
Pentru acviferul sub presiune, coeficientul de înmagazinare este exprimat în
funcţie de greutatea specifică a apei (γa )şi de coeficientul capacităţii elastice a
complexului apă – rocă β=[ (1−n ) ∙ βS+n∙ βa ]:S=γ a∙ [ (1−n ) ∙ βS+n ∙ βa ]..........
Ecuaţia de difuzivitate este obţinută din ecuaţia de continuitate şi din legea
Darcy.
Ecuaţia de continuitate a unui fluid de densitate constantă printr-un mediu
continuu având porozitatea n, este:
−¿ (ρ ∙ q⃗ )= ∂∂ t
(n ∙ ρ )
−∂∂x ( ρ∙qx )− ∂
∂ y (ρ ∙q y)−∂∂z ( ρ∙qz )= ∂
∂t(n ∙ ρ )
Mişcarea fluidului prin mediu poros caracterizat prin conductivitatea K (K x , K y ,K z )
este descrisă de legea Darcy
q⃗=−KgradH
Componentele debitului specific pe direcţiile unui sistem tridimensional Oxyz
sunt:
Dr. ing. Lacramioara Coarna Page 5
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
qx=−K x∂ H∂ x , q y=−K y
∂H∂ y , qz=−K z
∂ H∂ z
MODELUL NUMERIC
În condiţiile în care modelul analitic simplificat nu descrie cu fidelitate problema
de rezolvat, ecuaţiile diferenţiale pot fi aproximate numeric. În consecinţă, variabilele
continue sunt înlocuite cu variabile discrete care sunt definite pentru celulele sau
nodurile reţelei ataşate modelului.
Astfel, ecuaţiile diferenţiale continue care definesc sarcina piezometrică sau
concentraţia soluţiei oriunde în sistem sunt înlocuite printr-un număr finit de ecuaţii
algebrice care definesc sarcinile piezometrice sau concentraţiile în puncte specifice.
Sistemul de ecuaţii algebrice este rezolvat matricial.
În rezolvarea ecuaţiilor de curgere în mediul subteran sunt acceptate două clase
de metode numerice: metoda diferenţelor finite şi metoda elementului finit.
Fiecare din cele două clase de metode numerice prezintă o mare varietate de
subclase şi implementări alternarative.
Ambele metode necesită ca aria modelului (de interes) să fie subdivizată printr-o
reţea într-un număr de arii mai mici (celule sau elemente) care sunt asociate cu
punctele nodale (fie centrate, fie în colţuri).
Pe lângă cele două metode frecvent utilizate, mai există două metode al căror
principal avantaj este că pentru regiuni omogene pot furniza soluţii precise fără
discretizare:
- metoda ecuaţiilor integrale de margine:
- metoda elementelor analitice.
Astfel, dacă heterogeneitatea sistemului poate fi adecvat reprezentată folosind
doar câteva elemente foarte extinse, aceste metode pot fi foarte eficiente.
Pentru sisteme heterogene care trebuiesc descrise printr-un număr mare de
elemente, sunt recomandate metoda diferenţelor finite (MDF) şi metoda elementului finit
(MEF).
Până acum, pentru descrierea dinamicii apelor subterane au fost frecvent
utilizate MDF şi MEF.
Dr. ing. Lacramioara Coarna Page 6
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Metoda diferenţelor finite (MDF) aproximează prima derivată din ecuaţia
diferenţială cu derivare parţiale prin diferenţe de coeficienţi (diferenţa dintre valorile
variabilei independente în noduri adiacente, raportate la distanţa dintre două noduri şi la
două momente succesive de timp, referitor la durata unui pas de timp.
Metoda elementelor finite (MEF) utilizează funcţii fictive ale variabilei
dependente şi ale parametrilor pentru a evalua o formulare integrală echivalentă a
ecuaţiilor diferenţiale cu derivate parţiale.
În ambele abordări numerice, discretizarea dimensiunilor de spaţiu şi timp
permite obţinerea de valori continue pe limitele domeniului si reduce aflarea soluţiei
ecuaţiei diferenţiale cu derivate parţiale la aflarea soluţiei unui set de ecuaţii algebrice.
Aceste ecuaţii pot fi rezolvate folosind metode iterative sau matriciale.
Fiecare metodă are avantaje şi dezavantaje.
În general, MDF este mai simplă conceptual şi matematic şi este mai uşor de
programat. Acestea sunt manipulate, în mod specific, cu o reţea rectangulară, relativ
simplă.
MEF necesită, în general, utilizarea unor metode matematice mai sofisticate
pentru câteva probleme care pot fi mai exacte numeric decât MDF standard. Un avantaj
major al MEF este flexibilitatea reţelei de elemente finite care permite o cât mai bună
(exactă) aproximare spaţială a frontierelor neregulate ale acviferului şi/sau a
parametrilor din zonele acviferului, când sunt luate în considerare.
Cu toate acestea construcţia şi descrierea setului datelor de intrare sunt mai
dificile pentru o reţea de elemente finite neregulate decât pentru o reţea de elemente
finite rectagulare regulate.
Astfel, se recomandă utilizarea unui model preprocesor care permite generarea
reţelei, prezentarea schematică a numărului de noduri şi elemente ale reţelei şi
atribuirea coordonatele spaţiale ale fiecărui nod.
În figura de mai jos este prezentat în spaţiu 2D un acvifer care are limite
impermeabile şi prezintă o captare printr-un foraj şi reţelele de discretizare în cazul
MDF şi MEF. De asemenea, în figură este indicată ajustarea discretizării reţelei (reţea
mai fină) în zona captării.
Dr. ing. Lacramioara Coarna Page 7
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Reţeaua rectangulară aproximează graniţele acviferului (modelului) în mod
treptat, obţinându-se câteva noduri sau celule în afara limitei modelului. Elementele
reţelei triunghiulare, utilizate în MEF pot urmări cu mare fidelitate („îmbrăca”) graniţa
limitei modelului folosind un număr minim de noduri.
Metoda diferenţelor finiteModelul conceptual care descrie un mediu continuu, pentru care se presupune
că se cunosc proprietăţile acviferului, este transformat într-un mediu discret obţinut
printr-o reţea în care se cunosc proprietăţile acviferului în fiecare punct şi pe limitele
domeniului.
MDF presupune discretizarea domeniului de curgere printr-o reţea patrulateră.
Fiecărei reţele îi sunt atribuite nodurile în care se obţin valorile necunoscute ale sarcinii
piezometrice, prin rezolvarea ecuaţiei de curgere. De asemenea, în nodurile reţelei sunt
atribuite valorile cunoscute ale unor parametrii, cum ar fi: conductivitatea hiraulică (K),
transmisivitatea (T), capacitatea de înmagazinare (S). Nodurile reţelei ataşate
domeniului se pot afla, fie în centrul celulelor („block – centered”), fie la intersecţia
celulelor („mesh – centered„).
Alegerea uneia din cele două variante ale reţelei depinde condiţiile de margine
impuse modelului:
- reţeaua de tip „block – centered” este des utilizată cand se impun condiţii de
tip Neumann (flux impus);
- reţeaua de tip „mesh – centered” este des utilizată cand se impun condiţii de
tip Dirichlet (sarcină piezometrică impusă).
Dr. ing. Lacramioara Coarna Page 8
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Reţeaua poate fi regulată, având liniile şi coloanele perpendiculare. Reţeaua
poate fi pătratică dacă Dx = Dy. Adesea, este convenabil ca mărimea rândurilor şi
coloanelor reţelei să varieze, astfel că numărul de noduri este mai mare într-o anumită
zonă a reţelei (în zona forajelor, în zona cu variaţii mari ale lui h, faţă de vecinătăţi).
Mărimea variaţiei pasului celulelor Dx, sau Dy, de la o coloană sau un rând la cea / cel
adiacent este de preferat să fie de 30% - 50%.
Ecuaţia diferenţială a curgerii apelor subterane conţine termeni care reprezintă
derivatele în spaţiu şi timp ale variabilelor continue. MDF constă în aproximarea acestor
derivate (pantele curbelor) prin variaţii liniare discrete pe intervale discrete de spaţiu şi
timp. Dacă aceste intervale sunt suficient de mici, toate creşterile liniare vor reprezenta
o bună aproximare a suprafeţei curbilinii reale (hidrografului).
Considerând un foraj de observaţie într-un acvifer sub presiune, derivata sarcinii
piezometrice pe direcţia Ox, poate fi scrisă ca:
- diferenţe finite progresive:
dhdx
=h ( x+∆ x )−h ( x )
∆ x+O (∆ x2 )
- diferenţe finite regresive:
dhdx
=h ( x )−h ( x−∆ x )
∆ x+O (∆x2 )
Diferenţele finite se pot exprima, din dezvoltarea în serie Taylor ca:
h ( x+∆ x )=h ( x )+ dhdx
∆ x+ 12!
d2hd x2
∆ x+O (∆ x3 )
h ( x−∆x )=h (x )−dhdx
∆ x+ 12!
d2hd x2
∆ x−O (∆ x3 )
Adunând cele două relaţii obţinem:
h ( x+∆ x )+h ( x−∆ x )=2 ∙ h (x )+ d2hd x2
∆ x
Astfel, derivata de ordin doi a sarcinii piezometrice pe direcţia Ox se poate scrie ca:
d2hd x2
=h (x+∆x )+h ( x−∆ x )−2∙ h ( x )
∆ x2
Prin analogie, derivata de ordin doi a sarcinii piezometrice pe direcţia Oy, este:
Dr. ing. Lacramioara Coarna Page 9
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
d2hd y2
=h ( y+∆ y )+h ( y−∆ y )−2 ∙ h ( y )
∆ y2
În cazul unei reţele pătratice, pentru care se exprimă valorile sarcinii
piezometrice în noduri, aproximarea variaţiei sarcinii piezometrice pe direcţiile Ox şi Oy
se exprimă:
- prin diferenţe finite progresive
dhdx
=hi+1 , j−hi , j
∆ x+O (∆ x2 )
dhdy
=hi , j+1−hi , j
∆ y+O (∆ y2 )
- prin diferenţe finite regresive
dhdx
=hi , j−hi−1 , j
∆x+O (∆x2 )
dhdy
=hi , j−hi , j−1
∆ y+O (∆ y2 )
Derivatele de ordin doi ale sarcinii piezometrice h, pe direcţiile Ox şi Oy sunt:
d2hd x2
=hi+1 , j+hi−1 , j−2 ∙h i , j
∆ x2
d2hdy
=hi , j+1+hi , j−1−2∙ hi , j
∆ y2
Dr. ing. Lacramioara Coarna Page 10Δx
Δy
hi+1,jhi,j
hi,j+1
hi-1,j
hi,j-1
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Curgerea apelor subterane printr-un acvifer omogen şi izotrop, în regim staţionar,
este descrisă de ecuaţia:
d2hd x2
+ d2hd y2
=0
În condiţiile unei reţele pătratice (Δx = Δy), ecuaţia diferenţială a curgerii
subterane, poate fi aproximată prin:
hi+1 , j+hi−1 , j−2 ∙ hi , j
∆ x2+hi , j+1+hi , j−1−2 ∙ hi , j
∆ y2=0
Adică:
hi+1 , j+hi−1 , j+hi , j+1+hi , j−1−4 ∙ hi , j
∆x2=0
Sarcina piezometrică în nodul (i,j) al reţelei pătratice ataşată domeniului de
curgere al apelor subterane poate fi exprimată în funcţie de sarcinile piezometrice în
nodurile din imediata apropiere ale nodului (i,j):
hi , j=hi+1 , j+hi−1, j+hi , j+1+hi , j−1
4
Aproximarea derivatelor parţiale prin diferenţe finite este afectată de erori de
trunchiere şi erori de rotunjire. Erorile de trunchiere sunt datorate dezvoltării în serie
Taylor şi sunt cu atât mai mici cu cât reţeaua este mai fină (celulele sunt mai mici).
Erorile de rotunjire apar la rezolvarea sistemului de ecuaţii ataşat reţelei, sistem
de N ecuaţii (N este numărul de noduri ale reţelei) în care sarcina piezomterică în
fiecare nod al reţelei se calculează în funcţie de sarcinile piezometrice din nodurile
aflate în imediata vecinătate.
Metoda celulelorCurgerea apelor subterane printr-un acvifer neomogen şi anizotrop, în regim
staţionar, se analizează utilizând metoda celulelor. Prin această metodă se exprimă
necunoscuta ecuaţiei de curgere (sarcina piezometrică) în centrul fiecărei celule ale
reţelei de latură a, în funcţie de valorile lui H în nodurile învecinate.
Presupunând că, pentru fiecare celulă T = const.,
- debitul ce traversează fiecare celulă se exprimă ca:
Dr. ing. Lacramioara Coarna Page 11
hij
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Q=K ∙ I ∙ M ∙a→Q=T ∙ I ∙ a→Q=T ∙ ∆ Ha
∙a
- debitul din celula centrală este:
Q=q ∙a2
Astfel, pentru celulele din figura xxx, avem:
−HC−HV
a∙a∙T CV+
HE−HC
a∙a ∙T CE−
HC−H N
a∙a ∙T CN+
H S−HC
a∙a ∙ T CS=q ∙ a2
−HC (TCV +T CE+T CN+TCS )+T CV ∙HV +TCE ∙H E+T CN ∙H N+T CS ∙ HS=q ∙a2
C ii ∙HC−∑i
c ij ∙ H j=−qi ∙ a2
unde:
C ii=TCV +TCE+T CN+T CS
c={c ii ,daca i= jc ij , daca i≠ j – coeficient de transfer (de legătură, de influenţă), exprimată pe baza
transmisivităţii de pasaj.
Determinarea transmisivităţii de pasaj se face pe baza ecuaţiei de curgere între
două celule învecinate (spre exemplu celulele V şi C din figura XXX)
Dr. ing. Lacramioara Coarna Page 12
HCHV HE
HS
HN
HC
TCTV
HV HC
H*
HV
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
−T CV ∙ a ∙HC−HV
a =−TV ∙ a ∙H ¿−HV
a2
→T CV ∙ ( HC−HV )=2T V ∙ (H ¿−HV )
T CV
2∙ TV∙ (HC−HV )=(H ¿−HV )
−T CV ∙ a ∙HC−HV
a =−TC ∙ a ∙HC−H ¿
a2
→T CV ∙ (HC−HV )=2T C ∙ (HC−H ¿ )
T CV
2∙ TC∙ (HC−HV )=(HC−H ¿ )
Însumând relaţiile de mai sus obţinem:
( T CV
2∙ T V+
T CV
2 ∙T C) ∙ (H C−HV )=(HC−HV )→
TCV
2∙( 1TV
+ 1T C )=1→T CV=2
TV T C
T V+T C→T CV=
2
( 1T V+ 1T C )
Astfel, se obţine că:
- transmisivitatea de pasaj este media armonică a transmisivităţii din ochiurile
învecinate;
- transmisivitatea de pasaj este independentă de sarcina piezometrică a
fiecărei celule; considerând că H = constant pentru fiecare celulă a reţelei.
Metoda celulelor nu asigură continuitatea sarcinii piezometrice, dar asigură
continuitatea debitelor pe o latură, dacă se consideră debitul specific constant pe latura
care separă două celule învecinate
În condiţiile unei reţele pătratice ataşată domeniului de curgere şi utilizând metoda
celulelor, rezolvarea ecuaţiei diferenţiale care descrie curgerea apei subterane se
reduce la rezolvarea ecuaţiei matriciale:
[C ] {H }= {Q }
în care:
[C ] - matricea de permeabilitate (de legătură) care are următoarele
caracteristici:
- este o matrice simetrică C ij=C ji
Dr. ing. Lacramioara Coarna Page 13
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
- este o matrice tip bandă (singuii termeni nenuli corespund legăturii unui nod cu
nodurile învecinate)
{Q } – vectorul debit distribuit în ochiuri
Ecuaţia matricială se rezolvă pe baza unui sistem de ecuaţii, care este univoc
determinat considerând condiţiile de margine (contur / frontiere):
- Condiţii tip Dirichlet:
HC=H 0
- Condiţii tip Neumann:
q=0 - limita impermeabilă
q=q impus
Dacă se cunoaşte q, trebuiesc introduse ochiuri fictive, caracterizate prin Hf şi Tf.
Q=−T f ∙ a ∙HC−H f
a→Q=−T f ∙ (HC−H f )
Transmisivitatea de pasaj a unei celule fictive este:
T f=2
( 1T c+1T c )
→T f =22T C
→T f =TC
Astfel, debitul impus se exprimă prin:
Q=−T f ∙ (HC−H f )→Q=T C ∙ (H f−HC )
Prin metoda diferenţelor finite:
- mediul continuu se transformă într-un mediu discret;
- soluţia sistemului de ecuaţii Laplace (care descrie dinamica apei subterane)
se aproximează cu un sistem de ecuaţii algebrice;
- în cazul curgerii staţionare, matricea de permeabilitate este simetrică şi de tip
bandă iar coeficienţii acestei matrici sunt funcţie doar de transmisivitatea atribuită
fiecărei celule.
Reţeaua ataşată domeniului de curgere are mii de celule/noduri. Rezolvarea
ecuaţiei de curgere presupune obtinerea unei soluţii a acesteia în fiecare celulă / nod.
Aflarea soluţiei se bazează pe rezolvarea ecuaţiei preupunând o valoare iniţială
estimată a H în fiecare celulă / nod al reţelei. Prin condiţia de margine de tip Dirichlet (H
Dr. ing. Lacramioara Coarna Page 14
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
= Hi), H este fixată în nodurile de pe limita domeniului. Pentru nodurile în care nu se
impune condiţia de tip Dirichlet şi în nodurile interioare domeniului, H nu este fixată.
Ecuaţia dieferenţelor finite este rezolvată prin metode iterative. Pe baza valorilor
H fixate şi a estimărilor iniţiale, ecuţia de curgere în MDF se rezolvă pentru fiecare nod
al reţelei pe baza valorilor lui H din cel 4 noduri aflate în imediata vecinătate. Valoarea
sarcinii piezometrice este recalculată până ce diferenţa dintre valoarea estimată iniţial şi
valoarea calculată iterativ este minimă. Procesul iterativ se repetă până ce diferenţa
maximă dintre valoarea H calculată la doi paşi succesivi este cea dată prin criteriul de
convergenţă. Astfel, când soluţia este convergentă, ecuaţia este rezolvată. Cu cât
valoarea dată de criteriul de convergenţă este mai mică cu atât numărul de iteraţii este
mai mare şi timpul necesar atingerii soluţiei este mia mare.
Metode iterative Metoda Gauss – Seidel este metoda iterativă prin care valoarea lui h ij la un pas
al iteraţii se calculează pe baza sarcinii h din două noduri adiacente, calculate la acelaşi
pas al iteraţiei şi ale lui h din alte două noduri adiacente, calculate prin iteraţia
anterioară.
hi , jm+1=1
4 (hi−1 , jm+1 +h i , j−1
m+1 +hi+1 , jm +h1 , j+1
m )
în care: m –1: iteraţia anterioară, m+1: iteraţia curentă.
Metoda suprarelaxării succesiveAceastă metodă iterativă este o variantă a metodei Gauss – Seidel prin care
diferenţa valorilor lui h, calculate într-un nod prin două iteraţii succesive este cunoscut
ca reziduu. În timpul fiecărei iteraţii prin metoda Gauss – Seidel, reziduul se micşorează
până ce se atinge criteriul de convergenţă. La doi paşi succesivi ai metodei de
suprarelaxare se atribuie un factor care ajută la creşterea ratei de convergenţă. Acest
factor de suprarelaxare este cuprins între 1 şi 2, o valoare acceptabilă fiind determinată
prin calculul erorilor.
Astfel,
hi , jm+1=h i, j
m ++ f (hi , jm+1−h1 , j
m )
Dr. ing. Lacramioara Coarna Page 15
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
în care: (hi , jm+1−h1 , j
m ) – reziduul.
Tehnica matricialăPrin MDF şi MEF aproximarea soluţiei ecuaţiei de curgere conduce la un sistem
de ecuaţii algebrice format din ecuaţii pentru fiecare nod al reţelei ataşate domeniului de
curgere.
Sistemul de ecuaţii algebrice poate fi rezolvat numeric prin două metode: metoda
directă şi metoda iterativă.
în metoda directă, o secvenţă a operaţiilor este calculată doar odată cu
rezolvarea ecuaţiei matriciale, asigurând o soluţie exactă, dar afectată de erori de
rotunjire.
Prin metoda iterativă se ajunge la soluţie printr-un proces de aproximări
succesive. Aceasta implică o estimare iniţială a soluţiei şi apoi îmbunătăţirea acestei
estimări printr-un proces iterativ până ce este atins criteriul erorii (obţinerea unei erori
mici a rezultatului).
În ambele metode, trebuie atinsă convergenţa şi rata (viteza) de convergenţă.
Metodele directePrincipalele metode directe de rezolvare a sistemului de ecuaţii matriciale, la
care se reduce rezolvarea printr-o metodă numerică a ecuaţiei de curgere sunt:
- metoda determinanţilor;
- metoda eliminării succesive a necunoscutelor;
- metoda inversării matricilor.
Metodele directe prezintă două dezavantaje principale:
- necesită o resursă mare a calculatorului, memorie mare şi timp mare de
calcul pentru probleme complexe. Matricile de legătură conţin multe valori nule, ceea ce
reduce efortul de calcul.
Erorile de rotunjire se cumulează pentru anumite tipuri de matrici, datorită unui
număr mare de operaţii care conduc la rezolvarea acestora.
Metodele iterative
Dr. ing. Lacramioara Coarna Page 16
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
Prin aceste metode se evită calculul laborios al problemelor cu multe
necunoscute, care presupun matrice de dimensiuni mari.
Dintre numeroasele scheme dezvoltate, cele mai cunoscute sunt: metoda supra-
relaxării succesive, metoda explicită, metoda implicită.
Deoarece metodele iterative încep cu o estimare iniţială a soluţiei, eficeinţa
metodei depinde, cumva, de această estimare iniţială. Pentru accelerarea procesului
iterativ, se utilizează factori de accelerare şi relaxare. Din nefericire, definirea celor mai
bune valori ale acestor factori este, de obicei, o problemă de dependenţă. În plus,
abordarea iterativă necesită ca o eroare a toleranţei să fie specificată pentru finalizarea
procesului iterativ. O valoare optimă pentru toleranţă, care este utilizată pentru a evalua
când calculul iterativ converge către o soluţie, este o problemă de dependenţă.
Dacă această toleranţă este prea mare – procesul iterativ se poate opri înaintea
atingerii unei acurateţe numerice adecvate.
Dacă toleranţă este prea mică – procesul iterativ poate consuma resurse de
calcul foarte mari în atingerea preciziei numerice care poate fi de câteva ordine de
mărime mai mici decât precizia din câmpul datelor, sau procesul iterativ poate chiar să
rateze convergenţă.
Dintre metodele iterative cea mai cunoscută este metoda relaxării, care
presupune introducerea unei soluţii arbitrare {H 0 }.Astfel, ecuaţia de curgere devine:
[C ] {H 0 }−{Q }={R1 }
Unde {R1 } – restul corespunzător primului pas al iteraţiei.
La pasul următor:
[C ] {H 1 }− {Q }={R2 }
Astfel, se obţine o succesiune de soluţii {H }1 , {H }2 ,…, {H }n care converg către
soluţia exactă dacă şi numai dacă termenii diagonali ai matricii sunt pozitivi.
Se descompune matricea în două matrici:
[C ]=[D ]+ [ A ]
- [ D ] este o matrice diagonală
Dij={c ij , daca i= j0 , daca i≠ j
Dr. ing. Lacramioara Coarna Page 17
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
- [ A ] este o matrice cu termeni nuli pe diagonală:
Aij={0 , daca i= jc ij , daca i ≠ j
Ecuaţia de curgere, se poate scrie în formă matricială:
[ D ] {H 1}+ [ A ] {H 0 }={Q }
{H 1 }+[ D ]−1 [ A ] {H0 }=[ D ]−1 {Q }
{H 1 }=[D ]−1 ( {Q }−[ A ] {H 0 })
{H 2 }=[D ]−1 ( {Q }−[ A ] {H 1})Astfel, se obţin soluţiile ecuaţiei de curgere, prin iteraţii succesive, până când:
‖{H }n−{H }n−1‖≤|0|
REGIM DE CURGERE NESTAŢIONARÎn cazul curgerii nestaţionare, panta hidrografului (reprezentarea grafică a
nivelului piezometric în timp) în orice punct reprezintă derivata sarcinii piezometrice în
raport cu timpul şi poate fi aproximată ca:
∂h∂ t
=∆h∆ t
fie prin diferenţe finite progresive (referitor la sarcina piezometrică de la un moment
ulterior):
( ∂h∂ t )n∆ t
=hn+1−hn
∆ t
sau prin diferenţe finite regresive (referitor la sarcina piezometrică de la un moment
anterior):
( ∂h∂ t )n∆ t
=hn−hn−1
∆ t
Ecuaţia de curgere a apei subterane în regim nestaţionar, neconservativ
¿ (TgradH )+Q=S dHdt
poate fi scrisă în formă matricială ca:
Dr. ing. Lacramioara Coarna Page 18
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
[S ]{dHdt }+ [ M ] {H }={Q }
Pentru aproximarea ecuaţiei de curgere a apei subterane în regim nestaţionar
prin metoda diferenţelor finite, se consideră sarcina piezometrică în cinci noduri
învecinate şi la două momente de timp.
Analiza regimului nestaţionar de curgere al apelor subterane se poate face prin
două metode:
- metoda explicită (a) prin care se exprimă derivata spaţială a sarcinii
piezometrice la momentul n (valorile din noduri sunt cunoscute), iar derivata în raport cu
timpul prin diferenţe finite progresive faţă de valorile necunoscute ale sarcinii
piezometrice la momentul ulterior n+1.
Ecuaţiile diferenţiale finite explicite conduc la soluţii simple şi precise (de
încredere).
Aceste metode sunt asociate criteriilor de stabilitate. Dacă pasul de timp este
foarte mare, erorile numerice mici se pot propaga în erori numerice mari în stadiile de
calcul ulterioare.
- metoda implicită (b) prin care se exprimă derivata în raport cu timpul prin
diferenţe finite regresive ale sarcinii piezometrice la momentul n, faţă de valorile
Dr. ing. Lacramioara Coarna Page 19
b)a)
hi,j, n-1
hi,j-1,n
hi,j+1,n
hi+1,j,nhi-1,j,n hi,j,nhi,j+1,n
hi,j-1,n
hi-1,j,n hi+1,j,nhi,j,n
hi,j,n+1
timp
tn-1
tn+1
tn
i
j
SIMULAREA DINAMICII APELOR SUBTERANE - Note de curs
cunoscute ale sarcinii piezometrice la momentul n-1. Aceste valori sunt cunoscute fie
din condiţiile iniţiale specificate pentru prima treaptă de timp sau din soluţiile obţinute
pentru următorii paşi de timp. Derivatele spaţiale ale sarcinii piezomterice se
aproximează pentru momentul n, la care toate valorile sunt necunoscute. Astfel, pentru
fiecare nod al reţelei se scrie o ecuaţie diferenţială care are 5 necunoscute. Această
ecuaţie nu poate fi rezolvată în mod direct. Pentru reţeaua cu N noduri ataşată
modelului se obţine un sistem de N ecuaţii cu N necunoscute. Un astfel de sistem de
ecuaţii simultane împreună cu condiţiile de margine specificate (impuse) pot fi rezolvate
implicit. Cu toate că soluţiile implicite sunt mai complicate, ele au avantajul de a fi, în
general, necondiţionat stabile. Acestea implică faptul că soluţia va fi obţinută chiar dacă
nu este necesar să se exprime derivata care este aproximată. Soluţia obţinută prin
metoda implicită este mai precisă, dacă pasul de timp este mare în raport cu rata de
variaţie a sarcinii piezometrice.
Această metodă implicită este folosită în aproximarea ecuaţiilor de curgere
pentru cele mai accesibile modele de curgere.
Dr. ing. Lacramioara Coarna Page 20