Hidrológiai Közlöny 1986 (66. évfolyam)
4-5. szám - Székely Ferenc: Felszín alatti vizek konvektív kémiai tömegtranszportjának numerikus modellezése
256 HIDROLÓGIAI KÖZLÖNY 1986. 66. [ÉVFOLYAM 3. SZAM A vázolt fizikai modell alapján meghatározható a koncentrációmező számítására alkalmas matematikai modell, amelynél az alábbi egyszerűsítő feltételezéseket tesszük: — a transzportfolyamat izotermikus ill. nem függ a hőmérsékletváltozástól; — a szivárgási folyamat permanens; — a fedőréteg aktív vagy dinamikai pórustérfogatát hozzáadjuk a vízvezető réteg pórustérfogatához és a hidraulikai szintet egy összevont porozitással jellemezzük. Ekkor a függőleges átszivárgási időt elhanyagoljuk, a laterális szivárgási időt pedig kissé megnöveljük. A fentiek alapján az i-ik hidraulikai szint C t (*, y, t) [g/m 3] koncentrációeloszlását az alábbi parciális differenciálegyenlet rendszer írja le: [dCi(x, y, t)ldt+kCi(x, y, *)]X XMi(x, y)Ni(x, y)Ai(x, y) = = div [G'i(x, y, t)T t(x, y) grad H,(x, y)] -Ct(x, y, t)qi(x, y) + +Cf(:r, y, t)pi(x, y) + + Bi(x, y)C*^{x, y, t)AHi_ x(x, y) + + B i+ 1(x, y)C* + 1(x, y, t)AHi + l(x, y) (• = 1,2, ...,») (1) AH k{x, y) =H k(x, y) — Hi(x, y) (2) C k(x, y,t)=[ Ci(x>yt) haAm< 0 (3) ahol Hi, Mj, Ni, T u Bi — az i-edik hidraulikai szint permanens piezometrikus nyomása [m], vastagsága [m], porozitása, transzmisszibilitása [m 2/nap], a fedőréteg függőleges átszivárgási tényezője [l/nap]; pi, cji — a külső utánpótlódás vagy megcsapolás területegységre vonatkozó fajlagos hozama [m/nap]; C Ti — a p, beszivárgás koncentrációja [g/m 3]; X = (ln2)IT 1i 2 — a radioaktív bomlási állandó [l/nap], T i n felezési idő [nap]; — az egyensúlyi adszorpció dimenzió nélküli együtthatója (egy adott kőzettestben tárolt összes, valamint a vízfázisban oldott kémiai tömegek aránya); Gl—- a grad II vektorhoz csatlakozó nyomvonal (karakterisztika) koncentrációja [g/m 3]. A II 0 (x, y) és G 0 (x, y, t) függvények a felső peremfeltételt képező folyók, tározók vagy esetleg egy feljebb települő vízvezető réteg piezometrikus nyomását és koncentrációját adják meg. Az (1) — (3) egyenletrendszer megoldásához szükség van még az egyes vízvezető szintek elterjedési területének (geometriájának), valamint kezdeti koncentrációjának ismeretére. Megjegyezzük, hogy az (1)—(2) egyenletek alapján a permanens piezometrikus nyomáseloszlás is számítható azzal a módosítással, hogy az (1) egyenlet bal oldalát zérusnak vesszük, továbbá a jobb oldalon szereplő Ci, Gl, Ci és C" függvényeket elhagyjuk (1-gyel helyettesítjük). 3. A COMAÜ numerikus transzportmodell A COMAD-eljárás kidolgozásakor azt a célt tűztük ki, hogy egyrészt kiküszöböljük az ún. numerikus diszperzió batását és így lehetővé váljon az éles koncentrációfrontok konvektív mozgásának előrejelzése, másrészt pedig lehetőséget biztosítsunk a numerikus módszerek alkalmazásakor elkerülhetetlen tömegmérleg-hiba hatékony korrekciójához. Az (1) egyenlet szabatos numerikus megoldása a hagyományos differencia- vagy végeselem módszerekkel csak bizonyos feltételek mellett lehetséges. Egydimenziós szivárgás és függőlegesen zárt egyszintes tároló (B=p=q=0) esetében csak abban a triviális esetben kapunk pontos megoldást, ha adott At időlépcsők mellett az összes Ax szakaszra teljesül a Ax MNA=AtT\ bH\bx\ Courant-feltétel. Ekkor ugyanis Zlt időlépcsőként a koncentrációértékek a konvektív áramlás irányában egy csomóponttal eltolódnak, amit az explicit differenciaséma pontosan leír. Nyilvánvaló, hogy két- ill. háromdimenziós feladatoknál a lokális szivárgási sebesség értéke 2—3 nagyságrenden belül is változik, ezért a Courant-feltétel az összes csomópontra egyidejűleg nem teljesül. Ennek következtében a tároló kis szivárgási sebességgel jellemezhető részein a koncentráció a valóságosnál gyorsabban terjed, az éles frontok elmosódnak és kialakul az ún. longitudinális (hosszirányú) numerikus diszperzió (LND). Hasonló tulajdonságokkal rendelkeznek az ún. implicit differenciasémák is. Két- és háromdimenziós feladatok rácsmódszerrel történő megoldásakor valamely csomópont lokális áramvonala (karakterisztikája) általában nem megy át a szomszédos rácspontokon. Ilyen esetekben a beáramló víz koncentrációját keresztirányú interpolációval számítják, ami a koncentrációmezőnek az áramlás irányára merőleges kiegyenlítődéséhez, az éles határok elmosódásához vezet. Ezt a jelenséget nevezik keresztirányú vagy transzverzális numerikus diszperziónak (TDN). Konvektív transzportfeladatok numerikus megoldásakor hiba terheli nem csak a koncentráció eloszlását, hanem a rendszer egész tömegmérlegét is, mivel az előrejelzett és kezdeti tárolt tömegek különbsége eltér a számítási időszak alatti külső tömegáram időbeni integráljának és a fiziko-kémiai folyamatok hatására bekövetkező koncentrációváltozásoknak az összegétől. A vázolt nehézségek áthidalására kifejlesztett COMAD-eljárás lényegét tekintve differenciamód1. ábra. A COM AD transzporimodell A—D hálóelemei, 0—16 csomópontjai és 0 differenciaeleme (vonalkázott terület)