Post on 15-Jun-2015
transcript
Dispersia poluanşilor în apele subterane 49
Capitolul 3 MODELAREA MATEMATICĂ A CURGERII APELOR SUBTERANE 3.1 NOŢIUNI GENERALE PRIVIND MODELAREA MATEMATICĂ A CURGERII APELOR SUBTERANE
Un fluid newtonian este un fluid izotrop, la care presiunea nu depinde decât de variabilele de stare ρ şi T0, iar tensorul de vâscozitate are o formă liniară în funcţie de gradientul vitezei.
În mecanica şi în termodinamica fluidelor newtoniene, toate problemele de curgere se reduc la determinarea a şase necunoscute:
ρ - densitatea fluidului [ ML-3] ; p - presiunea [ ML-1T-2] ; T0 - temperatura [ T ] ; vx,vy,vz - componentele câmpului de viteze [ LT-1 ],
având la dispoziţie : • -cele trei ecuaţii ale sistemului Navier-Stokes , care exprimă, pentru un fluid
vâscos, principiul fundamental al dinamicii: ρ ∑ ⋅= amf ρ (3.1) • ecuaţia de continuitate, care exprimă conservarea masei pentru un volum fix
:
( )t
vdiv∂ρ∂
+⋅ρρ =0, (3.2)
• ecuaţia de stare a fluidului
ρ = ρ0 eβ (p-po) , (3.3)
• ecuaţia transportului conductiv şi convectiv al căldurii prin fluid.
50 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
În mediul poros se poate considera, în general, curgerea izotermă (dispare necunoscuta T0) Având în vedere specificul aspectelor fizice, ale curgerii în medii poroase, aceste ecuaţii pot fi înlocuite cu altele ce sunt obţinute pe baza unor cercetări experimentale. De exemplu sistemul de ecuaţii Navier-Stokes se înlocuieşte cu o ecuaţie ce descrie legea lui Darcy şi care permite o simplificare considerabilă a abordării matematice a problemei. Determinarea celor şase necunoscute, amintite mai sus, se face prin realizarea unui model al realităţii fizice. Modelele matematice utilizate în fizica solului sau în hidrogeologie se pot clasifica în două categorii : modele deterministe şi modele stohastice. În modelele deterministe, parametrii şi variabilele au o valoare perfect determinată şi rezultatul este unic. Modelele empirice stabilesc o relaţie între o caracteristica necunoscută a solului şi alte proprieţăţi ale acestuia fără să ia în considerare mecanismele fundamentale. Cele mai comune sunt modelele regresive care reprezintă corelaţii simple sau multiple între un parametru necunoscut şi celelalte caracteristici ale solului ( funcţii de pedotransfer ) . Modelele conceptuale se bazează pe concepte, adică pe o schemă de funcţionare incompletă voit, care simplifică realitatea. Modelele funcţionale se bazează pe o schematizare grosieră a realitaţii. Ele sunt simple din punct de vedere matematic, necesită un număr mic de date de intrare, sunt uşor de rezolvat şi sunt folosite, în special , pentru gestiunea resurselor . Modelele mecaniciste descriu procesele la scară macroscopică prin ecuaţii cu derivate parţiale. Aceste ecuaţii sunt deduse din legile fizice ce guvernează procesele de transfer (Darcy, Fick, Fourier, legea de continuitate ). Astfel de modele introduc un mare număr de parametrii, se rezolvă, în general , prin metode numerice şi trebuie să fie verificate prin încercări experimentale . Numărul mare de parametrii necesari limitează uneori folosirea modelelor în condiţiile de teren . În modelele stohastice variabilele de intrare şi parametrii sunt mărimi aleatoare, reprezentate prin funcţii de distribuţie de probabilitaţi. Rezultatele sunt, de asemenea, caracterizate de o funcţie de distribuţie. Modelele stohastice non-mecaniciste fac apel la o funcţie de transfer care transformă semnalul de intrare într-un semnal de ieşire ţinând seama, într-un mod global, de totalitatea proceselor care se desfăşoară în sistem. Modelele stohastice mecaniciste iau în consideraţie variabilitatea spaţială a datelor de intrare, luându-se drept funcţii de distribuţie de probabilitaţi. Aceste date sunt introduse în modelul mecanicist. Intoducând un mare număr de astfel de date se obţine o lege de distribuţie a variabilelor de ieşire. Modelele matematice se mai pot clasifica după modul de rezolvare: - modele analitice ;
Dispersia poluanşilor în apele subterane 51 - modele numerice ; şi după obiectivul lor : cercetare, gestiune, regularizare, educaţie, prognoză. 3.2. IPOTEZE SIMPLIFICATOARE ÎN MODELAREA MATEMATICĂ A CURGERII IN MEDII POROASE
În modelele mecaniciste se fac, în general, următoarele ipoteze simplificatoare :
• matricea poroasă este rigidă ( de multe ori este considerată un mediu omogen şi izotrop ),
• faza lichidă este incompresibilă, • faza gazoasă este continuă şi la presiunea atmosferică, • curgerea se face la temperatură constantă, • diferite mărimi care intervin în transfer ( flux, conţinut de apă, viteză … )
sunt reprezentate prin valori medii la scară macroscopică . 3.3. MODELAREA MATEMATICA A TRANSFERURILOR CE AU LOC ÎNTR-UN MEDIU POROS Transferurile de materie sau de energie într-un sol, indiferent de natura acestora ( apă, gaz, soluţii, căldură ), constau în suprapunerea a două procese :
• o mişcare descrisă printr-o lege dinamică ( mişcarea poziţiei particulelor în raport cu matricea solidă ),
• o variaţie a stocurilor în timp (acumulare sau pierdere). Această variaţie are loc datorită influenţelor externe (precipitaţii, evaporaţie, radiaţii), consumurilor locale (necesarul prelevat de rădăcini) sau schimburilor cu alte faze (îngheţ, evaporaţie, condensare).
Variaţiile stocului sunt descrise cantitativ prin legea conservării materiei (ecuaţia de continuitate). Deci descrierea globală a transferurilor se obţine prin asocierea unei legi dinamice cu ecuaţia de continuitate. 3.3.1. Legea dinamică exprimă faptul că mişcarea ( fluxul ) rezultă din acţiunea unei forţe motrice ( gradient de potenţial ). φ⋅= gradKJ
ρ (3.4)
unde J - flux sau densitate de flux; K - coeficient de transfer; φ - potenţial;
52 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
gradφ - forţa motrice care provoacă mişcarea.
Legile dinamice utilizate în mecanica mediilor poroase sunt:
• Legea lui Darcy care exprimă faptul că fluxul de apă este proporţional cu gradientul de potenţial hidraulic;
• Legea lui Fourier exprimă proporţionalitatea dintre fluxul de căldură şi temperatură ;
• Legea lui Fick traduce proporţionalitatea dintre fluxul de gaz sau de soluţie şi gradientul de concentraţie.
Se constată experimental că mişcarea apei într-un mediu poros poate fi produsă de existenţa unor gradienţi (diferiţi de gradientul de sarcină). Astfel, apa se deplaseză dinspre zona cu voltaj ridicat spre cea cu voltaj scăzut. Acest principiu a fost folosit pentru drenajul electrocinetic al solurilor puţin permeabile (Terzaghi şi Peck 1967). De asemenea apa se delpasează din zonele cu concentraţie mare spre cele cu concentraţie mică şi din zonele cu temperatură mare spre cele cu temperatură redusă. Se poate scrie o lege dinamică generalizată sub forma (3.5). Viteza în mediul poros este:
gradTKgradCKgradEKgradhKU 4321 −−−−=
ρ (3.5)
Coeficienţii de transfer iK pot fi scalari sau tensori., h este potenţialul
hidraulic, E este potenţialul electric, C este concentraţia, iar T este temperatura. Similar, alte fluxuri în mediul poros (flux de electricitate, de elemente în soluţie, de căldură) vor fi legate de aceşti gradienţi prin alţi coeficienţi. 3.3.2.Legea de conservare a materiei se exprimă sub forma :
iQJdivtE ∑+⋅−=
∂∂ ρ
(3.6)
unde E - concentaraţia volumică de element considerat ; J - fluxul sau densitatea de flux ; Qi - aporturi sau prelevări din sistem.
Ecuaţia de continuitate în mediul poros.
Ecuaţia de continuitate pentru curgerea unui fluid printr-un mediu poros este o ecuaţie diferenţială, cu derivate parţiale, care exprimă conservarea masei.
Dispersia poluanşilor în apele subterane 53 Fie v viteza reală a fluidului în porii mediului poros (viteză microscopică) şi ρ densitatea fluidului la această scară. Notăm cu n porozitatea punctuală (n = 1 într-un por şi n = 0 în particula solidă). Cantităţile microscopice sau medii în mediul poros ⟩⟨⟩ρ⟨⟩⟨ n,,v se pot defini fie prin integrare în spaţiu printr-o convoluţie (cu o funcţie de pondere m), fie printr-o definire probabilistică (prin speranţa matematică a marimilor v , ρ, n în punctul x considerat, pentru ansamblul de realizări posibile ale mediului). Ecuaţia de continuitate pentru mediul poros va fi [de Marsily 1981]:
[ ] [ ] 0=⟩⟨⋅⟩ρ⟨∂∂
+⟩⟨⋅⟩ρ⟨ nt
vdiv ρ (3.7)
unde < > are semnificaţia de medie . Această ecuaţie arată că într-un volum fix, suma fluxurilor masice care intră sau ies este egală cu variaţia masei conţinută în acest volum. Deşi se exprimă punctual, legea se stabileşte pentru un volum elementar, fix în spaţiu. Formula Ostrogradski, pentru o viteză oarecare se scrie:
∫∫∫ ∫∫ ⋅⋅−=⋅D S
dSnvdVvdiv ρρρ (3.8)
unde nρ este normala exterioară la S . Prin aplicarea formulei Ostrogradski mărimii [< ρ > < v >] este evident că
div [< ρ > < v >] reprezintă fluxul masic elementar care traversează suprafaţa dS a domeniului D. Viteza < v
ρ
> este o viteză medie fictivă, astfel integrarea se face pe toată suprafaţa S a domeniului D (nu numai a porilor). Masa de fluid conţinută în D nu este dV ci n dV pentru că fluidul este doar în pori (n este porozitatea
totală).
∫∫∫ρD
∫∫∫D
În continuare vom nota < v > = Uρ
viteza macroscopică medie în mediul poros sau viteza de filtraţie, < ρ >= ρ şi < n >= n. • Ecuaţia de continuitate macroscopică pentru cazul general ( fără sursă ) va fi :
[ ] [ ] 0=⋅ρ∂∂
+⋅ρ nt
Udivρ
(3.9)
• Dacă mediul poros este alimentat de o sursă exterioară, debitul masic primit sau
cedat din exterior fiind ( ρq ), ecuaţia de continuitate se scrie:
[ ] [ ] 0=⋅ρ±⋅ρ⋅∂∂
+⋅ρ qnt
Udivρ
(3.10)
54 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov • Dacă fluidul este incompresibil şi scheletul solid nedeformabil ecuaţia de
continuitate devine : 0=Udiv
ρ (3.11)
• Dacă fluidul este compresibil şi curgerea este permanentă 0=∂∂ρ
t,iar ecuaţia de
continuitate devine [ ] 0=⋅ρ Udiv
ρ. (3.12)
În cazul în care mediul poros este izotrop din punct de vedere al repartiţiei
porozitaţii, viteza microscopică medie, reală, într-o secţiune, este dată de relaţia: ρ
CnUvρ= , (3.13)
Viteza medie reală în pori este mai mare decât viteza de filtraţie (nC<1 este porozitatea cinematică). În cazul în care mediul este anizotrop se defineşte o porozitate cinematică de suprafaţă:
tiuniisecatotalarafatasup
eficaceporilorrafatasupnCS ⋅⋅⋅⋅⋅
= (3.14)
iar viteza real în pori va fi : ρă
SCn
Uvρ= (3.15)
Curgerea într-un mediu poros poate fi:
• uniformă (caracteristicile curgerii sunt invariabile în timp şi spaţiu), • permanentă ( constantă în timp ), • nepermanentă. • Din punct de vedere al regimului vitezelor, curgerea poate fi : o laminară (curgerea este lentă şi se desfăşoară în straturi paralele, fără
amestec de masă şi energie între ele), o turbulentă (curgere cu viteze mari, având loc transferul de masă şi energie
între straturi). Stabilirea regimului de curgere se face pe baza numărului Reynolds corespunzător:
υ⋅
=dURe (3.16)
U - viteza medie a apei (m/s),
Dispersia poluanşilor în apele subterane 55 d - diametrul porilor ( m ), υ - vâscozitatea cinematică ( m2/s). Dacă Re < 10 mişcarea este laminară, Re > 10 mişcarea este turbulentă. 3.4. CIRCULATIA APEI ÎNTR -UN MEDIU POROS SATURAT 3.4.1. Legea lui Darcy
Vom analiza curgerea într-un mediu poros saturat de lungime ∆x şi secţiune S , prin care curge un debit volumic Q,, constant în timp. O astfel de curgere poate fi realizată într-o instalaţie ca cea din figura 3.1.
Vom defini curgerea în mediul poros printr-un vector “ flux de curgere “ care
este debitul specific, q = QS
, sau viteza medie Darcy Uρ
. Această mărime reprezintă
media globală a fluxurilor microscopice într-un volum de sol suficient de mare în comparaţie cu dimensiunile porilor şi cu eterogenităţile microscopice. Darcy a stabilit experimental relaţia dintre debitul Q ( m3/s ) ce stăbate proba şi denivelarea ∆H (m) dintre cele două rezervoare :
xHSKQ
∆∆
⋅⋅−= . (3.17)
Relaţia (3.17) reprezintă legea lui Darcy pentru un mediu poros saturat. În această ecuaţie mărimile au următoarea semnificaţie:
∆H = H2-H1 = (z2 + h2 ) - (z1 + h1 ) < 0 (3.18) ∆H reprezintă pierderea de sarcină la traversarea probei,
IxH
=∆∆ . (3.19)
I se numeşte pantă hidraulică, respectiv pierderea de sarcină pe unitatea de lungime, în direcţia de curgere sau gradientul hidraulic şi reprezintă forţa motrice ce provoacă mişcarea. ∆x = x2 - x1>0, (3.20) K este conductivitatea hidraulică sau permeabilitatea hidrogeologică a mediului poros ([K]SI = LT-1),
UqSQ
== ( m/s ) (3.21)
56 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
Debitul specific, q ( debit prin unitatea de suprafaţă sau flux ) reprezintă volumul de apă scurs prin unitatea de suprafaţă în unitatea de timp. Acest flux are dimensiunile unei viteze (este viteza fictivă pe care ar avea-o apa dacă ar traversa toată suprafaţa S a solului). Viteza medie reală, microscopică, în pori va fi: ρ
CC nq
nUv
ρρ== , (3.22)
nC fiind porozitatea cinematică sau eficace. În literatură U se numeşte viteza Darcy sau viteza de filtrare iar q flux sau debit specific. 1 2 z2
h2 h1
z
∆H
H2 H1
z2
z1 ∆x
x1
S
Figura 3.1 Deducerea legii Darcy
Legea lui Darcy se poate scrie sub forma diferenţială:
sHKq
dd
−= sau dsdHK−=U ,
s fiind o direcţie oarecare . Într-un sistem tridimensional : gradHKq ⋅−=
ρ sau U gradHK ⋅−=ρ
x
x2(3.23)
(3.24)
Dispersia poluanşilor în apele subterane 57 Această lege arată că mişcarea apei se face în direcţia forţei motrice reprezentată de gradientul hidraulic, fluxul q fiind un vector perpendicular pe liniile echipotenţiale (H = ct.)
kzHj
yHi
xHH)k
zj
yi
x(HgradH
ρρρρρρ∂∂
+∂∂
+∂∂
=∂∂
+∂∂
+∂∂
=∇= . (3.25)
Conductivitatea hidraulică, K, este un tensor. 3.4.2. Limite de valabilitate ale legii lui Darcy
Legea lui Darcy este valabilă pentru regimurile de curgere laminară care au loc, de obicei, în nisipurile fine, silţuri şi argile. În nisipurile grosiere şi pietrişuri, vitezele cresc şi regimul devine turbulent. În acest caz relaţia dintre flux şi gradientul sarcinii nu mai este liniară ci de forma:
2UUgradH ⋅β+⋅α= . (3.26) αU reprezintă pierderile de sarcină datorate frecării vâscoase la pereţii matricei solide iar βU2 , pierderile datorate inerţiei fluidului (disipaţii de energie cinetică în pori, asemănătoare celor care apar la îngustarea unui tub). Se defineşte numărul Reynolds al mediului poros, adimensional:
υ⋅
=η
ρ⋅⋅≈
ηρ⋅⋅
=dUdUkURe (3.27)
unde mărimile au următoarea semnificaţie: U - viteza de filtrare ( m/s ) ;
k - rădăcina pătrată a permeabilităţii intrinsece (m) ; ρ - densitatea fluidului (kg/m3); η - vâscozitatea dinamică a fluidului (kg/ms) ; ν - vâscozitatea cinematică a fluidului (m2/s) ; d - diametrul mediu al particulelor sau diametrul eficace d10 ( m ). În practică se admite că legea lui Darcy este valabilă pentru numere Re mai mici decât o limită cuprinsă între 1 şi 10. În acest caz curgerea este pur laminară în interiorul porilor. Între 10 şi 100 începe un regim de tranziţie în care forţele de inerţie nu mai sunt neglijabile şi unde legea lui Darcy nu se mai aplică. Pentru Re > 100 regimul devine turbulent iar relaţia lui Darcy trebuie înlocuită cu o relaţie de forma (3.26). În practică curgerea rămâne laminară în majoritatea cazurilor de curgere în medii poroase, excepţie făcând regimul carstic şi zona din imediata apropiere a lucrărilor de captare. Sichardt recomandă o valoare limită pentru gradientul hidraulic (până la care este valabilă legea lui Darcy), în funcţie de conductivitatea hidraulică, K (m/s) :
58 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
KI
151
= (3.28)
Limita inferioară de valabilitate variază mult cu tipul de argilă. Astfel când vitezele sunt foarte mici ele nu mai sunt proporţionale cu gradientul sarcinii. Forţele de adsorbţie sunt predominante şi legea lui Darcy nu mai este valabilă. 3.4.3. Conductivitatea hidraulică şi permeabilitatea intrinsecă Conductivitatea hidraulică la saturaţie K, numită şi permeabilitatea hidrogeologilor, caracterizează posibilitatea solului de a lăsa să circule apa prin el.
dsdHqK −= (3.29)
Conductivitatea hidraulică este influenţată atât de proprietăţile mediului poros cât şi de cele ale fluidului. Un sol grosier (pietriş, nisip) lasă să circule apa mai uşor decât un sol argilos. Circulaţia apei va fi influenţată de structura solului şi de distribuţia porilor. Astfel influenţa mediului poros se defineşte printr-o mărime numită permeabilitate intrinsecă (ki). Această mărime se măsoară în (m2 ) şi reprezintă capacitatea unui mediu poros de a lăsa să circule un fluid oarecare. Ea este definită la scară macroscopică.
Dacă considerăm că adevăratele cauze ale deplasării unui fluid într-un mediu poros sunt gradienţii de presiune şi forţele exterioare (gravitaţionale în cazul de faţă), ( H = p/ρg + z), legea lui Darcy se poate exprima sub forma generală :
( gradzggradpk
U i ⋅⋅ρ+⋅η
−=ρ
) . (3.30)
qU ρρ= este viteza Darcy sau debitul specific şi este o mărime macroscopică iar
η, ρ, p vor fi valorile medii ale vâscozităţii, densităţii şi presiunii, ki este permeabilitatea intrinsecă, iar; η - vâscozitatea dinamică a fluidului . Dimensiunea permeabilităţii intrinseci este : [ki] = L2
[ ] [ ] [ ][ ] [ ]
2222
1113
1 LTLML
TLMTLLpS
Qki =⋅⋅⋅
⋅⋅⋅⋅=
⋅⋅η⋅
=−−
−−−
− (3.31)
Permeabilitatea intrinsecă se măsoară în DARCE = 10-12 m2 sau în
DARCY=0,987*10-12 m2. 1 DARCY este permeabilitatea unui mediu care sub diferenţa de presiune de
1 At (760 mm Hg ) pe un cm de lungime, lasă să curgă printr-o suprafaţă de1 cm2 un
Dispersia poluanşilor în apele subterane 59 debit de 1 cm3/s, pentru un fluid cu vâscozitatea dinamică de 1 centipoise (Bear, 1972 ). 1 MILIDARCY = 10-3DARCY
Presupunând fluidul incompresibil, relaţia (3.30) devine :
)zgp(gradk
U i ⋅⋅ρ+⋅η
−=ρ
. (3.32)
Dacă definim sarcina hidraulică totală (se neglijează termenul cinetic din cauza vitezelor foarte mici) ca fiind suma dintre poziţia faţă de un sistem de referinţă şi sarcina de presiune:
zg
pH +⋅ρ
= , (3.33)
gradHKgradHgk
U i ⋅−=⋅η
⋅ρ⋅−=
ρ (3.34)
Mărimea η
⋅ρ⋅=
gkK i se numeşte conductivitate hidraulică sau
permeabilitatea hidrogeologilor şi ţine seama atât de permeabilitatea intrinsecă a mediului (ki) cât şi de natura fluidului (densitate şi vâscositate). Dimensiunea conductivităţii hidraulice, în SI, este (m/s).
[ ] 111
232−
−−
−−
⋅=⋅⋅
⋅⋅⋅⋅= TL
TLMTLLMLK , [ K ]SI = m/s. (3.35)
Conductivitatea hidraulică depinde de temperatură (vâscositatea depinde de temperatură)
)t(K)t()t()t(K 1
2
12 ⋅
ηη
= (3.36)
K(t2) este corespunzătoare temperaturii t2 iar K(t1) temperaturii t1 . La temperatura de 200 C, petru o permeabilitate intrinsecă de 1 milidarcy permeabilitatea hidrogeologică este:
83
31510966,0
10002,181,91010987,0 −
−
−⋅=
⋅
⋅⋅⋅ m2/s
În tabelul 2.1 sunt date câteva valori ale conductivităţii hidraulice pentru diferite tipuri de soluri iar în tabelul 2.2 pentru diferite tipuri de roci. Materialele consolidate ( gresii, roci diverse, elemente carbonatate ) au valori ale lui K variabile în funcţie de porozitatea fisurală ( fisuri, canale de alteraţie sau de dizolvare a rocilor carbonatate ).
60 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
Tabelul 3.1. Conductivitatea hidraulică pentru diferite tipuri de sol
Natura solului Conductivitate hidraulică m/s
Conductivitate hidraulică m/zi
Sol argilos 10-7 - 10-6 0,01 - 0,1 Sol aluvionar de suprafaţă 10-6 - 10-5 0,1 - 1 Nisip fin 10-5 - 5*10-5 1 - 5 Nisip mediu 5*10-5 - 2,5*10-3 5 - 20 Nisip grosier 2,5*10-5 - 10-3 20 - 100 Pietriş > 10-3 > 100
Straturile de nisip sedimentare sau argilo-nisipoase au, datorită stratificaţiei o permeabilitate orizontală mai mare decât cea verticală. Mediile aluvionare sunt formate din straturi sau lentile alternative de nisip, pietriş şi argile. Pentru aceste medii curgerea va avea tendinţa de a urma direcţia cu permeabilitatea cea mai mare. Conductivitatea hidraulică trebuie considerată o proprietate tensorială. Se defineşte un tensor de ordinul doi prin regula transformării componentelor tensorului printr-o rotaţie a sistemului de coordonate. Într-un mediu stratificat, direcţiile paralele şi perpendiculare pe stratificaţii sunt direcţii principale ale curgerii, pentru care componentele tensorului se reduc la componentele diagonale. Dacă o matrice este simetrică , valorile sale proprii sunt distincte iar direcţiile proprii sunt ortogonale. Permeabilitatea intrinsecă, ik este o matrice cu 9 coeficienţi :
zzzyzx
yzyyyx
xzxyxx
i
kkkkkkkkk
k = cu kxy=kyx , kxz = kzx , kyz = kzy (3.37)
Dacă ]gradzggradp[ki ⋅⋅ρ+⋅η
−=ρ
U , componentele vitezei vor fi :
)gzp(
kypk
xpk
U xzxyxxx ⋅ρ+
∂∂
⋅η
−∂∂
⋅η
−∂∂
⋅η
−=
)gzp(
kypk
xpk
U yzyyxyy ⋅ρ+
∂∂
⋅η
−∂∂
⋅η
−∂∂
⋅η
−= (3.38)
)gzp(
kypk
xpk
U zzzyxzz ⋅ρ+
∂∂
⋅η
−∂∂
⋅η
−∂∂
⋅η
−=
Dispersia poluanşilor în apele subterane 61 Se vor deduce axele X,Y,Z, din primele, printr-o rotaţie astfel ca tensorul de permeabilitate să se reducă la componentele diagonalei principale. Matematic X,Y,Z,sunt direcţiile vectorilor proprii ai matricei ik . Fizic X,Y,Z sunt direcţiile pentru care curgerea este efectiv paralelă cu gradientul sarcinii ( în practică o direcţie este ortogonală la stratificaţie şi două paralele cu aceasta ). Aceste direcţii sunt numite direcţii principale de anizotropie . Tensorul ik devine :
zz
yy
xx
i
kk
kk
000000
= (3.39)
iar vitezele : xpk
U xxx ∂
∂⋅
η−=
ypk yy
y ∂∂
⋅η
−=U (3.40)
)gzp(
kzzz ⋅ρ+
∂∂
⋅η
−=U
În practică, în mediile sedimentare cu stratificaţii orizontale, se observă două permeabilităţi ( una verticală kzz şi una orizontală kxx=kyy ). Raportul de anizotropie kxx/kzz este în general cuprins între 1 şi 100. Un mediu poros este omogen atunci când permeabilitatea intrinsecă, într-o direcţie dată, este constantă. Dacă acest coeficient rămâne constant, oricare ar fi direcţiile la care ne referim, mediul se cheamă omogen şi izotrop; în caz contrar mediul este anizotrop (eterogen). În general mediile poroase naturale sunt neomogene şi anizotrope, permeabilitatea lor variind de la un punct la altul şi având în acelaşi timp proprietăţi direcţionale . Dacă raportul de anizotropie rămâne constant în tot domeniul mişcării mediul este ortotrop.
Având în vedere egalitatea : η
⋅ρ⋅=
gkK i aceleaşi aprecieri pot fi făcute în
legătură cu conductivitatea hidraulică a mediului poros. Diferenţierea între rocile permeabile şi impermeabile se face în mod arbitrar la o valoare a conductivităţii hidraulice de 10-9 m/s. Argilele sunt impermeabile în ciuda porozităţii totale mari (din cauza dimensiunilor mici ale porilor, porozitatea eficace este mică ).
62 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov Gresiile au o permeabilitate analoagă nisipului dacă nu sunt cimentate. Dacă gresiile sunt cimentate cu calcar acesta poate fi dizolvat de apele ce conţin CO2 şi permeabilitatea creşte. Tabelul 3.2 Conductivitatea hidraulică pentru diferite tipuri de roci Tipul de rocă sau de material
Conductivitate hidraulică K m/zi
Permeabilitate intrinsecă ki (m2)
Argile 10-7 - 10-3 10-19 - 10-15 Mâluri aluvionare (silţuri) 10-4 - 100 10-16- 10-12 Nisipuri 10-3 – 102 10-14 - 10-9 Pietrişuri 10-2 – 102 10-10 - 10-7 Şisturi argiloase , marne 10-8 - 10-4 10-20 - 10-16 Marne fracturate şi erodate
10-4 - 10-0 10-16 - 10-12
Gresie bine cimentată 10-5 - 10-2 10-17 - 10-14 Gresie friabilă 10-3 - 10-0 10-15 - 10-12 Sare 10-10 - 10-8 10-22 - 10-20 Anhidrite 10-7 - 10-6 10-19 - 10-18 Roci metamorfice nefracturate
10-9 - 10-5 10-21 - 10-17
Roci metamorfice fracturate
10-5 - 10-1 10-17 - 10-13
O unitate hidrogeologică este omogenă dacă proprietăţile sale hidraulice sunt
aceleaşi în orice punct. Eterogenitatea unei zone depinde de scara la care este analizat fenomenul. Se pune problema stabilirii unei valori medii a conductivitaţii (conductivitate hidraulică efectivă Ke) astfel încât să poată fi folosită în modelele aproximative. Pentru curgerea permanentă, cu un gradient hidraulic spaţial uniform, se folosesc următoarele reguli de mediere: 1. Mediu perfect stratificat ( N straturi de grosime li şi conductivitate hidraulică Ki
• Pentru curgerea paralelă cu stratul
∑∑=
=
⋅=
N
iN
ii
iie
l
KlK
1
1
(3.41)
• Pentru curgerea perpendiculară pe strat :
∑∑=
=
=N
iN
i i
i
ie
Kl
lK
1
1
(3.42)
Dispersia poluanşilor în apele subterane 63 2. Mediu eterogen , nestratificat , în care s-au făcut m măsurători :
• Modele bidimensionale Ke = KG = ( K1⋅K2⋅…Km) 1/m (3.43) • Modele tridimensionale Ke = KG = ( 1+σy
2 / 6 ) (3.44) unde σy
2 este varianţa logaritmilor naturali ai măsurătorilor conductivităţii. Dacă gradientul hidraulic nu este constant nu există reguli de mediere a
conductivităţii. 3.4.4. Transmisivitatea
z
O
M
Ux
x
Figura 3.2. Secţiune transversală prin stratul acvifer de grosime l
Dacă apa subterană circulă într-un strat de grosime l şi dacă dorim să calculăm fluxul printr-o suprafată transversală, pe direcţia de curgere, acesta este :
∫ ∫ ⋅=⋅⋅=l l
x dzUdznUlQ
0 0
ρρ (3.45)
nρ - normala la oz Ux- componenta vitezei în direcţia x. Presupunând că z este direcţia principală de anizotropie ( celelalte direcţii x,y sunt în planul stratului ), atunci în toate punctele M ale lui oz.
gradHKU M ⋅−=
ρ (3.46)
MK este tensorul conductivităţii în planul xy care trece prin M şi grad H este gradientul sarcinii în acest plan . Presupunând că grad H este constant pe direcţia Oz:
64 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
∫ ==l
M itatetransmisivTdzK0
(3.47)
Dacă mediul este izotrop ( K =constantă după direcţia z ) transmisivitatea este egală cu produsul dintre conductivitatea hidraulică şi grosimea stratului acvifer: T=K⋅l ( m2/s ) (3.48) 3.5. MODELAREA MATEMATICA A CURGERII IN ACVIFERE 3.5.1. Clasificarea acviferelor
Prin acvifer se înţelege o formaţiune geologică care conţine apă şi care permite, în condiţii normale, circulaţia unei cantităţi semnificative de lichid. (aqua = apă, ferre = a purta, a duce, phreatos = puţ). a) Pânza freatică sau acviferul cu suprafaţă liberă
Presupunem un sol poros, uniform şi permeabil. Apa de ploaie se infiltrează
şi saturează roca poroasă până la un anumit nivel, numit suprafaţă liberă. Numim pânza freatică sau acvifer freatic, zona saturată, aflată între
suprafaţa liberă şi stratul impermeabil de la bază. In pânza freatică apa circulă spre zona de izvorâre care cuprinde punctele
de cotă minimă ale topografiei (izvoare, râuri din reţeaua hidrografică de supafaţă).
Figura. 3.3. Schematizarea unei pânze freatice “de vale” In Figura 3.3 sunt reprezentate liniile de curent şi liniile de egală sarcină (echipotenţiale, curbe izopieze, curbe piezometrice, hidroizohipse).
Dispersia poluanşilor în apele subterane 65
Dacă permeabilitatea este izotropă, liniile de curent sunt ortogonale la liniile echopotenţiale.
Panta suprafeţei libere indică sensul circulaţiei în pânză. Apa circulă pe toată grosimea acviferului, vitezele fiind mai mari la suprafaţă decât la adâncime, întrucât traiectoriile sunt mai scurte pentru aceeiaşi diferenţă de sarcină.
In cursul anului, suprafaţa liberă a pânzei oscilează între valori maxime (aprilie-mai) şi valori minime (octombrie-noiembrie) datorită timpului de parcurgere a zonei nesaturate, de către apa rezultată din ploaie. In multe cazuri panta suprafeţei libere este mică, echipotenţialele sunt practic verticale, iar liniile de curent sunt practic paralele cu suprafaţa liberă (aproape orizontale). Excepţie face zona din imediata apropiere a regiunilor de izvorâre şi de alimentare (Figura 3.4).
Figura. 3.4. Pânză freatică cu curgere aproape orizontală
Pânza freatică din zonele aride are un carácter special. In cazul în care alimentarea pânzei din precipitaţii este slabă, suprafaţa liberă a pânzei este coborâtă iar aspectul suprafetei libere este ca cel din figura 3.5. Alimentarea pânzei se face din apele de suprafaţă.
Figura. 3.5. Pânză freatică alimentată din apele de suprafaţă
66 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
Pânzele aluvionare se desfaşoară în câmpiile aluvionare. Materialul aluvionar depus de râu (nisip, pietriş), este foarte permebil, legătura dintre apele de suprafaţă şi apele subterane fiind foarte puternică. În Figura 3.6 se observă că, pe anumite porţiuni, apa subterană alimentează pânza freatică. Suprafaţa liberă a pânzei freatice este foarte apropiată de suprafaţa solului şi de multe ori deasupra nivelului râului. Astfel de pânze se mai numesc şi subfluviale.
Legătura dintre râu şi pânza freatică poate fi împiedicată prin colmatarea fundului râului sau lacului.
Figura 3.6. Pânză freatică aluvionară
Dispersia poluanşilor în apele subterane 67
Pânzele nealimentate de râu şi suprapuse au la bază un strat impermeabil (de exemplu marne, argile) şi sunt cantonate într-un mediu poros foarte permeabil. Curgerea se face spre punctele cele mai coborâte, unde apar zone de izvorâre. Se poate întâmpla ca, sub stratul impermeabil să existe un alt strat permeabil (de exemplu, calcar) în care să se formeze o altă pânză freatică. Între cele două pânze pot exista legături verticale, prin drenanţă.
Figura 3.7. Pânze freatice suprapuse b. Pânza captivă sau acviferul sub presiune.
Un acvifer sub presiune (confinat) este limitat deasupra şi dedesubt prin formaţiuni impermeabile. Intr-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
Figura 3.8. Acviferul sub presiune
68 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
Nivelul apei dintr-un număr infinit de astfel de puţuri de observaţie reprezintă o suprafaţă imaginară, numită suprafaţă piezometrică.
Un acvifer artezian este un acvifer sub presiune pentru care nivelul suprafeţei piezometrice este deasupra suprafeţei solului. Un puţ practicat într-un astfel de acvifer va permite curgerea apei în sus, spre suprafaţa solului, fără pompare. Numele vine de la localitatea Artesia din nordul Franţei unde, in secolul XII, au fost construite astfel de puţuri.
Două acvifere pot comunica între ele, pe anumite porţiuni din suprafaţa lor, prin intermediul straturilor semimpermeabile ce le despart – (leaky acvifer). Astfel, pot exista acvifere sub presiune, alimentate prin drenanţă, dintr-un strat cu presiune mai mare.
Figura 3.9. Analogia între o pânză captivă şi un permeametru “U” c) Tipuri de acvifere în funcţie de condiţiile de margine.
Atât pentru acviferele cu suprafaţă liberă cât şi pentru acviferele sub presiune, se pot defini urmatoarele tipuri:
1. acvifere infinite, la care zona de alimentare se află la mare distanţă de
puţul de extracţie, 2. acvifere semiinfinite, la care una din zonele de alimentare se afla la
mare distanta iar in apropiere de puţul de extracţie se află o limită impermeabilă sau o altă zonă de alimentare,
Dispersia poluanşilor în apele subterane 69
3. acvifere tip bandă, la care zonele de alimentare sau limitele impermeabile se afla la mică distanţă de puţul de extracţie şi sunt în general paralele,
4. acvifere unghiulare, la care limitele amintite mai sus formează un unghi oarecare, cunoscut,
5. acvifere închise pe contur, la care limita de alimentare sau limita impermeabilă este reprezentată de un contur închis, bine determinat, cunoscut.
3.5.2. Rolurile unui acvifer
Sursă de apă reînnoibilă. Un acvifer poate fi alimentat din precipitaţii în funcţie de distribuţia şi de intensitatea ploii, de topografia terenului, de acoperirea cu vegetaţie şi de permeabilitatea solului. Apele de suprafaţă pot constitui condiţii de frontieră pentru apele subterane. Ele pot alimenta freaticul sau pot fi alimentate de acesta.
Există acvifere de adâncime, în care apa provine din înmagazinarea în ere îndepărtate, în condiţii climatice diferite de cele actuale. Acestea reprezintă surse nereînnoibile.
Rezervor de înmagazinare. Mari cantităţi de apă pot fi stocate într-un
acvifer freatic, folosind tehnici de încărcare artificială (pentru perioade mai lungi sau mai scurte).
Mediu conductor. Apa poate fi introdusă într-un punct din freatic şi captată,
prin pompare, într-un alt punct. Apa injectată va curge dinspre zona cu nivel ridicat, din zona de încărcare, spre puţul de pompare (cu nivel scăzut).
Filtru. Prin tehnica de încărcare artificială, un acvifer poate fi folosit ca filtru şi purificator pentru apa injectată, încărcată cu impurităţi. Astfel filtrarea va consta din: • reţinerea suspensiilor; • reacţii chimice (adsorbţia şi schimburi ionice la suprafaţa materiei solide, în
special în solurile argiloase); • amestecul apei poluate, injectate, cu apa din acvifer (datorită geometriei
traiectoriilor şi mecanismului dispersiei hidrodinamice).
Controlul curgerii de bază, se poate realiza (în izvoare şi râuri) prin controlul nivelului apei în acviferele care le alimentează.
70 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
3.5.3 Modelarea matematică a curgerii apei în acvifere Un model poate fi definit ca o versiune simplificată a sistemului real care
simulează aproximativ relaţia “excitaţie-răspuns” a acestuia. Simplificările sunt introduse sub forma unui set de ipoteze care exprimă
înţelegerea noastră privind natura sistemului şi comportarea sa. Ca urmare, nu va exista un model unic pentru un sistem acvifer dat. Fiecare set de ipoteze va duce la un alt model. Astfel, pot fi definite modelele matematice numerice, fizice, analogice.
Alegerea celui mai potrivit model, pentru un caz dat, depinde de obiectivele investigaţiei şi de resursele disponibile (timp, buget, computere).
Majoritatea modelelor exprimă bilanţul, cantităţilor considerate (masa de apă, masa de poluant, căldura).
Primul pas de modelare este construirea unui model conceptual al domeniului acvifer. Acesta constă dintr-un set de ipoteze care reduc problema reală şi domeniul real, la o versiune simplificată, acceptabilă din punct de vedere al obiectivelor şi al problemelor de management asociate.
Modelul conceptual se poate exprima printr-o formă matematică, numită model matematic. Acesta trebuie să conţină:
• Definirea geometriei domeniului considerat şi a frontierelor; • Ecuaţiile care exprimă bilanţul mărimilor considerate; • Ecuaţiile care descriu fluxurile cantităţilor extensive considerate (variabilele de
stare relevante); • Ecuaţiile constitutive care definesc comportarea materialelor (fluid, solid); • Condiţiile iniţiale care descriu starea sistemului la momentul iniţial; • Condiţiile pe frontiere, care descriu interacţiunea acviferului cu mediul
înconjurător.
Odată construit modelul matematic, în funcţie de variabilele de stare relevante, acesta trebuie rezolvat pentru cazuri de interes practic. (Sunt preferate metodele analitice, acestea dând soluţii general valabile). De cele mai multe ori metodele analitice nu pot fi folosite din cauza dificultăţii modelului. În acest caz se folosesc metode numerice, pentru care:
Dispersia poluanşilor în apele subterane 71 • Soluţia este dată prin valori discrete în timp şi spaţiu (nu prin valori continue); • Ecuaţiile cu derivate parţiale sunt înlocuite printr-un set de ecuaţii algebrice care
conţin valori discrete ale variabilelor de stare; • Soluţia este obţinută folosindu-se un set de valori numerice ale coeficienţilor
caracteristici acviferului.
O etapă importantă a modelării o reprezintă calibrarea modelului. Aceasta presupune estimarea parametrilor caracteristici acviferului astfel încât rezultatele modelului să coincidă cu cele măsurate. Modelul astfel obţinut trebuie validat.
Validarea se face prin compararea valorilor obţinute cu ajutorul modelului cu valorile măsurate sau, pentru cazuri simple, cu cele obţinute analitic. Astfel se pot înlătura erorile rezultate din aproximarea numerică. Modelarea unui caz concret de curgere într-un acvifer va conţine:
1. Descrierea problemei. 2. Obiectivele modelării. 3. Ipoteze de calcul. 4. Modelul conceptual. 5. Ecuaţiile matematice ale modelului. 6. Coeficienţii şi parametrii modelului. 7. Modelul numeric. Programul utilizat. 8. Calibrarea şi estimarea parametrilor. 9. Simulări ale fenomenelor. 10. Studiul sensibilităţii modelului.
72 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
Figura3.10. Schema logică a unui model numeric hidrogeologic (Peck, 1988)
Dispersia poluanşilor în apele subterane 73 3.5.4. Ecuaţia de difuzivitate în 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.
Vom considera o prismă dintr-un acvifer cu suprafaţă liberă, de înălţime h, aflată între substratul impermeabil b(x,z) ş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ă (ipoteza lui Dupuit).
Figura 3.11. Schematizarea unui volum elementar dintr-o pânză freatică
Presupunem că tensorul permeabilităţii admite verticala ca şi 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 ‚şi ieşirilor, astfel:
74 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
a) Fluxul masic care intră în unitatea de timp prin cele două feţe perpendiculare pe
Ox este:
)z,y,dxx(xJ)z,y,x(xJxJ +−= , (3.49)
∫∫ ⋅+⋅⋅ρ−⋅⋅⋅ρ=)t,y,x(H
)y,x(b
)t,y,x(H
)y,x(b
dz)z,y,dxx(U xdydz)z,y,x(U xdyJ x . (3.50)
Ux este componenta vitezei de filtraţie (viteza lui Darcy) după direcţia x. Rezultă
[ ]dxHb dzxU
xdyxJ ∫
∂∂
⋅⋅ρ−= . (3.51)
Conform legii lui Darcy viteza în direcţia x este:
( )xHz,y,xKU xxx ∂
∂−= , (3.52)
iar fluxul masic după direcţia x devine:
⋅
∂∂
⋅∂∂
⋅⋅ρ= ∫H
b xxx dzxHK
xdydxJ (3.53)
(în relaţia anterioară xH
∂∂ nu depinde de z).
În acelaşi fel vom calcula fluxul prin suprafeţele perpendiculare pe Oy:
( )( )
( ) (( )
( )dzz,dyy,xUdxdzz,y,xUdxJ
t,y,xH
y,xb yt,y,xH
y,xb yy ⋅+⋅⋅ρ−⋅⋅⋅ρ= ∫∫ ) (3.54)
dydzUy
dxJH
byy
⋅
∂∂
⋅⋅ρ−= ∫ (3.55)
Conform legii lui Darcy viteza în direcţia y este:
yH)z,y,x(KU yyy ∂
∂−= (3.56)
(unde yH
∂∂ nu depinde de z), iar fluxul masic după direcţia y devine:
Dispersia poluanşilor în apele subterane 75
∂∂
∂∂
⋅⋅⋅ρ= ∫ dzyHK
ydydxJ
H
b yyy (3.57)
Fluxul masic care intră în unitatea de timp prin cele patru feţe perpendiculare pe Ox respectiv pe Oy, este
∂∂
∂∂
+
∂∂
∂∂
⋅⋅ρ=+= ∫∫ dzKyH
ydzK
xH
xdydxJJJ
H
b yyH
b xxyx (3.58)
b) Pentru a se ţ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 ( , w
fiind debitul prelevat pe unitatea de suprafaţă a pânzei freatice
)dydxw ⋅⋅⋅ρ
⋅ tLL2
3
• fie debitul masic Qm de fluid prelevat în element dintr-o sursă exterioară (pozitiv dacă este prelevat şi negativ dacă este injectat), 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:
( ) dydxbHndM d ⋅⋅−⋅ρ= (3.59) nd este porozitatea de drenaj (nu porozitatea totală).
Variaţia acestei mase în timp, va fi:
dydxtHn
dtdM
d ⋅⋅∂∂
⋅ρ= (3.60)
Bilanţul intrărilor şi ieşirilor, ţinând seama de conservarea masei va fi
reprezentat de ecuaţia:
=
∂∂
∂∂
+
∂∂
∂∂
⋅⋅ρ ∫∫ dzKyH
ydzK
xH
xdydx
H
b yyH
b xx
dydxwdydxt
Hnd ⋅⋅⋅ρ±⋅∂
∂⋅ρ= (3.61)
76 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov Cum ρ ≠ 0, dxdy ≠ 0, rezultă
( )( )
( )( )
=⋅±
∂∂
∂∂
+
∂∂
∂∂
∫∫ wdzKyH
ydzK
xH
xt,y,xH
y,xb yyt,y,xH
y,xb xx tHnd ∂
∂= . (3.62)
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 de toată verticala, ecuaţia
difuzivităţii devine:
( ) ( )t
HnwbHKyH
ybHK
xH
x dyyxx ∂∂
=±
−
∂∂
∂∂
+
−
∂∂
∂∂ (3.63)
sau
( ty,x,hb-H =∂
∂=±
∂∂
∂∂
+
∂∂
∂∂ unde,
tHnwhK
yH
yhK
xH
x dyyxx ) . (3.64)
2. Pentru un mediu anizotrop, ţinând seama de definirea transmitivităţii:
h.K)bH(KdzKT xxxx
H
bxxxx =−== ∫ ,
h.K)bH(KdzKT yyyy
H
byyyy =−⋅== ∫
ecuaţia difuzivităţii devine:
tHnw
yHT
yxHT
x dyyxx ∂∂
=±
∂∂
∂∂
+
∂∂
∂∂ . (3.65)
Vom presupune că transmisivitatea variază puţin cu sarcina h, adică variaţiile
lui h sunt neglijabile faţă de (H – b), de exemplu mai mici de 10% sau 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%.
Dispersia poluanşilor în apele subterane 77 3. Dacă mediul este izotrop, transmisivitatea este constantă (Txx = Tyy = T), ecuaţia
difuzivităţii devine:
tH
Tn
Tw
yH
xH d
∂∂
=±
∂∂+
∂∂
2
2
2
2. (3.66)
Această ecuaţie cu derivate parţiale este de tip parabolic şi este asemănătoare cu ecuaţia căldurii. 4. Dacă stratul b(x,y) este orizontal şi luăm b(x,y) = 0, ca plan de referinţă pentru
potenţial, H(x,y,t) – b(x,y) = h(x,y,t) = H(x,y,t), (3.67) iar ecuaţia difuzivităţii devine:
th
nwyhhKyx
hhKx dyyxx ∂∂
=±
∂∂
⋅⋅∂∂
+
∂∂
⋅⋅∂∂ . (3.68)
5. Dacă mediul este izotrop şi uniform, Kxx =K yy =K, ecuaţia difuzivităţii devine:
th
nwyhh
yK
xhh
xK d ∂
∂=±
∂∂
⋅∂∂
⋅+
∂∂
⋅∂∂
⋅ (3.69)
sau
th
Kn
Kw
yhh
yh
yh
xhh
xh
xh d
∂∂
⋅=±∂∂⋅+
∂∂
⋅∂∂
+∂∂⋅+
∂∂
⋅∂∂
2
2
2
2, (3.70)
respectiv
th
Kn
Kwh d
∂∂
=±∇2222 . (3.71)
6. Dacă regimul este permanent, 0=∂∂
th şi ecuaţia difuzivităţii devine::
Kwh 222 ±=∇ . (3.72)
7. Dacă pânza nu este alimentată la suprafaţă (w = 0) şi regimul este permanent, ecuaţia difuzivităţii va fi:
02
2
2
2=
∂∂
+∂∂
yh
xh (3.73)
78 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov 3.5.5. Ecuaţia de difuzivitate în cazul acviferelor sub presiune Ecuaţia de difuzivitate 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 σ. Vom prezenta o deducere a ecuaţiei de difuzivitate, dată de Raudkivi (1976). Deducerea se face pentru un caz particular, al unui mediu izotrop
, iar tasările se fac în mod special după direcţia z (cele orizontale sunt neglijabile în raport cu cele verticale).
KKKK zzyyxx ===
Ecuaţia de continuitate pentru un volum paralelipipedic , având porozitatea “n” , saturat cu apă (cu densitatea ρ), este:
dzdydxdV ⋅⋅=
( ) ( ) ( ) ( )zyxnt
dzdydxUz
Uy
Ux zyx ∆⋅∆⋅∆⋅ρ⋅
∂∂
=⋅⋅
⋅ρ
∂∂
+⋅ρ∂∂
+⋅ρ∂∂
− , (3.74)
unde Ux, Uy, Uz sunt componentele vitezei Darcy. Dacă mediul poros este incompresibil (n=const) şi elementul de volum considerat nu îşi schimbă mărimea:
( ) ( ) ( )t
nUz
Uy
Ux zyx ∂
ρ∂⋅−=⋅ρ
∂∂
+⋅ρ∂∂
+⋅ρ∂∂ . (3.75)
Dacă mediul poros este compresibil, porozitatea mediului poate varia în spaţiu şi timp.
Randkivi (1976) face ipoteza că variaţiile dimensiunilor sunt neglijabile în comparaţie cu . Astfel:
y,x ∆∆z∆
( ) ( ) yxt
zntnzz
tnzyxn
t∆⋅∆⋅
∂ρ∂
⋅∆⋅+∂∂
⋅∆⋅ρ+∆∂∂
⋅ρ⋅=∆⋅∆⋅∆⋅ρ⋅∂∂ (3.76)
Încercăm să evaluăm toţi termenii din ecuaţia (3.76), introducând noţiunile de compresibilitate a solului şi a apei. Astfel se ştie că modulul de elasticitate al solului, Es, este legat de coeficientul de compresibilitate al solului, α [M-1LT2], prin relaţia:
( )
∆∆σ
−=α
=
zzd
dE zz
s1 , (3.77)
Dispersia poluanşilor în apele subterane 79
)
unde σzz este componenta verticală a efortului de presiune intergranulară (ML-1T-2).
( ) ( zzdzzd σ⋅∆⋅α−=∆ , (3.78)
( ) ( )t
zzt
zz
∂σ∂
⋅∆⋅α−=∆∂∂ , (3.79)
( ) ( )t
znzt
n zz
∂σ∂
⋅∆⋅ρ⋅⋅α−=∆∂∂
⋅ρ⋅ . (3.80)
Pentru a evalua termenul tnz
∂∂
⋅∆⋅ρ ţinem seama de ipoteza că volumul
particulelor solide din elementul de volum de mediu poros rămâne constant. Schimbarea volumului elementar se face datorită modificării volumului porilor. ( ) constzyxn =∆⋅∆⋅∆⋅−1 . , (3.81) Dacă , ., constx =∆ consty =∆
( ) ( ) ( )t
ntz
zn
tn zz
∂σ∂
⋅α⋅−−=∂∆∂
⋅∆−
=∂∂ 11 , (3.82)
( ) ( )t
zntnz zz
∂σ∂
⋅∆⋅α⋅−⋅ρ−=∂∂
⋅∆⋅ρ 1 . (3.83)
Pentru a evalua termenul t
zn∂ρ∂
⋅∆⋅ vom ţine seama de conservarea masei
de fluid. constVV =∆⋅ρ=∆⋅ρ 00 . , (3.84)
0=∆⋅ρ+∆⋅ρ )V(dVd , (3.85)
( )
VVdd
∆∆
⋅ρ−=ρ . (3.86)
Dacă se ţine seama de compresibilitatea fluidului:
dpVdV
Ew ⋅−==β
1 , (3.87)
unde p este presiunea fluidului, β coeficientul de compresibilitate al apei, iar modul de elasticitate al apei.
wE
80 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
( ) dpVVdd ⋅β⋅ρ=
∆∆
⋅ρ−=ρ , (3.88)
tp
t ∂∂
⋅ρ⋅β=∂ρ∂ . (3.89)
După cum s-a arătat în [Marinov, 2000, &7.1] efortul total aplicat complexului solid-lichid se descompune în eforturi efective şi presiune neutrală.,
constpzz =+σ . în timp.
Deci tp
tzz
∂∂
−=∂σ∂ . (3.90)
În aceste condiţii termenul din dreapta ecuaţiei (3.76) devine:
( )
( )
( )[ ]
βα
+⋅β⋅∆⋅ρ⋅∂∂
=β⋅+−⋅α+⋅α⋅∆⋅ρ⋅∂∂
=
=∂∂
⋅ρ⋅β⋅⋅∆+∂σ∂
⋅∆⋅α⋅−⋅ρ−∂σ∂
⋅∆⋅ρ⋅⋅α−=
=∆⋅∆
∂ρ∂
⋅∆⋅+∂∂
⋅∆⋅ρ+∂∆∂
ρ⋅
nztpnnnz
tp
tpnz
tzn
tzn
yxt
zntnz
tzn
zzzz
1
1 (3.91)
Ecuaţia de continuitate (3.74) devine:
( ) ( ) ( )
βα
+⋅β⋅ρ⋅∆⋅∆⋅∆⋅∂∂
−=∆⋅∆⋅∆
⋅ρ
∂∂
+⋅ρ∂∂
+⋅ρ∂∂ nzyx
tpzyxU
zU
yU
x zyx . (3.92)
Vom exprima, în ecuaţia de continuitate, presiunea p în funcţie de sarcina
piezometrică
zph +γ
= (z nu este funcţie de x, y, t). (3.93)
Astfel zghgp ⋅⋅ρ−⋅⋅ρ= , (3.94)
( )[ ]g
pgxx
hgzhgxx
hgx
zgx
hgxhg
xp
⋅ρ⋅⋅
∂ρ∂
+∂∂
⋅⋅ρ=−⋅∂ρ∂
+∂∂
⋅⋅ρ=∂ρ∂
⋅⋅−∂ρ∂
⋅⋅+∂∂
⋅⋅ρ=∂∂ (3.95)
ρ⋅
∂ρ∂
+∂∂
⋅⋅ρ=∂∂ p
yyhg
yp , (3.96)
( )ρ
⋅∂ρ∂
+
−
∂∂
⋅⋅ρ=∂ρ∂
⋅−⋅+
−
∂∂
⋅⋅ρ=∂∂ p
zzhg
zzghg
zhg
zp 11 , (3.97)
Dispersia poluanşilor în apele subterane 81
ρ⋅
∂ρ∂
+∂∂
⋅⋅ρ=∂∂ p
tthg
tp , (3.98)
dpd ⋅ρ⋅β=ρ , (3.99)
tp
t ∂∂
⋅ρ⋅β=∂ρ∂ , (3.100)
ρ
⋅∂ρ∂
+∂∂
⋅⋅ρ⋅ρ⋅β=∂∂
⋅ρ⋅β=∂ρ∂ p
xxhg
xp
x , (3.101)
∂ρ∂
⋅ρ
+∂∂
⋅⋅ρ⋅ρ⋅β=∂∂
⋅ρ⋅β=∂ρ∂
yp
yhg
yp
y , (3.102)
∂ρ∂
⋅ρ
+
−
∂∂
⋅⋅ρ⋅ρ⋅β=∂∂
⋅ρ⋅β=∂ρ∂
zp
zhg
zp
z1 , (3.103)
∂ρ∂
⋅ρ
+∂∂
⋅⋅ρ⋅ρ⋅β=∂ρ∂
tp
thg
t . (3.104)
Din relaţiile (3.101), (3.102), (3.103), (3.104) rezultă:
xp
xhg
x ∂ρ∂
⋅ρ
⋅ρ⋅β+∂∂
⋅⋅ρ⋅β=∂ρ∂ 2 , (3.105)
( )xhgp
x ∂∂
⋅⋅ρ⋅β=⋅β−⋅∂ρ∂ 21 , (3.106)
xh
pg
x ∂∂
⋅⋅β−⋅ρ⋅β
=∂ρ∂
1
2, (3.107)
yh
pg
y ∂∂
⋅⋅β−⋅ρ⋅β
=∂ρ∂
1
2 , (3.108)
)zh(
pg
z1
1
2−
∂∂
⋅⋅β−⋅ρ⋅β
=∂ρ∂ , (3.109)
82 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov
tp
t ∂∂
⋅ρ⋅β=∂ρ∂ , (3.110)
iar
tpp
thg
tp
thg
tp
∂∂
ρ⋅β⋅ρ
+∂∂
⋅⋅ρ=∂ρ∂
⋅ρ
+∂∂
⋅⋅ρ=∂∂ ,
( )thgp
tp
∂∂
⋅⋅ρ=⋅β−⋅∂∂ 1 ,
th
pg
tp
∂∂
⋅⋅β−
⋅ρ=
∂∂
1 , (3.111)
Dacă se ţine seama de ordinul de mărime al modulului de elasticitate: Ew = 2,07 GPa, respectiv β , 1210 Nm10829,4 −−⋅≈
148010
101
19
9≈
⋅−=
⋅β− p,p şi se poate neglija.
Primul termen al ecuaţiei de continuitate (3.74) devine:
BAz
Uy
Ux
Uz
Uy
Ux
Uzyx
zyx +=∂ρ∂
⋅+∂ρ∂
⋅+∂ρ∂
⋅+
∂
∂+
∂
∂+
∂∂
ρ (3.112)
B A
şi ţinând seama de relaţia Darcy pentru componentele vitezei:
zhKU
yhKU
xhKU
z
y
x
∂∂
⋅−=
∂∂
⋅−=
∂∂
⋅−=
(3.113)
hKzh
yh
xhKA 2
2
2
2
2
2
2∇⋅⋅ρ−=
∂∂
+∂∂
+∂∂
⋅⋅ρ−= , (3.114)
Dispersia poluanşilor în apele subterane 83
∂∂
−
∂∂
+
∂∂
+
∂∂
⋅β−⋅⋅ρ⋅β
−
=
−
∂∂
⋅⋅β−⋅ρ⋅β
⋅∂∂
⋅−
+
∂∂
⋅⋅β−⋅ρ⋅β
⋅∂∂
⋅−+
∂∂
⋅⋅β−⋅ρ⋅β
⋅∂∂
⋅−=
zh
zh
yh
xh
pKg
zh
pg
zhK
yh
pg
yhK
xh
pg
xhKB
2222
2
22
1
11
11
(3.115)
Ecuaţia de continuitate (3.74) devine (3.116), respectiv:
( )
th
pgn
zh
zh
yh
xh
pKg
zh
yh
xhK
∂∂
⋅β⋅−
⋅ρ⋅
βα
+⋅β⋅ρ−=
=
∂∂
−
∂∂
+
∂∂
+
∂∂
⋅β−⋅⋅ρ⋅β
−
∂∂
+∂∂
+∂∂
⋅⋅ρ−
1
1
2222
2
2
2
2
2
2
11
1≈
⋅β− p, iar termenul al doilea din partea stângă a ecuaţiei (3.116) este
neglijabil. Astfel:
thgn
zh
yh
xhK
∂∂
⋅⋅
βα
+⋅β⋅ρ−=
∂∂
+∂∂
+∂∂
⋅⋅ρ− 22
2
2
2
2
2, (3.117)
sau
th
K
gn
zh
yh
xh
∂∂
⋅⋅
βα
+⋅β⋅ρ=
∂∂
+∂∂
+∂∂
2
2
2
2
2
2 . (3.118)
Vom numi coeficient specific de înmagazinare, mărimea Ss având expresia:
[ 1−=
βα
+⋅β⋅⋅ρ LSng s ]- (3.119)
Ecuaţia de difuzivitate pentru cazul unui acvifer sub presiune având K=const. va fi:
th
KS
zh
yh
xh s
∂∂
⋅=∂∂
+∂∂
+∂∂
2
2
2
2
2
2 . (3.120)
84 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov În [Marsily 1981] este dată o demonstraţie pentru cazul general, în care tensorul conductivităţilor hidraulice nu are toate componentele constante. În acest caz ecuaţia de difuzivitate pentru acviferele sub presiune, se scrie:
( ) qthSgradhKdiv s +
∂∂
⋅=⋅ . (3.121)
Se definesc mărimile:
∫ ⋅=b
a xxxx dzKT este transmisivitatea după direcţia x, (3.122)
∫ ⋅=b
a yyyy dzKT este transmisivitatea după direcţia y, (3.123)
∫ ⋅=b
a s dzSS , coeficientul de înmagazinare al acviferului sub presiune (3.124)
∫ ⋅=b
adzqQ , (3.125)
unde a este cota culcuşului (cota stratului impermeabil de la baza acviferului sub presiune), iar b, cota coperişului (cota stratului impermeabil de deasupra acviferului sub presiune).
Presupunând Kxx, Kyy, Ss, constante pe toată înălţimea M acviferului, se pot calcula: Txx = Kxx ⋅M - transmisivitatea după direcţia x , Tyy = Kyy ⋅M - transmisivitatea după direcţia y ,
S = Ss ⋅M =
βα
+⋅β⋅⋅⋅ρ nMg . (3.126)
S este denumit coeficientul de înmagazinare al acviferului sub presiune (adimensional) şi variază între 10-3 şi 10-5.
Presupunând culcuşul şi coperişul paralele (M=const), ecuaţia de difuzivitate devine:
( ) QthSgradhTdiv +
∂∂
⋅=⋅ . (3.127)
Dacă mediul poros este izotrop, transmisivitatea, T este constantă după toate direcţiile, iar ecuaţia de difuzivitate devine:
TQ
th
TS
zh
yh
xh
+∂∂
⋅=∂∂
+∂∂
+∂∂
2
2
2
2
2
2 . (3.128)
Raportul S/T este numit difuzivitatea acviferului, iar Q este aportul de debit din exterior (+Q dacă este injectat sau –Q dacă este pompat).
Dispersia poluanşilor în apele subterane 85
Dacă mişcarea este permanentă
=
∂∂ 0
th şi Q = 0 ecuaţia de difuzivitate
devine:
00 2
2
2
2
2
22 =
∂∂
+∂∂
+∂∂
=∇zh
yh
xhsau,h . (3.129)
3.5.6. Anizotropia
Terenurile anizotrope se pot studia ca şi cele izotrope printr-o schimbare de coordínate. Dacă cele trei direcţii principale de anizotropie sunt x, y, z , iar Kx, Ky, Kz sunt cele trei conductivităţi hidraulice după aceste direcţii, legea lui Darcy se scrie:
xhKU xx ∂
∂⋅−=
yhKU yy ∂
∂⋅−= (3.130)
zhKU zz ∂
∂⋅−=
iar ecuaţia de difuzivitzte este:
thS
zhK
yhK
xhK szyx ∂
∂⋅=
∂∂
⋅+∂∂
⋅+∂∂
⋅ 2
2
2
2
2
2 . (3.131)
Făcând o schimbare de coordonate:
xKKx
x
' ⋅= yKKy
y
' ⋅= zKKz
z
' ⋅= (3.132)
unde K este un coeficient având dimensiunile conductivităţii hidraulice, se obţine
xh
KK
dxdx
xh
xh x
'' ∂∂
⋅=⋅∂∂
=∂∂ , (3.133)
2
2
2
2
xh
KK
dxdx
xxh x
'' ∂∂
⋅=⋅
∂∂
∂∂
=∂∂
xh , (3.134)
iar ecuaţia de difuzivitate se va scrie:
th
KS
zh
yh
xh s
''' ∂∂
⋅=∂∂
+∂∂
+∂∂
2
2
2
2
2
2 . (3.135)
86 Modelarea matematică a curgerii apelor subterane. Anca Marina Marinov Fie o ecuaţie Laplace ordinară în noul sistem de axe. Trebuie remarcat faptul că odată cu anizotropia dispare ortogonalitatea dintre echipotenţiale şi liniile de curent din sistemul de coordonate reale x, y, z, dar aceasta există în sistemul x’, y’, z’. Componentele vitezelor în noul sistem sunt:
''x x
hKU∂∂
⋅−=
''y y
hKU∂∂
⋅−= (3.136)
''z z
hKU∂∂
⋅−=
Se deduce că:
'x
xx U
KK
U ⋅=
'y
yy U
KK
U ⋅= (3.137)
'z
zz U
KKU ⋅=
Dacă se calculează fluxul Q’ al vectorului U’ la traversarea unei suprafeţe oarecare Σ’
(∫∫ΣΣ
⋅⋅⋅+⋅+⋅=σ⋅⋅=''
dvduJUJUJUdnUQ ''z
''y
''x
'''321 ) , (3.138)
fiind cosinuşii directori ai normalei la şi u, v coordonate parametrice
oarecare ale suprafeţei .
'iJ 'Σ
'Σ Dacă se caută fluxul vectorului U prin suprafaţa omoloagă definită prin aceleaşi coordonate parametrice u, v, se obţin următoarele relaţii între Jacobienii (cosinuşii directori) J
Σ
1, ai celor două suprafeţe Σ şi Σ'J1’
( )( ) 1
2
1 JKK
Kv,uDz,yDJ
zy
''' ⋅
⋅== , ( )
( )v,uDz,yDJ =1 ,
înlocuind U cu U'x x în integrala (3.138) se vede că:
QKKK
KQzyx
' ⋅⋅⋅
=3
(3.139)
sau
Dispersia poluanşilor în apele subterane 87
'zyx QK
KKKQ ⋅
⋅⋅= 3 . (3.140)
Acestea sunt relaţiile dintre debitele din sistemul anizotrop şi sistemul izotrop echivalent. Pentru ca aceste debite să fie identice este suficient ca relaţia dintre conductivităţile hidraulice, în cele două sisteme, să fie : 3 zyx KKKK ⋅⋅= . (3.141) Deci trecerea de la un mediu anizotrop la unul izotrop se face calculând un coeficient de conductivitate hidraulică echivalent, de forma (3.141), făcând toate calculele hidraulice în noul sistem şi trecând apoi în vechiul sistem cu relaţia (3.140).