Vanadis – model rozprzestrzeniania się zanieczyszczeń w atmosferze
W rozwiązaniu zastosowana jest metoda Eulera, w której układ współrzędnych i siatka elementów skończonych związane są ze źródłem zanieczyszczenia.
Zanieczyszczenie przemieszcza się względem siatki z prędkością vx, vy, vz . Przyjęty sposób dyskretyzacji obszaru za pomocą elementów ośmiowęzłowych i trójkątnych zobrazowany jest na poniższym rysunku
Dyskretyzacja obszaru trójwymiarowego (a) za pomocą elementów ośmiowęzłowych oraz obszaru dwuwymiarowego (b) za pomocą elementów trójkątnych
W wyniku odpowiedniego podziału obszaru na elementy uwzględniona zostaje topografia terenu. Rozwiązanie równania dyfuzji metodą elementów skończonych w przestrzeni trójwymiarowej polega na określeniu pola stężeń w węzłach elementów.
Warunki brzegowe i parametry modelu
Równanie dyfuzji rozwiązywane jest przy przyjęciu na części powierzchni obszaru warunku brzegowego pierwszego rodzaju S(x,y,z)=0
Warunek ten zakłada całkowity zanik zanieczyszczenia w dużej odległości od źródła emisji. Natomiast przyjęcie warunku brzegowego różnego od zera (S(x,y,z)>0) równoznaczne jest z uwzględnieniem tła zanieczyszczenia.
Na powierzchni ziemi zanieczyszczenie jest pochłaniane; gęstość pochłanianego strumienia zanieczyszczeń jest proporcjonalna do współczynnika wnikania (prędkość suchego osiadania).
W całej objętości zanieczyszczenie ulega rozpadowi z czasem połowicznego zaniku ln2/P (współczynnik zaniku P może być funkcją stężenia zanieczyszczenia, funkcją położenia i czasu oraz warunków meteorologicznych, umożliwia to uwzględnienie reakcji chemicznych, opadów i wszystkich tych czynników, które wpływają na rozpad, bądź też powstawanie zanieczyszczeń). W wybranych elementach umieszczone jest źródło zanieczyszczenia o objętościowym natężeniu emisji Qv
Numeryczne rozwiązanie równania transportu masy napotyka na trudności związane z występowaniem członu konwekcji, który może powodować do rozwiązań oscylacyjnych. Otrzymanie zadowalających wyników może być osiągnięte poprzez zmniejszenie wymiarów siatki elementów lub przez zastosowanie metody reszt ważonych z niesymetrycznymi funkcjami wagi. W tym rozwiązaniu zastosowałem metodę Galerkina.
Obliczanie wartości elementów składających się na układ równań liniowych jest realizowane metodą Gaussa. Układ równań rozwiązany jest metodą gradientów sprzężonych.
Przykładowe obliczenia (dane z pliku konfiguranyjnego vanadis.cfg):
Q[kg/s] K[m^2/s] P[1/s] V1[m/s] V2[m/s] V3[m/s] ALPHA[m/s]
*********************************************************************************
0 1 2 3 4 5 6 7 8
012345678901234567890123456789012345678901234567890123456789012345678901234567890
*********************************************************************************
+001.00E+00 +040.00E+00 +07.00E-04 +00.00E+00 -00.00E+00 +01.50E+00 +01.00E-03
Q[kg/s] K[m^2/s] P[1/s] V1[m/s] V2[m/s] V3[m/s] ALPHA[m/s]
*********************************************************************************
0 1 2 3 4 5 6 7 8
012345678901234567890123456789012345678901234567890123456789012345678901234567890
*********************************************************************************
+001.00E+00 +040.00E+00 +07.00E-04 +00.00E+00 -00.00E+00 +00.00E+00 +01.00E-03
Porównanie modelu trójwymiarowego z modelem Pasquilla
W celu porównania modeli przeprowadziłem obliczenia rozprzestrzeniania się zanieczyszczeń wokół źródła emisji.
Rozwiązanie Pasquilla opracowałem w oparciu o „Wytyczne obliczania stanu zanieczyszczenia powietrza atmosferycznego Ministerstwo Ochrony Środowiska, Zasobów Naturalnych i Leśnictwa, 1993” autorstwa J. Iwanka.
W rozwiązaniu Pasquilla przyjęte są następujące dane:
· natężenie emisji 1 g/s
· położenie źródła emisji X=0, Y=0, H=60 m
· typ równowagi atmosfery 4
· średnia prędkość wiatru 1.55 m/s
· czas połowicznego zaniku zanieczyszczenia 1000 s
· szorstkość terenu 1.5 m
Rozkład zanieczyszczeń otrzymany modelem Pasquilla w płaszczyźnie XY dla Z=0 (a) oraz w płaszczyźnie YZ dla X=900m (b); vy=0
Rozkład zanieczyszczeń otrzymany metodą elementów skończonych w płaszczyźnie XY dla Z=0 (a) oraz w płaszczyźnie YZ dla X=900m (b); vy=0