3.Modelarea matematică a curgerii în acvifere

3.1 Definirea şi clasificarea acviferelor

Acviferul este un depozit de roci permeabile capabile să înmagazineze cantităţi importante de apă şi să permită mişcarea acestora prin ele.

Acviferul poate fi ecranat la partea inferioară de un strat impermeabil iar suprafaţa liberă işi poate modifica poziţia în timp. Acest tip de acvifer se numeşte acvifer cu suprafaţă liberă. Dacă acviferul se găseşte între două strate impermeabile el se numeşte acvifer captiv sau sub presiune (fig. 3.1)

Acviferele pot fi mărginite de un strat semi-impermeabil numit acvitard sau de un strat perfect impermeabil denumit acvifug.

Figura 3.1: Acviferul cu suprafaţă liberă şi acviferul sub presiune

Acviferele cu suprafţă liberă, denumite şi acvifere freatice, se reîncarcă de obicei direct de la suprafaţă, din apa de precipitaţie sau din alte mase de apă (râuri sau lacuri) conectate hidraulic cu acviferul.

Acviferul freatic (1) sau pânza freatică este mărginită de un strat impermeabil la bază (2) şi de suprafaţa liberă (3) ( fig. 3.2.).

Nivelul pânzei freatice poate fi vizualizat (măsurat) prin realizarea unui foraj în care apa se va ridica până la nivelul suprafeţei libere.

Figura 3.2: Acviferul cu suprafaţă liberă

Un acvifer sub presiune (confinat) este limitat deasupra şi dedesupt prin formaţiuni impermeabile. Într-un puţ care străbate un acvifer sub presiune, nivelul apei se ridică deasupra stratului impermeabil superior (coperiş). Acest nivel indică sarcina piezometrică din centrul puţului. Nivelul apei dintr-un număr infinit de astfel de puţuri de observaţie reprezintă o suprafaţă imaginară, numită suprafaţă piezometrică (fig. 3.3.).

Figura 3.3 : Acviferul sub presiune

Suprafaţa liberă în cazul acviferelor freatice şi suprafaţa piezometrică în cazul acviferelor sub presiune pot fi suprafeţe oarecare în spaţiu.

Poziţia suprafeţei libere sau a suprafeţei piezometrice se poate stabili faţă de un nivel de referinţă şi poate fi descrisă matematic printr-o funcţie H(x,y).

Dacă pe acestă suprafaţă se trasează curbe H(x,y)=C=constant, pentru fiecare valoare a constantei C se obţine o curbă numită hidroizohipsă.

Totalitatea acestor curbe formează o hartă cu hidroizohipse (izopieze).

Din punct de vedere a condiţiilor de margine, atât acviferele cu suprafaţă liberă cât şi cele sub presiune se pot clasifica în următoarele tipuri:

Acvifere infinite, la care zona de alimentare se află la mare distanţă de puţul de extracţie; Acvifere semiinfinite, la care una din zonele de alimentare se află la mare distanţă iar în

apropiere de puţul de extracţie se află o limită impermeabilă sau o altă zonă de alimentare; Acvifere tip bandă, la care zonele de alimentare sau limitele impermeabile se află la mică

distanţă de puţul de extracţie şi sunt în general paralele; Acvifere unghiulare, la care limitele amintite mai sus formează un unghi oarecare,

cunoscut; Acvifere închise pe contur, la care limita de alimentare sau limita impermeabilă este

reprezentată de un contur închis, bine determinat, cunoscut.

Caracterul conservativ sau neconservativ al mişcării apelor subterane rezultă din condiţiile de alimentare a acviferului prin culcuş sau coperiş.

Acviferul conservativ nu işi modifică curgerea în funcţie de condiţiile pe frontierele superioară şi inferioară. De exemplu, un acvifer sub presiune care se desfăşoară între două straturi impermeabile se poate considera conservativ. Dacă el comunică, prin drenanţă, cu alte acvifere din zonă se va numi neconservativ. Un acvifer freatic este, de obicei, neconservativ pentru că este alimentat de la suprafaţa solului din precipitaţii. În lipsa precipitaţiilor el se poate considera conservativ.

3.2 Ecuaţia de difuzivitate pentru acviferele cu suprafaţă liberă

O pânză freatică liberă este un mediu poros care nu este saturat decât sub o anumită cotă. Deasupra acestei cote mediul este nesaturat. În acest caz se poate neglija compresibilitatea apei (ρ=ct) şi a mediului poros (n=ct). Toate variaţiile sarcinii vor antrena o mişcare a suprafeţei libere. 2

Vom considera o prismă dintr-un acvifer cu suprafaţă liberă, de înălţime h, aflată între substratul impermeabil b(x,y) şi suprafaţa liberă H(x,y,t). Vom presupune că în această pânză cu suprafaţă liberă toate vitezele sunt orizontale şi paralele între ele, pe aceiaşi verticală.

2 13de Marsily, G., “Hydrogeologie Quantitive”, Masson Editeur, Paris, 1981

Figura 3.4: Schematizarea unui volum elementar dintr-o pânză freatică

Presupunem că tensorul permeabilităţii admite verticala ca una din direcţiile principale. Vom considera sarcina H(x,y) ca necunoscută, independentă de z. Deci H reprezintă sarcina pe verticală şi în particular, cota suprafeţei libere a pânzei. Dacă nu există gradienţi de sarcină verticali, din legea lui Darcy rezultă că nu există componente verticale ale vitezei.

Vom alege axele x, y ca fiind cele două direcţii principale de anizotropie, în plan.

Ecuaţia de continuitate pentru prisma dx, dy, (H-b) se poate exprima făcând bilanţul intrărilor si ieşirilor, astfel:

a) Fluxul masic care intră în unitatea de timp prin cele două feţe perpendiculare pe Ox este:

Jx=Jx(x,y,z)–Jx (x+dx, y, z) (3.1)

J x =ρdy⋅ ∫b(x,y )

H (x,y,t )

U x( x,y,z )dz− ρdy ∫b(x,y )

H (x,y,t )

U x( x+dx,y,z )dz

(3.2)

Ux= viteza Darcy după direcţia x

Rezultă :

J x=− ρdy ∂∂ x [∫

b

H

U x dz ]dx

(3.3)

Conform legii lui Darcy, viteza în direcţia x este:

U x=−K xx ( x,y,z ) ∂ H∂ x

(3.4)

Iar fluxul masic după direcţia x devine:

J x =ρdxdy ∂∂ x [∫

b

H

K xx∂ H∂ x

dz ] (3.5)

În acelaşi mod se va calcula fluxul prin suprafeţele perpendiculare pe Oy.

J y=ρdxdy ∂∂ y [∫

b

H

K yy∂ H∂ y

dz ] (3.6)

Fluxul masic care intră în unitatea de timp prin cele patru feţe perpendiculare pe Ox respectiv pe Oy, este:

J=J x +J y =ρdxdy { ∂∂ x [∂ H

∂ x∫b

H

Kxx dz ]+ ∂∂ y [ ∂ H

∂ y∫b

H

K yy dz]} (3.7)

b) pentru a ţine seama de schimburile pânzei freatice cu exteriorul, se introduce în bilanţul masic:

Fie debitul masic prin unitatea de suprafaţă a pânzei freatice ( ρwdxdy )

, w fiind debitul prelevat pe unitatea de suprafaţă a pânzei freatice;

Fie debitul masic Qm de fluid prelevat în element dintr-o sursă exterioară (pozitiv dacă este injectat în acvifer şi negativ dacă este pompat, de exemplu printr-un puţ)

c) variaţia masei elementului considerat se va face prin ridicarea sau coborârea nivelului suprafeţei libere (apa este liberă să se ridice în stratul permeabil);

d) masa de apă gravitaţională conţinută în elementul de volum este:

dM=ρnd ( H−b ) dxdy (3.8)

nd= porozitatea de drenaj.

Variaţia acestei mase în timp va fi:

dMdt

=ρnd∂ H∂ t

dxdy (3.9)

Bilanţul intrărilor şi ieşirilor, ţinând cont de conservarea masei va fi reprezentat de ecuaţia:

ρdxdy { ∂∂ x [∂ H

∂ x∫b

H

K xxdz ]+ ∂∂ y [ ∂ H

∂ y∫b

H

K yy dz ]}=ρnd∂ H∂ t

dxdy± ρwdxdy

(3.10)

Cum ρ≠0

, dx , dy≠0

, rezultă:

∂∂ x [∂ H

∂ x∫

b (x,y)

H (x,y,t )

K xx dz ]+ ∂∂ y [∂ H

∂ y∫

b(x,y )

H (x,y,t )

K yy dz ]±w=nd∂ H∂ t

(3.11)

Această ecuaţie se numeşte ecuaţia de difuzivitate a pânzei freatice cu suprafaţă liberă. Ea este neliniară în H.

Cazuri particulare:

1. Dacă conductivităţile Kxx şi Kyy sunt constante pe toată verticala, ecuaţia difuzivităţii devine:

∂∂ x [∂ H

∂ xKxx h]+ ∂

∂ y [ ∂ H∂ y

K yy h]±w=nd∂ H∂ t

, în care H−b=h ( x,y,t )

(3.12)

2. Pentru un mediu anizotrop, ţinând seama de definirea transmisivităţii:

T xx=∫b

H

Kxx dz =K xx ( H−b )=K xx h

(3.13)

T yy=∫b

H

K yy dz=K yy( H−b )=K yy h

(3.14)

Ecuaţia difuzivităţii devine:

∂∂ x (T xx

∂ H∂ x )+ ∂

∂ y (T yy∂ H∂ y )±w=nd

∂ H∂ t

(3.15)

Vom presupune ca transmisivitatea variază puţin cu sarcina h, adică repartiţia verticală a lui K este astfel încât variaţiile lui h nu antrenează o variaţie a transmisivităţii, T, mai mare de 10%.

Dacă mediul este izotrop, transmisivitatea este constantă (T xx =T yy=T

), iar ecuaţia difuzivităţii devine:

(∂2 H∂ x2

+ ∂2 H∂ y2 )±w

T=

nd

T∂ H∂ t

(3.16)

3. Dacă stratul b(x,y) este orizontal şi luăm b(x,y)= 0, ca plan de referinţă pentru sarcina hidraulică ,

H ( x,y,t )−b( x,y )=h( x,y,t )=H ( x,y,t ) (3.17)

iar ecuaţia difuzivităţii devine:

∂∂ x (K xxh

∂h∂ x )+ ∂

∂ y (K yyh∂h∂ y )±w=nd

∂ h∂ t

(3.18)

4. Dacă mediul este izotrop şi uniform, Kxx=Kyy=K, ecuaţia difuzivităţii devine:

K ∂∂ x (h

∂ h∂ x )+K ∂

∂ y (h∂ h∂ y )±w=nd

∂h∂ t

(3.19)

Sau:

∂h∂ x

∂ h∂ x

+h∂2 h∂ x2

+ ∂ h∂ y

∂ h∂ y

+h∂2 h∂ y2

± wK

=nd

K∂h∂ t

(3.20)

Respectiv:

∇2h2±2wK

=2nd

K∂h∂ t

(3.21)

5. Dacă regimul este permanent,

∂h∂ t

=0

şi ecuaţia difuzivităţii devine:

∇2 h2=±2wK

(3.22)

6. Dacă pânza nu este alimentată la suprafaţă (w=0) şi regimul este permanent, ecuaţia difuzivităţii va fi:

∂2h∂ x2

+ ∂2h∂ y2

=0

(3.23)

3.3 Ecuaţia de difuzivitate pentru acviferele sub presiune

Se consideră un acvifer sub presiune (fig. 3.2) de grosime M(x,y), care se gaseste între două straturi impermeabile a(x,y) şi b(x,y). Mişcarea în acvifer se face datorită diferenţei de sarcină piezometrică între capetele acviferului (hA-hB). Suprafaţa piezometrică h se găseşte deasupra limitei superioare a acviferului. Dacă se face un foraj în acviferul sub presiune, nivelul apei din puţ se ridică deasupra nivelului coperişului şi reprezintă sarcina piezometrică din acvifer.2

2 13de Marsily, G., “Hydrogeologie Quantitive”, Masson Editeur, Paris, 1981

Figura 3.5 : Acviferul sub presiune

În cazul acviferului sub presiune atât apa cât şi matricea poroasă sunt supuse la eforturi de compresiune care pot modifica porozitatea mediului.

Ecuaţia de difuzivitate pentru acviferul sub presiune se obţine din expresia ecuaţiei de continuitate pentru un volum elementar reprezentativ, pentru care se iau în considerare atât compresibilitatea apei cât şi a scheletului solid. Compresibilitatea scheletului solid se ia în calcul prin modificarea porozităţii n a mediului poros, datorită acţiunii unui efort σ.

Raudkivi dă o demonstraţie pentru deducerea ecuaţiei de difuzivitate într-un caz particular, al unui mediu izotrop Kxx=Kyy=Kzz=K, considerând că tasările se fac în mod special dupa direcţia z.

Ecuaţia de continuitate pentru un volum paralelipipedic dV=dxdydz, având porozitatea n, saturat cu apă (cu densitatea ρ), este:

−[ ∂∂ x ( ρU x )+ ∂

∂ y ( ρU y )+ ∂∂ z (ρU z) ]dxdydz= ∂

∂ t(nρΔxΔyΔz )

(3.24)

în care Ux, Uy, Uz sunt componentele vitezei Darcy.

Dacă mediul poros este compresibil, porozitatea mediului poate varia în spaţiu şi timp.

Introducând noţiunile de compresibilitate a solului şi apei (coeficientul de compresibilitate al solului, α ; coeficientul de compresibilitate al apei, β≈4,829 * 10-10 m2N-1; modulul de elasticitate al solului, Es; modulul de elasticitate al apei Ew=2,07 GPa), ecuaţia de continuitate devine:

[ ∂∂ x ( ρU x )+ ∂

∂ y ( ρU y )+ ∂∂ z (ρU z )] ΔxΔyΔz=−∂ p

∂ tΔxΔyΔzρβ [n+

αβ ]

(3.25)

Se introduc în ecuaţia de continuitate:

Presiunea p, exprimată în funcţie de sarcina piezometrică h:

p=ρgh−ρgz

(3.26)

Vitezele Darcy după cele trei direcţii de curgere:

U x=−K∂ h∂ x

; U y=−K

∂h∂ y

; U z=−K

∂ h∂ z

(3.27)

Se obţine o ecuaţie de forma:

∂2h

∂ x2+∂2h

∂ y2+∂2 h

∂ z2=

ρβ (n+αβ ) g

K∂ h∂ t

(3.28)

Mărimea:

Ss =ρgβ(n+αβ )

; [Ss]SI=[L-1] (3.29)

se numeşte coeficient specific de înmagazinare.

Ecuaţia de continuitate pentru cazul unui acvifer sub presiune având conductivitate hidraulică constantă va fi:

∂2 h∂ x2

+ ∂2 h∂ y2

+ ∂2 h∂ z2

=Ss

K∂ h∂ t

(3.30)

Pentru cazul general în care tensorul conductivităţilor hidraulice nu are toate componentele constante, ecuaţia de continuitate pentru acviferele sub presiune se scrie sub forma:

div ( ¯̄K⋅gradh )=Ss∂h∂ t

+q (3.31)

în care:

h este sarcina piezometrică într-un punct din acvifer,

¯̄K este tensorul conductivităţii hidraulice,

q reprezintă debitul injectat sau pompat raportat la grosimea acviferului,

Ss este coeficientul specific de înmagazinare.

Se definesc transmisivităţile acfiverului după direcţiile principale de anizotropie:

T xx=∫a

b

Kxx dz

- transmisivitatea după direcţia x (3.32)

T yy=∫a

b

K yy dz

- transmisivitatea după direcţia y (3.33)

şi coeficientul de înmagazinare al acviferului sub presiune:

S=∫a

b

Ss dz

(3.34)

unde a este cota culcuşului iar b cota coperişului.

Presupunând Kxx, Kyy, Ss constante pe toată înălţimea M a acviferului, se pot calcula transmisivităţile:

T xx =K xx M (3.35)

T yy =K yy M (3.36)

şi coeficientul de înmagazinare al acviferului sub presiune S

S=Ss M=ρgMβ (n+αβ )

(3.37)

având valori cuprinse între 10-3 şi 10-5.

Dacă culcuşul şi coperişul sunt paralele (M=const), ecuaţia de difuzivitate devine:

div ( ¯̄T⋅gradh)=S∂ h∂ t

+Q (3.38)

în care:

- ¯̄T

este tensorul transmisivităţii;

- Q este aportul de debit din exterior (+Q daca este injectat sau –Q daca este pompat).

Dacă mediul poros este izotrop, transmisivitatea T este constantă după toate direcţiile, iar ecuaţia de difuzivitate se poate scrie sub forma:

∂2h∂ x2

+ ∂2h∂ y2

+∂2 h∂ z2

= ST

∂h∂ t

+ QT

(3.39)

Raportul S/T este numit difuzivitatea acviferului.

Dacă mişcarea este permanentă (∂h∂ t )=0

şi Q=0 ecuaţia de difuzivitate:

∂2h∂ x2

+ ∂2h∂ y2

+∂2 h∂ z2

=0

(3.40)

descrie o mişcare potenţială.

Rezolvarea unei probleme de curgere într-un acvifer sub presiune constă în integrarea ecuaţiei de difuzivitate exprimată sub una din formele (3.38 – 3.40), în condiţii de unicitate specifice cazului analizat.

Condiţiile de unicitate sunt reprezentate prin condiţiile iniţiale şi prin condiţiile pe frontiera domeniului.

Condiţiile iniţiale constau în specificarea valorilor sarcinii piezometrice, respectiv ale poziţiei suprafeţei piezometrice, la momentul iniţial:

h (x,y,z,t o)=ho (x,y,z ) (3.41)

Condiţiile pe frontieră pentru acviferul din figura 3.2 sunt:

Pentru x=0, h (0, y,z,t ) =hA ( y,z,t )

; (3.42)

Pentru x=L, h ( L,y,z,t )=hB ( y,z,t )

(3.43)

Condiţiile pe frontieră pot fi şi de alte tipuri. De exemplu, se poate specifica fluxul printr-o frontieră. Pe o frontieră impermeabilă se poate pune condiţia ca debitul prin acea suprafaţă să fie nul.

Integrarea ecuaţiei de difuzivitate se poate face prin metode numerice (metoda diferenţelor finite sau metoda elementului finit) ori prin metode analitice, în cazuri simplificate.

De exemplu, dacă mişcarea este permanentă, iar acviferul este conservativ şi izotrop, se poate considera conductivitatea hidraulică constantă, ecuaţia de difuzivitate corespunzătoare este (3.40) şi curgerea se poate analiza folosind metoda mişcărilor potenţiale plane.

Soluţia ecuaţiei de difuzivitate este reprezentată de funcţia h(x,y,z,t) în toate punctele domeniului, la diferite momente de timp. Ea dă poziţia suprafeţei piezometrice a acviferului faţă de un nivel de referinţă. Având această funcţie se pot face diferite calcule legate de analiza curgerii în acvifer:

Se pot trasa linii de egal nivel al suprafeţei piezometrice faţă de nivelul de referinţă (hidroizohipse);

Se pot trasa liniile de curent;

Se poate calcula viteza în orice punct din domeniu, timpul de parcurgere al acviferului în lungul unei linii de curent şi debitul printr-o suprafaţă normală la direcţia de curgere, debitul intrat în acvifer prin limitele de alimentare.

3.4 Metode numerice de rezolvare a ecuaţiei de difuzivitate

Ecuaţia de difuzivitate atât pentru acviferul cu suprafaţă liberă cât şi pentru acviferul sub presiune este o ecuaţie cu derivate parţiale de ordinul 2, neliniară în general în cazul acviferelor anizotrope.

Integrarea în condiţii de unicitate a unei astfel de ecuaţii ( condiţii iniţiale şi condiţii pe frontieră) se poate face în general numeric şi în anumite cazuri particulare analitic.

În general se pot găsi anumite ipoteze hidraulice care pot simplifica problema tridimensională transformând-o într-o problemă bidimensională. De exemplu pentru un mediu izotrop în cazul mişcării permanente se poate folosi metoda mişcărilor potenţiale plane.

Astfel, dacă ecuaţia de difuzivitate se poate exprima sub forma Δh=0

,

∂2h∂ x2

+ ∂2h∂ y2

=0

, folosirea principiului suprapunerilor mişcărilor potenţiale permite rezolvarea unor probleme de curgere în acvifere.

Alte ipoteze hidraulice ce pot fi folosite sunt:

1. mişcarea în acvifer se face astfel încât liniile de curent sunt aproximativ paralele (ipoteza Depuit) dacă panta hidraulică a suprafeţei piezometrice respectiv a suprafeţei libere este foarte mică.

2. mişcarea spre un puţ de pompare este aproximativ radială.

În astfel de situaţii există posibilitatea abordării analitice a problemei de infiltraţii. În majoritatea situaţiilor reale mediul poros este anizotrop existând de cele mai multe ori direcţii preferenţiale de anizotropie.

În astfel de situaţii metodele analitice nu pot fi folosite şi trebuie folosite metode numerice, respectiv metoda diferenţelor finite sau a elementului finit.

3.4.1 Metoda diferenţelor finite

Cea mai veche şi folosită tehnică numerică este metoda diferenţelor finite care constă în înlocuirea derivatelor parţiale din ecuaţia de difuzivitate prin diferenţe finite:

∂h∂ x

= limΔx→0

ΔhΔx

≈ ΔhΔx

(3.44)

O astfel de aproximaţie este cu atât mai corectă cu cât pasul de spaţiu Δx este mai mic.

În locul unei funcţii continue h care satisface ecuaţia de curgere într-un număr infinit de puncte, soluţiile numerice asigură valori într-un număr finit de puncte numite noduri.

Ecuaţia de curgere este considerată în fiecare nod şi toate derivatele parţiale sunt înlocuite prin diferenţe finite utilizând valorile sarcinii din noduri. Astfel se obţine o ecuaţie algebrică pentru fiecare nod. Toate aceste ecuaţii în diferenţe finite formează un sistem de ecuaţii care poate fi rezolvat, rezultând valoarea sarcinii hidraulice în fiecare nod.

Dacă problema este nepermanentă se face o discretizare a derivatei în timp alegându-se un număr limitat de intervale de timp, sarina în noduri fiind calculată la aceste valori de timp. Soluţia obţinută este o soluţie aproximativă.

Eroarea rezultată între soluţia exactă şi soluţia numerică este numită eroare de trunchiere. Dacă această eroare este suficient de mică, soluţia obţinută prin metoda diferenţelor finite este considerată convergentă. Obţinerea unei astfel de situaţii presupune alegerea unui număr mare de noduri în domeniul de curgere.

Cu cât sistemul de ecuaţii în diferenţe finite este mai mare cu atât integrarea acestuia este mai dificilă, apărând erori. Astfel de erori apar datorită tehnicii de aproximare a soluţiei sistemului de ecuaţii sau din cauza capacităţii finite a computerului. Diferenţa dintre soluţia exactă şi soluţia aproximativă a setului de ecuaţii este numită eroare de rotunjire. Când această eroare este suficient de mică soluţia este considerată stabilă. Cu cât numărul de noduri este mai mare cu atât obţinerea stabilităţii este mai greu de realizat. Astfel, condiţiile de obţinere a unei soluţii convergente intră în conflict cu condiţiile de obţinere a unei soluţii stabile.

În practică va exista întotdeauna o conexiune între erorile de trunchiere şi erorile de rotunjire care va depinde de tipul de model şi de performanţele sistemului de calcul.

Soluţiile practice presupun folosirea unor date din teren pentru calibrarea modelelor.

Acurateţea aproximaţiilor în diferenţe finite este foarte importantă în model. Există diferite posibilităţi de aproximare. De exemplu considerăm dezvoltarea în serie Taylor a funcţiei sarcinii hidraulice h(x):

h ( x+Δx )=h ( x )+ Δxdhdx

+( Δx )2

2 !d2 h

dx2+

( Δx )3

3!d3h

dx 3+…

(3.45)

unde:

h(x) este valoarea sarcinii hidraulice într-un nod (x);

h(x+Δx) este sarcina hidraulică într-un nod aflat la o distanţă Δx de celalalt punct.

O aproximaţie în diferenţe finite rezultă pentru prima derivată:

dhdx

=h ( x+ Δx )−h ( x )

Δx+ Δx (… )

(3.46)

Când derivata este aproximată prin primul termen din partea dreaptă a ecuaţiei (3.46) se obţine o aproximaţie în diferenţe finite forward:

dhdx

≈h ( x+ Δx )−h (x )

Δx (3.47)

Eroarea făcută este dată de termenul al doilea din ecuaţia (3.46), este proporţională cu Δx şi este cu atât mai mică cu cât Δx este mai mic.

Dacă considerăm o dezvoltare în serie Taylor luând un pas înapoi (x-Δx):

h ( x−Δx )=h ( x )−Δxdhdx

+( Δx )2

2 !d2 h

dx2−

( Δx )3

3!d3 h

dx3+…

(3.48)

rezultă pentru prima derivată o aproximaţie în diferenţe finite de tip backward:

dhdx

=h ( x )−h ( x−Δx )

Δx+Δx (… )

(3.49)

dhdx

≈h ( x )−h ( x−Δx )

Δx (3.50)

Înlocuind h(x) din relaţia (3.45) în ecuaţia (3.48) se obţine:

h ( x+Δx )−h ( x−Δx )=2 Δxdhdx

+2 ( Δx )3

3 !d3 h

dx3+…

(3.51)

Rezultă:

dhdx

=h ( x+ Δx )−h ( x−Δx )

2 Δx+ ( Δx )2 (… )

(3.52)

Prin trunchiere rezultă o aproximaţie în diferenţe finite centrată, având o eroare de trunchiere proporţională cu (Δx)2 care este mult mai mică decât cea obţinută prin aproximaţii forward şi backward:

dhdx

≈h ( x+ Δx )−h (x−Δx )

2 Δx (3.53)

Dacă ecuaţiile (3.45) şi (3.48) sunt adunate se obţine relaţia:

h ( x+Δx )−h ( x−Δx )=2 h ( x )+ 2 ( Δx )2

2 !d2 h

dx 2+( Δx )4 (… )

(3.54)

din care prin trunchiere se poate obţine aproximarea în diferenţe finite pentru derivata de ordinul 2:

d2 hdx2

≈h ( x+Δx )−2h ( x )+h ( x−Δx )

( Δx )2 (3.55)

Eroarea de trunchiere este proporţională cu (Δx)4.

Aproximaţiile în diferenţe finite centrate prezintă probleme la frontiera domeniului deoarece în astfel de puncte nu pot fi luaţi paşi forward şi backward. În aceste puncte trebuie impuse condiţii pe frontieră specifice fiecărei probleme.

Prin considerarea mai multor termeni din dezvoltările în serie Taylor se obţin soluţii aproximative mai apropiate de cea reală dar sistemul de ecuaţii rezultat va fi mai complicat şi erorile de rotunjire vor fi mai mari.

Cele mai utilizate în practică sunt aproximaţiile prin diferenţe finite centrate.

3.4.2 Soluţii ale ecuaţiei de difuzivitate pentru cazul mişcării permanente uniforme (mediu poros omogen şi izotrop).

Pentru cazul în care ecuaţia de difuzivitate poate fi reprezentată prin ecuaţia lui Laplace în două dimensiuni:

∂2h∂ x2

+ ∂2h∂ y2

=0

(3.56)

domeniul de curgere este împărţit printr-o reţea de puncte cu pasul Δx= Δy ca în figura 2.6.

Figura 3.6 : Prezentarea reţelei de discretizare

Nodurile sunt considerate situate la intersecţia liniilor reţelei. Un nod din această reţea poate fi notat prin numărul liniei în direcţia x, respectiv y sub forma (i,j) iar sarcina hidraulică într-un astfel de nod se notează h(i,j), respectiv:

hi , j=h (x i , y j ) (3.57)

Frontiera domeniului de curgere este aproximată prin segmente drepte în lungul liniilor reţelei aşa cum este arătat în figura 3.7, astfel încât nodurile interioare pot fi diferenţiate de nodurile frontierei (nodurile exterioare nu sunt luate în considerare).

Figura 3.7 : Reprezentarea frontierei reale şi frontierei aproximative

Un nod interior tipic are patru noduri vecine (figura 3.8)

Figura 3.8 : Reprezentarea unui nod din reţea

Derivata de ordinul 2 în raport cu x din ecuaţia (3.56) poate fi aproximată printr-o schemă în diferenţe finite centrate ca în ecuaţia (3.55), rezultând:

∂2h∂ x2

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

( Δx )2 (3.58)

Pentru derivata în raport cu y rezultă:

∂2 h∂ y2

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

( Δx )2 (3.59)

Înlocuid aceste derivate în relaţia (3.56) rezultă o ecuaţie algebrică cu diferenţe finite:

hi−1 , j−2 hi , j+hi+1, j

( Δx )2+

h i , j−1−2 hi , j+hi , j+1

( Δx )2=0

(3.60)

hi−1, j+hi , j−1−4 hi , j+hi+1 , j+h i , j+1=0 (3.61)

Această ecuaţie se aplică în toate nodurile interioare dar nu poate fi scrisă în nodurile aflate pe frontieră. Din această cauză condiţiile pe frontieră trebuie să fie transformate în ecuaţii algebrice.

Fie nodul (i,j) aşezat pe o frontieră. În acest punct se poate pune o condiţie de sarcină impusă:

hi , j=h0 ( xi , y j ) (3.62)

Dacă se pune condiţie de flux impus pe frontieră, datorită aproximării frontierei curbilinii prin segmente drepte în lungul liniilor reţelei, fluxul prin frontiera reală trebuie să fie transformat în flux prin frontiere atât în direcţia x cât şi în direcţia y .

Figura 3.9 : Reprezentarea descompunerii fluxului pe frontiera domeniului

Fluxul unitar prin frontieră (viteză Darcy) se poate scrie:

qxo=−K

∂h∂ x

(3.63)

q yo=−K

∂ h∂ y

(3.64)

Presupunem cunoscut debitul în direcţia x: qx

0 ( xi , y j ). Acesta poate fi exprimat în funcţie

de valorile sarcinii în punctele (i+1,j) şi (i,j) cu relaţia (3.65) (diferenţă forward) sau cu relaţia (3.67) (diferenţă centrată).

qxo ( xi , y j )=−K

hi+1 , j−hi , j

Δx (3.65)

qxo ( xi , y j )=−K

hi+1 , j−hi−1 , j

2 Δx (3.66)

Aproximaţia făcută cu relaţia (3.66) este mai bună dar nu se cunoaşte valoarea sarcinii hidraulice din punctul (i-1,j). Pentru a rezolva această problemă se consideră un nod imaginar (i-1,j) şi din relaţia (3.66), dacă se cunoaşte valoarea debitului intrat prin frontieră în punctul (i,j), se determină valoarea sarcinii hidraulice în punctul (i-1,j) cu relaţia:

hi−1, j=h i+1 , j+2 qxo (x i , y j ) Δx /K

(3.67)

Dacă se introduce valoarea sarcinii din nodul imaginar (i-1,j) în ecuaţia (5.18) se obţine o ecuaţie similară care ţine seama de debitul prin frontieră în punctul (i,j):

hi , j−1−4hi , j+2hi+1, j+h i , j+1=−2qxo (x i , y j ) Δx / K

(3.68)

Ecuaţia (3.68) se poate scrie în toate punctele (i,j) de pe frontiera aproximativă pentru ca nu sunt necesare valori ale sarcinii în punctele (i-1,j). La fel se poate proceda după direcţia y.

În cazul în care debitul prin frontieră este zero, din ecuaţia (3.67) rezultă că sarcina în punctul imaginar (i-1,j) este egală cu sarcina din punctul (i+1,j).

Aprecierea condiţiilor pe frontieră este o problemă relativ complicată care necesită cunoaşterea condiţiilor fizice reale din domeniu.

Ecuatia (3.61) se poate scrie în toate punctele (i,j) din domeniu mai putin în cele aflate în imediata apropiere a frontierelor. În aceste puncte (din apropierea frontierelor) se vor scrie ecuaţii de forma (3.62) dacă condiţiile sunt de tipul sarcină impusă sau de tipul (3.68) dacă condiţiile sunt de tipul debit impus pe frontieră.

Astfel se vor obtine M sisteme de N ecuaţii cu N necunoscute, unde M este numărul de noduri după direcţia y (j=1:M) iar N este numărul de noduri după direcţia x (i=1:N).

O metoda simplă de rezolvare a unei astfel de probleme este metoda relaxării. Aceasta este o metodă iterativă în care valoarea dintr-un nod (i,j) se obţine ca medie aritmetică a valorii sarcinii din nodurile vecine:

hi , j=hi−1 , j+h i , j−1+hi+1, j+h i , j+1

4 (3.69)

Această ecuaţie se poate scrie în toate nodurile din reţea mai puţin în cele aflate pe frontiere. Pentru nodurile aflate pe frontieră se foloseşte relaţia:

hi , j=2q x

0 ( xi , y j ) Δx+hi , j−1+2hi+1, j+h i , j+1

4 K (3.70)

Se porneşte de la valori aproximative ale sarcinii în tot domeniul şi se obţine pentru

iteraţia m+1 valoarea sarcinii hi , j

m+1

în funcţie de valorile sarcinii în iteraţia m:

hi , jm+1=

hi−1 , jm +hi , j−1

m +hi+1 , jm +hi , j+1

m

4 (3.71)

Oprirea calcului iterativ se face punând condiţia ca diferenţa dintre valoarea sarcinii în iteraţia m şi cea din iteraţia m+1 să fie mai mică decât o valoare impusă, numită toleranţă.

Această metodă este numită metoda iterativă Jacobi:

max|hi , jm+1−hi , j

m |<TOL (3.72)

O metodă mai rapidă, numită metoda iterativă Gauss-Seidel, presupune determinarea sarcinii hidraulice în nodul (i,j) în iteraţia m+1, ţinându-se seama de valorile sarcinii hidraulice în nodurile (i+1,j) şi (i,j+1) în iteraţia anterioară m şi de valorile sarcinii hidraulice din nodurile (i-1,j) şi (i,j-1) la iteraţia m+1, obţinându-se o schemă semiimplicită:

hi , jm+1=

hi−1 , jm+1 +hi , j−1

m+1 +hi+1 , jm +hi , j+1

m

4 (3.73)

O altă metodă, numită suprarelaxare (figura 3.10), presupune calcularea valorii sarcinii în iteraţia m+1 ţinându-se seama de un factor de suprarelaxare ω:

hi , jm+1= (1−ω )hi , j

m +ω (h i−1, jm+1 +hi , j−1

m+1 +hi+1, jm +hi , j+1

m ) /4 (3.74)

Factorul ω este cuprins între 0 şi 2.

Figura 3.10 : Metoda suprarelaxării

Pentru a vedea cât de bună este metoda folosită se poate proceda astfel: se variază valoarea lui Δx (de exemplu se dublează) şi se calculează diferenţa maximă dintre rezultatele obţinute. O altă metodă este de a compara soluţia numerică cu soluţii analitice existente în cazuri mai simple.

Pentru cazul în care mediul nu este uniform, metoda prezentată anterior este mai complicată. Presupunem ca ecuaţia de difuzivitate pentru un acvifer este:

∂∂ x [T ( x , y ) ∂ H

∂ x ]+ ∂∂ y [T ( x , y ) ∂ H

∂ y ]+R ( x , y )=0

(3.75)

Transmisivitatea şi termenul sursă variază în funcţie de poziţie. Se presupun cunoscute poziţia şi debitul unor puţuri de pompare. Domeniul de curgere este acoperit cu o reţea de noduri (i,j) pentru care nodurile interne sunt reprezentate în figura 3.11.

Figura 3.11 : Reprezentarea nodurilor interne

O metoda de rezolvare a problemei poate fi ca în loc de a obţine o aproximaţie în diferenţe finite pentru ecuaţia de curgere (3.74) să se facă un bilanţ al debitului intrat şi ieşit din domeniul aflat în jurul nodului (i,j).

Scriind că debitele de apă intrate în volumul elementar paralelipipedic de dimensiuni Δx: Δx:1 sunt egale cu cele ieşite, rezultă:

Q1+Q2+Q3+Q4+R i , j ( Δx )2−W i , j=0 (3.76)

în care:

Q1, Q2, Q3, Q4 sunt debitele pe cele patru feţe ale paralelipipedului elementar;

Ri,j reprezintă încărcarea medie pe suprafaţa (Δx)2;

Wi,j debitul pompat din toate puţurile situate în pătratul considerat.

Ecuaţia de bilanţ va fi:

Q1+Q2+Q3+Q4+R i , j ( Δx )2−W i , j=0 (3.77)

Debitul Q1 se calculează folosind relaţia lui Darcy:

Q1=U(i−1

2, j)⋅A

( i−12

, j)=U

( i−12

, j)⋅Δx⋅1=−K

ΔHΔx

⋅Δx⋅1=−T( i−1

2, j)

H i , j−H i−1, j

ΔxΔx=

¿−T(i−1

2, j) (

H i , j−H i−1 , j ) (3.78)

La fel sunt gândite debitele Q2, Q3, Q4.

Ti-1/2,j reprezintă transmisivitatea în punctul (i-1/2, j) şi se poate obţine ca o medie aritmetică a valorilor din nodurile alăturate dar este recomandabil să se folosească o medie armonică:

Ti−1

2, j=

2 T i , j⋅T i−1, j

T i , j+T i−1, j (3.79)

Înlocuindu-se valorile debitelor în ecuaţia (3.75) se poate scoate valoarea sarcinii hidraulice într-un punct (i,j) cu o soluţie de forma:

H i , j=

T( i−1

2, j)

H i−1 , j+T( i , j−1

2 )H i , j−1+T

( i+ 12

, j)H i+1, j+T

( i , j+ 12 )

H i , j+1+Ri , j ( Δx )2−W i , j

T( i−1

2, j)

+T( i , j−1

2 )+T

(i+ 12

, j)+T

(i , j+12 )

(3.80)

Pentru stabilirea condiţiilor pe frontieră se procedează analog ca în cazul precedent. Condiţiile pe frontieră pot fi: sarcină impusă sau debit impus. Pentru cazul cu debit impus (fig.3.12), dacă în punctul (i,j) aflat pe frontieră se impune un debit cunoscut Q1, sarcina în punctul (i,j) se va calcula cu relaţia:

H i , j=

Q1+T(i , j−1

2 )H i , j−1 /2+T

(i+ 12

, j)H i+1 , j+T

(i , j+ 12)

H i , j+1/2+Ri , j ( Δx )2/2−W i , j

T(i , j−1

2 )/2+T

( i+ 12

, j)+T

(i , j+ 12)/2

(3.81)

Figura 3.12 : Condiţiile pe frontieră pentru un debit impus Q1.

3.4.3 Soluţii pentru curgerea nepermanentă

Fie un caz simplu de curgere nepermanentă, unidimensională, în lungul axei x descrisă de ecuaţia:

S∂ H∂ t

=T∂2 H∂ x2

(3.82)

Se consideră un pas de spaţiu Δx şi un pas de timp Δt ca în figura 2.13

Figura 3.13 : Reprezentarea paşilor de timp şi spaţiu

Sarcina hidraulică într-un nod i la un moment dat t se scrie:

H it=H (x i , t )

(3.83)

O aproximaţie în diferenţe finite centrată pentru derivata în timp în nodul i este dată de:

∂ H∂ t

=H i

t+ Δt−H it−Δt

2 Δt (3.84)

Considerând pentru ecuaţia (3.82) o aproximare prin diferenţe difinte centrată în timp (2.84) şi în spaţiu (3.58), rezultă ecuaţia:

SH i

t+Δt−H it−Δt

2 Δt=T

H i+1t −2H i

t+H i−1t

( Δx )2 (3.85)

Întrucât sunt necesari trei paşi de timp t-Δt, t, t+Δt se preferă folosirea unei scheme în diferenţe finite de tip forward:

∂ H∂ t

=H i

t+ Δt−H it

Δt (3.86)

care conduce la ecuaţia (3.86) din care se poate calcula sarcina hidraulică în toate nodurile reţelei la momentul t+Δt în funcţie de valorile la momentul anterior t, de pasul de timp şi de caracteristicile acviferului.

SH i

t+Δt−H it

Δt=T

H i+1t −2H i

t+H i−1t

( Δx )2 (3.87)

H it+ Δt=[ TΔt

S ( Δx )2 ]H i−1t +[1−

2TΔt

S ( Δx )2 ]H it+[ TΔt

S ( Δx )2 ]H i+1t

(3.88)

Folosirea acestei ecuaţii este foarte simplă. Pornind de la momentul de timp t=0 se pot calcula în orice punct valorile sarcinii hidraulice la momentul Δt. Această metodă se numeşte metodă în diferenţe finite explicită. Pentru obţinerea de rezultate bune pasul de timp trebuie să fie foarte mic (pentru micşorarea erorii de trunchiere). O condiţie suficientă pentru convergenţă este ca coeficienţii din partea dreapta a ecuaţiei (3.88) să fie cuprinşi între 0 şi 1, ceea ce presupune un pas de timp:

Δt≤S ( Δx )2

2T (3.89)

Rezultă că timpul de calcul este foarte mare pentru cazuri practice reale.

O altă aproximaţie poate fi de tip implicit:

SH i

t+Δt−H it

Δt=T

H i+1t+ Δt−2 H i

t+ Δt+H i−1t+ Δt

( Δx )2 (3.90)

rezultând ecuaţia:

[ −TΔt

S ( Δx )2 ]H i−1t +Δt+[1+ 2TΔt

S ( Δx )2 ]H it +Δt+[ −TΔt

S ( Δx )2 ]H i+1t+Δt=H i

t

(3.91)

în care necunoscutele sunt valorile sarcinii la timpul t+Δt în toate nodurile reţelei. Astfel, scriindu-se ecuaţia (3.91) în toate nodurile reţelei se obţine un sistem de ecuaţii, în general tridiagonal. Pasul de timp poate fi luat în general de două ori mai mare decât în cazul schemei explicite.

Se pot alege scheme de tipul:

SH i

t+Δt−H it

Δt=T

H i+1t+ Δt/2−H i

t+ Δt /2+H i−1t+ Δt /2

( Δx )2 (3.92)

în care sarcina hidraulică la momentul t+ Δt/2 este dată de relaţia:

H it+ Δt /2=( H i

t+H it+ Δt ) /2

(3.93)

rezultând o ecuaţie asemănătoare cu (3.91):

[−TΔt

2S ( Δx )2 ]H i−1t+ Δt+[1+

TΔt

S ( Δx )2 ]H it+ Δt+[−TΔt

2 S ( Δx )2 ]H i+1t+ Δt=[TΔt

2S ( Δx )2 ]H i−1t +

+[1−TΔt

S ( Δx )2 ]H it+[TΔt

2S ( Δx )2 ]H i+1t

(3.94)

Această schemă este numită metoda în diferenţe finite centrată sau metoda Crank-Nicholson, fiind asemănătoare cu metoda implicită.

Ecuaţiile (3.91) sau (3.94) pot fi scrise în toate nodurile reţelei, mai puţin în nodurile aflate pe frontiere i=1 şi i=N. În aceste puncte se impun condiţiile pe frontieră de tip sarcină impusă sau de tip debit impus. De exemplu dacă pe frontiera din stânga se cunoaşte o valoare de debit impus Q1, ecuaţia (5.50) devine:

[1+ 2TΔt

S ( Δx )2 ]H1t +Δt+[−2TΔt

S ( Δx )2 ]H2t +Δt=H1

t +2Q1 Δt

SΔx (3.95)

Unde Q1 este debitul de intrare în domeniu pe întreaga grosime a acviferului.

Toate ecuaţiile în diferenţe finite combinate pot fi scrise în formă matricială:

A⋅H=b (3.96)

unde:

H este vectorul necunoscutelor (sarcina în noduri la timpul t+Δt);

b este un vector al coeficienţilor din partea dreaptă a ecuaţiei (3.94)

Soluţia poate fi obţinută prin metoda relaxării care necesită însă un timp de calcul mare. Având în vedere forma trigiagonală a matricii A, sistemul (3.96) se poate integra uşor (algoritmul Thomas):

a i H i−1+bi H i+c i H i+1=d i (3.97)

e1=−c1

b1 :

f 1=d1

b1

for i=2 to n : ei= - ci/(bi+aiei-1) : fi=(di-aifi-1/)/(bi+aiei-1)

Hn=fn

for i=n to 2 : Hi+1=eiHi+fi (3.98)