DYSPERSJA ZANIECZYSZCZEŃ CHEMICZNYCH W OŚRODKU POROWATYM I JEJ SYMULACJA NA MODELACH CYFROWYCH

Krzysztof Mazurek, Andrzej Pawuła

[Biuletyn Instytutu Kształtowania Środowiska, 12, 19-25, 1980, Warszawa]

               Pojęcie dyspersji zanieczyszczeń stosowane jest do ogólnego określenia zjawiska rozprzestrzeniania się substancji  chemicznych w badanym ośrodku, a także do parametrycznego zdefiniowania jednego z czynników procesu wodnej migracji zanieczyszczeń w ośrodku porowatym, a mianowicie dyspersji hydrodynamicznej. Proces migracji. zanieczyszczeń uwarunkowany jest kompleksem takich czynników, jak: przepływ wody, dyspersja hydrodynamiczna, dyfuzja molekularna, infiltracja, sorpcja i wymiana jonowa. Aktywność poszczególnych czynników jest przy tym  uzależniona od własciwości fizykochemicznych ośrodka  i substancji zanieczyszczających. Jak wykazują bowiem badania eksperymentalne Phannkucha [11] oraz Biggara i Nielsena [3], dyfuzja molekularna odgrywa istotne znaczenie jedynie przy bardzo małych prędkościach wody w ośrodku słabo przepuszczalnym. Dla nasyconego wodą ośrodka porowatego, w warunkach filtracji zgodnych z prawem Darcy, istotne znaczenie posiada przede wszystkim dyspersja hydrodynamiczna /mechaniczna/ zwiazana z porowatościa ośrodka i zróżnicowaniem prędkości wody w przestrzeni porowatej i krętością elementarnych strug wody. Dyfuzja molekularna ujawnia  się natomiast przy bardzo małych prędkościach porowych, występujących w utworach półprzepuszczalnych, a także w tej części strefy nienasyconej, która nie bierze aktywnego udziału w przepływie wody.

              Parametryczne określenie roli pozostałych czynników fizykochemicznych jest znacznie trudniejsze z uwagi na wielkie zróżnicowanie zanieczyszczeń  chemicznych oraz dynamikę i skalę tych zjawisk w ośrodku gruntowym. Podejmowane są jednak próby szczegółowszego rozpoznania i opisu matematycznego tych zjawisk [2, 4, 13]. W opisie matematycznym dyspersji hydrodynamicznej korzysta się z analogii zjawiska dyfuzji. Wyprowadzając równanie dyspersji należy oprzeć się na II prawie Ficka dla dyfuzji, w postaci:

 dC/dt = D^m div grad C                                                 /1/

gdzie: dC/dt - pochodna cząstkowa koncentracji molowej C, względem czasu t;  D^m - macierz współczynnika dyfuzji,

 Przy formułowaniu równania /1/ przyjęto następujące założenia:

 -  brak reakcji chemicznych,

 -  zaniedbywalnie mała szybkość konwekcji,

 -  izotermiczność i izobaryczność,

 -  brak zależności współczynnika od stężenia.

Ze względu na to, że w wielu przypadkach przepływ odgrywa istotną rolę dla zjawiska dyspersji, autorzy przyjęli następujące założenia dodatkowe: na całym obszarze przepływ odbywa się zgodnie z prawem Darcy'ego, a więc jest laminarny i liniowy. Przy tych założeniach zjawisko dyspersji hydrodynamicznej można opisać poniższym układem równań różniczkowych cząstkowych:

 div k grad H = S dH/dt + Q

  v = - k grad H                                                             

 div D grad C - div (v/n × C) = dC/dt × R + Q'              /2/

 gdzie:

dH/dt - pochodna cząstkowa potecjału hydraulicznego H, względem czasu t

dC/dt - pochodna cząstkowa koncentracji molowej C, względem czasu t

k  -  tensor współczynnika filtracji,

D  -  tensor współczynnika dyspersji,

S  -  współczynnik pojemności wodnej,

R - współczynnik opóznienia wskutek adsorpcji,

v -  prędkość filtracji (prędkość fikcyjna wg. hipotezy  Darcy),

n - porowatość efektywna,

Q - wydatek lub pobór wody

Q' - natężenie przepływu masy zanieczyszczenia na jednostkę objętości.

 Tensor współczynnika dyspersji /jak i tensor współczynnika filtracji/ po przyjęciu uproszczeń można wyrazić w postaci macierzy kwadratowej drugiego stopnia, zawierającej tylko współczynniki zgodne co do kierunku z przyjętymi. osiami układu współrzędnych. Wtedy układ równań /2/ można przedstawić w postaci równania różniczkowego cząstkowego:

D_xxd²C/dx² + D_yy d²C/dy² + C (dH/dt + Q) + (v_x/ n)grad C + (v_y/n)grad C = R dC/dt + Q'      /3/

gdzie:

d²C/dx² - pochodna cząstkowa drugiego rzędu koncentracji masy C, względem osi układu x, y

dH/dt - pochodna cząstkowa potecjału hydraulicznego H, względem czasu t

dC/dt - pochodna cząstkowa koncentracji masy C, względem czasu t

H = H(x, y, t),

C = C (x, y, t),

D_xx, D_yy  - współczynniki dyspersji podłużnej i poprzecznej,

v_x, v_y      - składowe wektora prędkości filtracji Darcy.

 Zakłada się przy tym, że ośrodek jest anizotropowy i jednorodny, a współczynnik dyspersji podłużnej jest zawsze zgodny z kierunkiem największego gradientu hydraulicznego w danym punkcie.

 Dla równania /3/ określamy warunki poczatkowe i brzegowe:

Warunki poczatkowe

H(x, y, 0) = Ho(x, y);               (x, y) Î W u G                /4/

C(x, y, 0) = Ho(x, y);                (x, y) Î W u G

 gdzie:

Ho i C - funkcje dane,

W -  wnętrze rozpatrywanego obszaru,

G - brzeg rozpatrywanego obszaru;

 Warunki brzegowe

-   I rodzaju   H(x, y, t)/G = H1(x, y, t);       0 =< t = < T                            

                    C(x, y, t)/G = C1 (x, y, t);       0 =< t = < T         /5/

 gdzie:              H1 i C1  funkcje dane;

- II rodzaju   Q(x, y, t)/G = Qo (x, y, t);       0 =< t =< T                               

                        Q'(x, y, t)/G = Q'o (x, y, t);  0 =< t =< T         /6/

 gdzie: Qo i Q'o funkcje dane.

 Warunki I rodzaju określają stałe wartości potencjału hydraulicznego oraz koncentracji (moga być one również określone dla punktów położonych wewnątrz obszaru W). Natomiast warunki II rodzaju określają  stały wydatek przepływu i koncentracji  w kierunku normalnej do brzegu G.

              Równanie /3/ z warunkami granicznymi /4/, /5/ i /6/ rozwiązano stosując metodę różnic skończonych. W tym celu cały rozpatrywany obszar W wraz z brzegiem G pokryto kwadratową siatką o boku A. Dla każdego z bloków o powierzchni A²  określa się parametry wchodzące  w skład równania /3/, przyjmując za punkt obliczeniowy środek bloku. W tak zdyskretyzowanym obszarze, pochodne występujące w równaniu /3/ przybliżono różnicami skończonymi, otrzymując  dla każdego bloku obliczeniowego /I/ następujący schemat różnicowy:                                             

                                                    N

                                          W       I       E

                                                    S

R C_I ^t2 = R C_i ^t1  + Dt { C_W^t2 [D_xx/A² + 1/n_W (IV_WI-V_W)/2A] + C_E^t2 [D_xx/A² + 1/n_E (IV_EI-V_E)/2A] + C_N^t2 [D_yy/A² + 1/n_N (IV_NI-V_N)/2A] + C_S^t2 [D_yy/A² + 1/n_S (IV_SI-V_S)/2A] +    C_i ^t2 [-2 D_xx/A² - 2 D_yy/A² - 1/n_N (IV_NI+V_N)/2A) - 1/n_S (IV_SI+V_S)/2A) - 1/n_W (IV_WI+V_W)/2A) - 1/n_E (IV_EI+V_E)/2A] - C_I'}

gdzie:

-         V_N, V_S, V_W, V_E  - prędkości filtracji z oczka centralnego I do oczek oznaczanych N, S, E, W,

-        IV_NI, IV_SI, IV_WI, IV_EI - wartości bezwględne prędkości filtracji

-          n_N, n_S, n_W, n_E - średnie wartości porowatości efektywnej pomiędzy oczkami centralnymi i odpowiednio oczkami N, S, W, E siatki modelu o boku A.

Pochodną czasową występującą po prawej stronie równania /7/ można aproksymować jednym ze schematów:

-         schemat prosty (ekspilicité),

-         schemat uwikłany (implicité),

-         schemat mieszany (metoda Cranka - Nicholsona) [1].

Ze względu na swoją stabilność dla dowolnego kroku czasowego i bezwarunkową zbieżność, wybrano metodę  implicité (metoda uwikłana).

Po utworzeniu dla każdego bloku obliczeniowego schematu różnicowego /7/, otrzymuje się układ równań liniowych (ilosć równań jest równa ilości bloków modelu). Do rozwiązania tego układu przyjęto metodę iteracyjną Gaussa - Seidla z nadrelaksacją.

              Opisany powyżej sposób postępowania został zastosowany przez autorów przy tworzeniu programu symulacyjnego "DYSKON". Zadaniem tego programu jest symulacja przepływu wód podziemnych w jednej warstwie wodonośnej, z jednoczesnym uwzględnieniem dyspersji zanieczyszczeń chemicznych.

              Możliwości programu są szerokie i pozwalają na symulację:

-         przepływu ustalonego,

-         przepływu nieustalonego ze stałym zmiennym krokiem czasowym,

-         dyspersji zanieczyszczeń chemicznych przy założeniu przepływu ustalonego zarówno przy zmiennym, jak i stałym kroku czasowym,

-         przepływu i dyspersji zanieczyszczeń przy przepływie nieustalonym zarówno przy zmiennym, jak i stałym kroku czasowym.

W trakcie przygotowywania danych do modelowania wymagane jest ustalenie warunków  poczatkowych i brzegowych:

    -   poczatkowego potencjału hydraulicznego oraz początkowych wartości koncentracji dla każdego bloku modelu,

-         stałych wartości potencjału hydraulicznego na brzegu modelu (wa­runek I rodzaju) lub stałych wartości przepływu dla bloków brzegowych (warunek II rodzaju),

-         warunków brzegowych j.w. dla bloków położonych wewnątrz modelu,

-         wartości wydatku lub zasilania dla bloków wewnętrznych modelu (wydatek może być zmienny w czasie),

-         stałych wartości   koncentracji dla bloków, w których znajdują się ogniska emitujące stały ładunek zanieczyszczeń (warunek I rodzaju) lub wielkości ładunku przepływającego w określonym czasie przez blok brzegowy (warunek II rodzaju).

Ponadto wymagane jest ustalenie wielkości parametrów hydraulicznych występujących w równaniu dla każdego bloku modelu:

-         współczynników filtracji dla kierunków pokrywających się z osiami współrzędnych (blokami siatki),

-         współczynników dyspersji poprzecznej i podłużnej,

-         współczynników pojemności wodnej,

-         współczynników opóźnienia wskutek adsorpcji.

              W wyniku działania programu otrzymujemy rozkład potencjałów hydraulicznych i koncentracji zanieczyszczeń dla każdego bloku w wybranym momencie czasu symulacji (dla kroków czasowych). Wyniki otrzymane z symulacji mogą posłużyć do ustalania przybliżonej prognozy przepływu zanieczyszczeń, jak również do wstępnego oszacowania zasięgu oddziaływania ognisk zanieczyszczeń na wody podziemne.

              Ogólne tendencje światowe wskazują na celowość kontynuowania prac rozwijających modelowe badania migracji zanieczyszczeń. Wiaże się to z łatwością modelowania nawet dużych i geometrycznie skomplikowanych obszarów oraz dużą łatwością wprowadzania i korekcji danych wejsciowych. Minusem stosowania metod modelowania cyfrowego dla dużych obszarów może być wysoki koszt i czasochłonność pracy maszyny. Można jednak starać się skrócić czas obliczeń przez wybór optymalnego kroku siatki, kroku czasowego oraz optymalnego doboru procesu iteracyjnego.

              Wybrany przez autorów sposób rozwiązania problemu, przy zastosowaniu metody uwikłanej, zapewnia dużą stabilność rozwiązania, natomiast zastosowanie metody iteracyjnej Gaussa - Seidla z nadrelaksacją gwarantuje szybką zbieżność do rozwiązania dokładnego i mały błąd metody (przy prawidłowym wyborze przybliżenia początkowego). Dla usprawnienia procesu obliczeń, współczynnik nadrelaksacji wyznaczany jest dynamicznie w programie. Otrzymane wyniki mogą więc być obarczane jedynie błędami wynikającymi z wadliwego rozpoznania hydrogeologicznego i dyskretyzacji, jak również z przyjęcia błędnej wielkości występujących współczynników.

              Autorzy wskazują również na celowość kontynuowania prac w kierunku modelowania dyspersji w powiązaniu z procesami przemian zanieczyszczeń w ośrodku gruntowym. Prace te powinny być jednak powiązane z pracami. teoretycznymi i eksperymentalnymi, mającymi na celu wyznaczenie parametrów zjawisk towarzyszących dyspersji hydrodynamicznej. Pozwoliłoby to na uszczegółowienie modelu matematycznego zjawiska migracji zanieczyszczeń.

 Bibliografia

      [1] ANGEL E., BELLMAN R.: Dynamic programming and partial differential equations. New York--London: Academic Press 1972

[2] BABU D.K.: Infiltration analysis and perturbation methods. Adsorption with exponential diffusivity. Water Resources Research, 1976 T.12 nr 1 s. 89--93

[3] BIGGAR J., NIELSEN D.: Diffusion effects in miscible displacement occurring in saturated and unsaturated porous materials. Journal Geophysical Research of USA, 1964 Z. 65 nr 9

[4] CEARLOCK D.E.: A system approach to management of the ground water basin. The National Ground Water Quality Symposium in Denver Col., 1971/08/

[5] CHRUSZCZYŃSKI A., KAŹMIERCZAK-WIJURA Z., MAZUREK K., MOSSOR A., PAWUŁA A., SCHMIDT E.: Metodyka użytkowania programu SIMONE w systemie cyfrowym Qdra 1300. 1976, 43 s. maszyn. IKŚ. Zakład Ujęć i Ochrony Wód Podziemnych, Poznań

[6] COLLATZ L.: Metody numeryczne rozwiązywania równań różniczkowych. Warszawa: PWN 1960

[7] EMSELLEM Y.: Construction de modeles mathématiques en hydrogéologie. CIG Fontainebleau 1971

[8] FRIED J.J.: Etudes théorétiques et méthodologiques de la dispersion en milieu naturel. These de Doctorat d'Etat es Sciences Physiques. Bordeaux 1972

[9] LESSI J.: Simulation numérique de la propagation d'un polluant dans un milieu poreux saturé. These de Docteur - ingénieur. Université Louis Pasteur, Strasbourg 1976

[10] PAWUŁA A.: Sprawozdanie ze stażu naukowego we Francji w zakresie badań nad migracją zanieczyszczeń chemicznych do wód pod­ziemnych oraz konstrukcji modeli matematycznych w hydrogeologii (cz. II). 1976, 19 s. maszyn. IKŚ, Poznań

[11] PHANNKUCII H.0.: Contribution a l'étude des déplacements de fluides miscibles dans un milieu poreux. Revue de l'Institut Français du Pétrole, 1963, no 2--18

[12] RALSTON A.: Wstęp do analizy numerycznej. Warszawa: PWN 1971

[13] TARDY Y.: Equilibres solutions - minéraux dans les aquiferes superficiels modeles de simulation thermodynamique. Centre Nationale des Recherches Scientifiques 1975.

__________________________________________________________________________

tekst artykulu w edytorze WORD: Biuletyn.doc

adres do korespondencji: pawula@main.amu.edu.pl