46
Komputerowa symulacja
rozprzestrzeniania
zanieczyszczeń w atmosferze
Jacek Piekarski, Marcin Lubierski
Katedra Techniki Hydro-Szlamowej
i Utylizacji Odpadów
Politechnika Koszalińska
1. Wstęp
Atmosfera otaczająca powierzchnię Ziemi dzięki swojej naturze posia-
da zdolność zmniejszania i kumulowania stę\
enia oraz przemieszczania emito-
wanych zanieczyszczeń, które następnie ulegają licznym przemianom fizyko-
chemicznym. Bezpośredni wpływ na te przemiany ma zarówno cyrkulacja
obejmująca kontynent, a nawet całą planetę, jak i nakładające się nań ró\
no-
rodne systemy prądów powietrznych mniejszej skali. Znajomość procesów
zachodzących w atmosferze, warunków meteorologicznych panujących na da-
nym obszarze jest niezbędna dla prawidłowego przewidywania skutków emisji
zarówno ze z
ródeł istniejących jak i projektowanych [6,10].
Na potrzeby prognozowania rozprzestrzeniania stÄ™\
eń zanieczyszczeń
opracowano szereg modeli matematycznych analityczno-empirycznych. Modele
wykorzystywane są przy tworzeniu systemów alarmowych oraz do specyfikacji
obszarów zagro\
eń na wypadek wystąpienia du\
ych emisji losowych. Mo\
na je
stosować do symulacji występujących stę\
eń w atmosferze, co jest wskazane ze
względu na wysokie koszty odpowiednich pomiarów. Dlatego te\
znajÄ…c para-
Jacek Piekarski, Marcin Lubierski
metry emitora, warunki meteorologiczne i terenowe, z wystarczającą dokładno-
ścią mo\
na określić stę\
enia zanieczyszczeń w emisji.
Mo\
na wyró\
nić trzy podstawowe kategorie modeli: Eulera, Lagrange'a
oraz Gaussa. W modelu Eulera rozpraszanie się zanieczyszczeń jest opisywane
względem nieruchomego układu związanego z ziemią. Zmienne w przestrzeni
i czasie, pole prędkości przedstawione jest w ka\
dym punkcie badanego obsza-
ru, w ustalonym na powierzchni ziemi układzie współrzędnych kartezjańskich
lub sferycznych. Eulerowskie pole prędkości obrazują izogony, izotachy albo
linie prÄ…du. Gdy jest ono stacjonarne, to linie prÄ…du pokrywajÄ… siÄ™ z trajekto-
riami cząstek. W takim polu następstwo poruszających się cząsteczek powietrza
jest widziane przez nieruchomego obserwatora jako przenoszone przez wiatr.
W modelu Lagrange'a zmiany stÄ™\
enia zanieczyszczenia sÄ… opisywane
w układzie związanym z przemieszczającymi się masami powietrza. Obserwa-
tor poruszajÄ…cy siÄ™ wraz z powietrzem utrzymuje kontakt z tymi samymi czÄ…st-
kami powietrza w zadanym okresie czasu. W modelu gaussowskim stÄ™\
enie
zanieczyszczenia wzdłu\
osi smugi ma statystyczny rozkład Gaussa. Formuła ta
mo\
e być wyprowadzona przy zastosowaniu metody Eulera lub Langrage`a.
Jednak modele gaussowskie w swej najbardziej uproszczonej postaci winne być
postrzegane jako specjalne uproszczenie metody Langrange`a.
Obok trzech podstawowych modeli (Eulera, Lagrange`a i Gaussa), ist-
nieją równie\
inne określające rozprzestrzenianie się zanieczyszczeń w atmos-
ferze [2,3].
Proste modele statystyczne stosowane głównie w praktyce in\
ynierskiej
wykorzystujÄ… modele regresyjne, poczÄ…wszy od prostych tablic rozdzielczych
i regresji wielokrotnej do analizy ciągów chronologicznych. Po zainstalowaniu
sieci pomiarowej monitoringu, pojawiły się w Polsce proste modele progno-
styczne dyspersji zanieczyszczeń, stosowane w skali lokalnej.
Opracowany przez polskich matematyków model stochastyczny progno-
zuje pole stÄ™\
eń tlenków azotu, ulegających w atmosferze zło\
onym przemianom.
Poniewa\
charakterystyczną cechą tlenków azotu jest losowe powstawanie oraz
zanikanie, dlatego trajektorie cząstek NOx są nieciągłe. Niniejsze zjawiska opisy-
wane są w dwóch przestrzeniach probabilistycznych. Jedna dotyczy procesów
mikroskopowych, natomiast druga procesów makroskopowych.
Lagrangeowskie modele rozprzestrzeniania czÄ…stek sÄ… obliczeniowo
kosztowne, dlatego w przypadku większych obszarów, wykorzystuje się modele
hybrydowe lagrangeowsko eulerowskie (modele dyfuzyjno-stochasty-czne).
W otoczeniu punktowego z
ródła emisji obliczenia przeprowadzane są modelem
LPD (lagrangeowski model cząstek). Po dostatecznie długim czasie przebiegu,
gdy smuga ma du\
e rozmiary w porównaniu z przestrzennym krokiem całko-
666 VI Ogólnopolska Konferencja Naukowa
Komputerowa symulacja rozprzestrzeniania zanieczyszczeń w atmosferze
wania, cząstki włączają się do objętościowego pola emisji modelu eulerowskie-
go i zanikajÄ….
Na podstawie RozporzÄ…dzenia Ministra Åšrodowiska z dnia 5 grudnia
2002 r., do określania stanu zanieczyszczeń powietrza dla z
ródeł istniejących
stosuje się model gaussowski (formułę Pasquilla)[8].
Pasquill przystosował uproszczone równanie dyfuzji zanieczyszczenia
gazowego w poruszającym się ośrodku gazowym, do zagadnień rozprzestrze-
niania się zanieczyszczeń wyrzucanych do atmosfery z kominów. Równanie
dyfuzji w postaci Pasquilla ma postać [4,6,7,9]:
ëÅ‚ ëÅ‚
"S d Sy2 öÅ‚ "2S d Sz2 öÅ‚ "2S
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚
= Å"ìÅ‚ ÷Å‚Å" + Å"ìÅ‚ ÷Å‚ Å"
"t dt 2 "y2 dt 2 "z2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
(1)
gdzie:
S  stÄ™\
enie zanieczyszczenia w punkcie recepcyjnym,
t  czas wędrówki zanieczyszczenia od chwili wyrzutu do osiągnięcia
punktu recepcyjnego,
Sy  współczynnik dyfuzji Pasquilla wzdłu\
osi 0Y w układzie współ-
rzędnych t, y, z,
Sz  współczynnik dyfuzji Pasquilla wzdłu\
osi 0Z w układzie współ-
rzędnych t, y, z.
Równanie opisuje rozprzestrzenianie się zanieczyszczeń dla dowolnego
punktu emisji, ustalonego w czasie. Mo\
e to być z
ródło punktowe, liniowe lub
powierzchniowe  zale\
y to od przyjętych warunków brzegowych:
z
ródło emisji o masowym natę\
eniu E jest umieszczone w punkcie (z
ródło
punktowe) o współrzędnych t = 0, y = 0 i z = H,
na powierzchni podło\
a nie zachodzi pochłanianie zanieczyszczeń, tj. dla
z = 0 strumień zanieczyszczeń jest zerowy (ilość zanieczyszczeń docierają-
cych do podło\
a jest równa ilości zanieczyszczenia od podło\
a),
stÄ™\
enie zanieczyszczeń w odległości nieskończenie du\
ej od z
ródła (tzn. t,
y oraz z dÄ…\
ą do nieskończoności) ma wartość zerową.
VI Ogólnopolska Konferencja Naukowa 667
Jacek Piekarski, Marcin Lubierski
Równanie spełniające warunki brzegowe ma postać:
-y2
(z-Hp)2 -(z+Hp)2
-
2Sy2 2Sz2
2Sz2
e Å" (e + e )
S = Eg Å"
0,0063Å" Usr Å"Sy Å"Sz
(2)
gdzie:
S  stÄ™\
enie substancji gazowej w punkcie (x, y, z), [µg/m3],
Eg  maksymalna emisja substancji gazowej z z
ródła, [mg/s],
z  wysokość, dla której oblicza się stę\
enie substancji zanieczyszczajÄ…-
cej w powietrzu, [m],
Hp  efektywna wysokość emitora, [m],
Usr  średnia prędkość wiatru w warstwie powietrza od z = h do z = Hp,
[m/s],
h  geometryczna wysokość emitora, [m],
Sy  współczynnik poziomej dyfuzji atmosferycznej, [m],
Sz  współczynnik pionowej dyfuzji atmosferycznej, [m].
WykorzystujÄ…c zjawisko, \
e główny ruch zanieczyszczenia z prędkością
wiatru Usr odbywa się wzdłu\
osi x i zawsze występuje określona relacja po-
między tymi wielkościami, w wzorze (2) w miejsce zmiennej t podstawia się
współrzędną x.
Smuga gazów wylotowych wydobywa się z punktu pozornej emisji A
znajdującego się na wysokości Hp [m]  rysunek 1. Efektywna wysokość emi-
tora Hp stanowi sumę geometrycznej wysokości komina h [m] i tzw. wysokości
wyniesienia Dh [m]. Wartość wyniesienia Dh zale\
na jest od prędkości wylo-
towej gazów v [m/s], emisji ciepła Q [kJ/s] i prędkości wiatru Uh [m/s] na wy-
sokości wylotu z emitera. W zale\
ności od ilości emitowanego ciepła Q wyró\
-
niamy dwa rodzaje wyniesienia: termiczne oraz dynamiczne. Wartość wynie-
sienia termicznego obliczana jest na podstawie wzoru CONCAWE, który został
opracowany dla kominów energetycznych. Zgodnie z Rozporządzeniem Mini-
stra Åšrodowiska z dnia 5 grudnia 2002 r. [8], granicÄ™ jego zastosowania okre-
ślono dla wartości emisji ciepła powy\
ej 24000 [kJ/s]. Formuła Hollanda (obli-
cza wysokość wyniesienia dynamicznego gazów) ma zastosowanie, gdy ilość
emisji ciepła jest mniejsza od 16000 [kJ/s]. Uwzględnia ona ró\
nicę pomiędzy
prędkością wiatru na wysokości emitora Uh a prędkością gazów odlotowych v.
W przypadku, gdy wielkość emisji ciepła mieści się w przedziale od 16000 [kJ/s]
do 24000 [kJ/s] stosowane są obie formuły, pomniejszone o ró\
nicę pomiędzy
668 VI Ogólnopolska Konferencja Naukowa
Komputerowa symulacja rozprzestrzeniania zanieczyszczeń w atmosferze
wartością emisji ciepła Q a wartościami granicznymi. Powy\
sze wzory na wynie-
sienie gazów wylotowych zakładają pionowy wyrzut gazów do atmosfery. Stąd
przy wyrzucie poziomym, w wyniku zadaszenia czy wyrzutni poziomej nale\
y
przyjąć zerową wartość wyniesienia.
Rys. 1. Graficzna interpretacja formuły Pasquilla (opis w tekście) [6]
Rys. 1. The graphic interpretation of Pasquill formula (description in text) [6]
Wymiar poprzeczny smugi zanieczyszczeń stanowią półosie elipsy
o powierzchni, przez którą przepływa 90% masy wydalanej z punktowego z
ró-
dła emisji, w tym tak\
e 90% określonego zanieczyszczenia  rysunek 1. Przy
takim zało\
eniu półosie elipsy wynoszą 2,15 Sy i 2,15 Sz. Współczynniki dyfu-
zji Sy i Sz mają wymiar długości. Są one proporcjonalne do wymiarów po-
przecznych rozszerzającej się i niezakłóconej powierzchnią ziemi smugi. Wraz
ze wzrostem odległości od z
ródła emisji wartości współczynników dyfuzyjnych
rosną ze względu na rozpraszanie się smugi. Istnieje kilka metod obliczania ich
wartości, najczęściej stosowaną w Polsce jest metoda Nowickiego. Charakter
zmian współczynników opisany jest przez wyra\
enia [1,4,8]:
Sy = A Å" xa
(3)
VI Ogólnopolska Konferencja Naukowa 669
Jacek Piekarski, Marcin Lubierski
b
Sz = B Å" x
(4)
gdzie:
A, B,  współczynniki empiryczne, [m],
a, b  współczynniki empiryczne, [ ].
Wykładniki potęgowe a i b charakteryzują dynamikę rozpraszania smu-
gi zanieczyszczeń w atmosferze. Maleją wraz ze zmniejszaniem się intensyw-
ności turbulencji  dotyczy to szczególnie wykładnika b decydującego o kącie
pionowego rozwarcia smugi. Zale\
ne są od stanu równowagi atmosfery, przy
czym nie są one opisane w postaci reguły matematycznej, lecz ich wartość do-
biera siÄ™ na podstawie tablic [4,6].
Parametry A i B zale\
ne są od warunków terenowych, w których pro-
wadzony był pomiar współczynników dyfuzji atmosfery. Stosunek pomiędzy
efektywną wysokością emitera a aerodynamiczną szorstkością terenu jest kryte-
rium opisującym wpływ czynników terenowych na wielkość dyfuzji. Uwzględ-
niony jest równie\
wpływ stanu równowagi atmosfery poprzez wykładnik me-
teorologiczny m dobierany podobnie jak wykładniki a i b [4]:
Hp
A = 0,088Å" (6m-0,3 +1- ln )
Z0
(5)
Hp
B = 0,38m1,3 Å" (8,7 - ln )
Z0
(6)
gdzie:
m  współczynnik empiryczny, zale\
ny od stanu atmosfery, [ ],
Hp  efektywna wysokość emitora, [m],
Z0  aerodynamiczny współczynnik szorstkości terenu, [m].
W modelu uwzględnione są (rysunek 2): prędkość wiatru na wysokości
emitora Uh [m/s], średnia prędkość wiatru w warstwie powietrza od z = h do
z = Hp Usr [m/s] oraz prędkość wiatru na wysokości anemometru Ua [m/s].
Prędkość wiatru Uh, obliczana jest w zale\
ności od wysokości emitora h, pręd-
kości wiatru na wysokości anemometru Ua oraz stanu atmosfery (parametru
meteorologicznego m). Natomiast wartość średniej prędkości Usr zale\
y od
efektywnej wysokości emitora Hp, wysokości emitora h, prędkości wiatru na
wysokości anemometru Ua oraz stanu atmosfery (parametr m).
670 VI Ogólnopolska Konferencja Naukowa
Komputerowa symulacja rozprzestrzeniania zanieczyszczeń w atmosferze
Z
A
Usr
U
Uh
anemometr
Ua
X
Rys. 2. Rozkład prędkości wiatru wraz z wysokością w modelu Pasquilla [6]
Rys. 2. The distribution of wind velocity of along with the height in Pasquill model [6]
Aerodynamiczny współczynnik szorstkości terenu Z0 ma wymiar dłu-
gości i uwzględnia wpływ pokrycia powierzchni na intensywność rozpraszania
zanieczyszczeń w atmosferze. Elementy o zró\
nicowanej wysokości mają bez-
pośredni wpływ na poziomy ruch mas powietrza, rozdzielając je na poszcze-
gólne strumienie i powodując zawirowania turbulencyjne.
Algorytm umo\
liwiajÄ…cy obliczenie stÄ™\
enia zanieczyszczeń S, współczyn-
ników dyfuzji Sy i Sz, wysokości Hp pozornego punktu emisji A oraz pozostałych
niezbędnych parametrów emitora i meteorologicznych przedstawia rys. 3.
2. Opis programu
Na podstawie algorytmu przedstawionego na rysunku 3 wykonano pro-
gram komputerowy  Atmo słu\
Ä…cy do symulacji rozprzestrzeniania siÄ™ zanie-
czyszczeń w atmosferze. Program pracuje w środowisku Windows
95/98/XP/NT. Składa się z dwóch podstawowych modułów projektowych
umo\
liwiających wprowadzenie danych wstępnych i wyprowadzenie wyników
obliczeń oraz dwóch modułów pomocniczych słu\
Ä…cych wizualizacji przepro-
wadzonej symulacji poprzez przedstawienie zmian wartości współczynników
dyspersji oraz stÄ™\
enia w zale\
ności od zmiany wartości odległości receptora
od emitera. Aplikacja posiada opcjÄ™ zapisu do pliku oraz odczytu z pliku
wprowadzonych danych poczÄ…tkowych. Pliki skojarzone z aplikacjÄ… majÄ… roz-
szerzenie  *.zan . Aplikacja generuje raport końcowy oraz daje mo\
liwość
wydruku wprowadzonych danych początkowych i wyników obliczeń.
VI Ogólnopolska Konferencja Naukowa 671
Dh
Hp
h
14
Jacek Piekarski, Marcin Lubierski
S
h, d, v, T, Z0, Ua, T0, ATM, x, y, z, Eg
T N
278,76·d2·v·(T - T0)
h > 300 Uh = 21,429·Uam Uh < 0,5 Q =
T
T
N
Uh = Ua·(0,071·h)m Uh = 0,5
N
DhH = 0 v > 0,5Uh
Hollanda
T
N
1,5 Å" v + 0,00974Å" Q N
Dh=DhH Q > 16000
DhH = v < Uh
Uh
T
CONCAWE
T
1,5Å" v + 0,00974Å"Q v - 0,5Uh
1,126 Å" Q0,58 N
DhH = -
DhC = Q > 24000
Uh 0,5Uh
Uh0,7
T
Hollanda + CONCAWE
T
24000- Q 1,126Å"Q0,58 Q -16000
Dh = DhH Å" + Å" Usr < 0,5 Usr = 0,5
8000 Uh0,7 8000
N
Hp=h+Dh
Hp T Hp
< 10 =10
Z0 Z0
T
N
Hp = h Usr = Ua·(0,071·h)m
Hp T Hp
N
> 1500 =1500
Z0
Z0
N
Ua Å"(Hp1+m - h1+m)
Hp > 300 Usr =
(Hp - h)(m +1)Å"14m N
T
Hp öÅ‚
N A = 0,088Å"ëÅ‚6m-0,3 +1- ln
ìÅ‚ ÷Å‚
h < 300 Usr = 21,429·Uam
Z0
íÅ‚ Å‚Å‚
T
Hp öÅ‚
îÅ‚3001+m - h1+ m
Å‚Å‚
Ua
ìÅ‚
Usr = Å" + (Hp - 300)Å"300m śł B = 0,38Å" m1,3 Å"ëÅ‚8,7 - ln Z0 ÷Å‚
íÅ‚ Å‚Å‚
(Hp - h)Å"14m ïÅ‚ 1+ m
ðÅ‚ ûÅ‚
2 2
Å‚Å‚
öÅ‚ ëÅ‚ öÅ‚ Sy = A·xa
y2 îÅ‚ ëÅ‚ (z - Hp) (z + Hp)
÷Å‚ ÷łśł
exp(- ) Å" expìÅ‚- + expìÅ‚-
ïÅ‚
2Sy2 ïÅ‚ ìÅ‚ 2Sz2 ÷Å‚ ìÅ‚ 2Sz2 ÷łśł
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
ðÅ‚ ûÅ‚
S = Eg Å"
Sz = B·xb
0,0063Å" Usr Å"Sy Å"Sz
m, a, b, Uh, Q, Dh, Hp, Usr, A, B, Sy, Sz, S K
Rys. 3. Algorytm umo\
liwiajÄ…cy obliczenie stÄ™\
enia zanieczyszczeń S (opis w tekście)
Rys. 3. The algorithm enabling calculation of pollutants concentration S (description in
text)
672 VI Ogólnopolska Konferencja Naukowa
Komputerowa symulacja rozprzestrzeniania zanieczyszczeń w atmosferze
Rys. 4. Wprowadzanie danych
Rys. 4. Data input
Po wybraniu z menu  Zanieczyszczenia opcji  Dane wstępne (F5) na
ekranie monitora wyświetlone zostaje okno przedstawione na rysunku 4.
Moduł umo\
liwiajÄ…cy wprowadzenie danych projektowych podzielony
jest na cztery części charakteryzujące: parametry techniczne emitera, parametry
meteorologiczne, lokalizację receptora oraz tzw. wartości inne. Do parametrów
technicznych zalicza się: wysokość emitera h [m], średnicę wewnętrzną d [m]
wylotu przewodu emitującego substancje, prędkość v [m/s] i temperaturę T [K]
gazów na wylocie z emitera oraz maksymalną emisję Eg [mg/s] substancji ga-
zowej. EmisjÄ™ substancji zanieczyszczajÄ…cych ze z
ródeł projektowanych ustala
siÄ™ szacunkowo na podstawie wskaz
ników emisji tych\
e substancji, charaktery-
stycznych dla procesu technologicznego lub przez analogiÄ™ do emisji ze z
ródła
istniejÄ…cego. W przypadku z
ródeł istniejących i z
ródeł projektowanych o przy-
jętej technologii i rozwiązaniach technicznych, emisję substancji zanieczysz-
VI Ogólnopolska Konferencja Naukowa 673
Jacek Piekarski, Marcin Lubierski
czających ustala się na podstawie wyników pomiarów, zało\
eń technologicz-
nych oraz wielkości produkcji.
Parametry meteorologiczne to: prędkość wiatru Ua [m/s] mierzona na
wysokości anemometru, średnia temperatura powietrza T0 [K] odnosząca się
do okresu obliczeniowego oraz stan równowagi atmosfery ATM [ ]. Lokaliza-
cja receptora podawana jest w formie odległości od emitera: x [m] (równoległa
z kierunkiem wiatru) i y [m] (prostopadła do kierunku wiatru) oraz wysokości,
dla której określa się stę\
enie substancji w terenie z [m]. Parametry inne to
współczynnik aerodynamicznej szorstkości terenu Z0 [m], przy pomocy którego
uwzględnia się wpływ pokrycia powierzchni na intensywność rozpraszania
zanieczyszczeń w atmosferze. Jego wartość wyznacza się na podstawie map
topograficznych.
Podczas edycji poszczególnych wartości parametrów obliczeniowych
program automatycznie dokonuje sprawdzenia ich poprawności, pod względem
przynale\
ności do określonych zakresów, w których równania opisujące zmiany
stÄ™\
enia zanieczyszczeń w atmosferze są słuszne. Błędy sygnalizowane są od-
powiednim komunikatem ukazujÄ…cym siÄ™ w oknie modalnym.
Po prawidłowym wprowadzeniu danych wstępnych nale\
y wybrać
z menu  Zanieczyszczenia opcję  Wyniki obliczeń (F6)  rysunek 5, program
zacznie realizować algorytm przedstawiony na rysunku 3.
Na podstawie wartości stanu równowagi atmosfery ATM [ ], która
zmienia siÄ™ w granicach od 1 do 6, program dobiera parametry meteorologiczne
m [ ], a [ ] oraz b [ ], które są u\
ywane w dalszych obliczeniach. W zale\
ności
od wysokości h [m] emitora oblicza się prędkość wiatru Uh [m/s] na wysokości
wylotu z emitora. W przypadku, gdy obliczona wartość będzie mniejsza zosta-
nie przyjęta Uh = 0,5 [m/s]. Emisja ciepła Q [kJ/s] z emitera obliczana jest
z zale\
ności od średnicy d [m] emitora, prędkości v [m/s] gazów na wylocie
z emitera oraz temperatur T0 [K] otoczenia i T [K] gazów na wylocie z emitera.
Następnie na podstawie wartości emisji ciepła Q [kJ/s] z emitera program do-
biera i wyświetla nazwę formuły, według której następnie oblicza wartość pa-
rametru wyniesienia Dh [m] gazów odlotowych. Dla Q < 16000 [kJ/s] zastoso-
wana zostaje formuła Hollanda, w której w zale\
ności od stosunku prędkości
wiatru Uh [m/s] na wysokości emitera do prędkości v [m/s] gazów na wylocie z
emitera jest obliczane Dh, jako funkcja emisji ciepła Q, średnicy d emitora oraz
prędkości v gazów. Dla Q > 24000 [kJ/s] stosowana jest formuła CONCAWE,
gdzie wartość wyniesienia Dh jest obliczana na podstawie wartości emisji
Q i prędkości wiatru Uh. W przypadku, gdy Q będzie zawierać się w przedziale
wartości od 16000 do 24000 [kJ/s] zastosowane zostaną obie formuły, odpo-
wiednio pomniejszone o ró\
nicę pomiędzy wartością emisji Q a wartościami
granicznymi. Suma wartości parametrów wysokości emitera h [m] oraz wynie-
674 VI Ogólnopolska Konferencja Naukowa
Komputerowa symulacja rozprzestrzeniania zanieczyszczeń w atmosferze
sienia Dh [m] gazów odlotowych stanowi efektywną wysokość emitera Hp [m].
W zale\
ności od efektywnej wysokości Hp emitora oraz wysokości h emitora
program oblicza średnią prędkość wiatru Usr [m/s] w warstwie, której wyso-
kość waha się od z = h do z = Hp. Jej wartość jest obliczana na podstawie efek-
tywnej wysokości emitera Hp, wysokości emitera h oraz prędkości wiatru Ua.
Program, podobnie jak dla Uh, przyjmie Usr = 0,5 [m/s], gdy średnia prędkość
wiatru będzie mniejsza. Następnie dokonywana zostaje weryfikacja ilorazu
efektywnej wysokości emitera Hp [m] oraz współczynnika aerodynamicznej
szorstkości terenu Z0 [m], który powinien zawierać się w przedziale wartości
od 10 do 1500.
Rys. 5. Wyprowadzenie wyników symulacji komputerowej
Rys. 5. Output results of computer simulation
VI Ogólnopolska Konferencja Naukowa 675
Jacek Piekarski, Marcin Lubierski
W kolejnym etapie realizacji algorytmu (rysunek 3), program
oblicza współczynniki A [ ] oraz B [ ], które są zale\
ne od współczynnika
meteorologicznego m oraz stosunku efektywnej wysokości Hp emitora
i współczynnika aerodynamicznej szorstkości terenu Z0. Na podstawie obliczo-
nych współczynników A i B, dobranych wcześniej współczynników empirycz-
nych a i b oraz odległości x [m] receptora od emitora określane zostają współ-
czynniki dyspersji poziomej Sy [m] i pionowej Sz [m]. Ostatnim parametrem,
wyznaczanym przez program jest stÄ™\
enie jedno-godzinne substancji gazowej
S [µg/m3], zale\
ne od wartości emisji Eg [mg/s] gazów, średniej prędkości Usr
wiatru, odległości y [m] emitora od receptora, wysokości z [m] receptora, efek-
tywnej wysokości emitora Hp oraz współczynników dyspersji Sy i Sz.
Rys. 6. Graficzne przedstawienie zmiany współczynnika dyspersji w zale\
ności od
odległości
Rys. 6. Graphic represaentation of dispersion coefficient change in dependence on
distance
676 VI Ogólnopolska Konferencja Naukowa
Komputerowa symulacja rozprzestrzeniania zanieczyszczeń w atmosferze
Dodatkowymi elementami programu są dwa moduły graficzne. Pierw-
szy jest uruchamiany poprzez wybór z menu  Zanieczyszczenia opcji  Dys-
persja (F8) i przedstawia (rysunek 6) rozpraszanie smugi zanieczyszczeń, czyli
zmiany wartości współczynników dyspersji poziomej Sy [m] oraz pionowej Sz
[m] w funkcji zmiany odległości x [m] (równoległej do kierunku wiatru) od
receptora od emitera (przyjęto x " <0 m; 10000 m>). Drugi moduł (rysunek 7)
uruchamiany przy pomocy polecenia  StÄ™\
enie (F9) z menu  Zanieczyszcze-
nia obrazuje graficznie rozkład stę\
enia Sx [µg/m3] substancji gazowej w osi
wiatru, na powierzchni terenu w odległości x [m] od emitora.
Rys. 7. Rozkład stę\
enia zanieczyszczeń w osi wiatru i na powierzchni terenu
Rys. 7. Distribution of pollutant concentration in axis of wind and on the terrain surface
Uzupełniającym elementem programu jest mo\
liwość generowania
raportu w postaci tabelarycznej, podzielonej na cztery kolumny  rysunek 8,
które zawierają informacje odnośnie nazwy parametru, symbolu, jednostki oraz
wartości parametru. Raport podzielony jest pod względem informacyjnym na
VI Ogólnopolska Konferencja Naukowa 677
Jacek Piekarski, Marcin Lubierski
dwie części. Część pierwsza zawiera dane początkowe (wejściowe), natomiast
w części drugiej przedstawione są wyniki obliczeń. Raport mo\
na wydrukować
lub zapisać w postaci pliku: tekstowego (*.txt), w formacie arkusza kalkulacyj-
nego Excel (*.csv) oraz w formacie html (*.html). Program umo\
liwia odczyt
oraz zapis wprowadzonych danych początkowych oraz współczynników obli-
czeniowych do pliku tekstowego (*.koc).
Rys. 8. Raport symulacji komputerowej
Rys. 8. Report of computer simulation
678 VI Ogólnopolska Konferencja Naukowa
Komputerowa symulacja rozprzestrzeniania zanieczyszczeń w atmosferze
3. Przykładowe obliczenia komputerowe
Dane wejściowe (Input data)
Wysokość emitera (od powierzchni terenu) h = 44,0 [m]
Średnica wewnętrzna wylotu przewodu emitujące- d = 1,6 [m]
go substancje
Prędkość gazów na wylocie z emitera v = 14,4 [m/s]
Temperatura gazów na wylocie z emitera T = 493,80 [K]
Emisja NO2 Eg = 43,9 [mg/s]
Prędkość wiatru (na wysokości anemometru) Ua = 13,3 [m/s]
Åšrednia temperatura powietrza dla okresu oblicze- T0 = 270,13 [K]
niowego
Stan równowagi atmosfery ATM = 4 [ ]
Odległość receptora od emitera (równoległa do kie- x = 500 [m]
runku wiatru)
Odległość receptora od emitera (prostopadła do y = 500 [m]
kierunku wiatru)
Wysokość dla której określa się stę\
enie substancji z = 100 [m]
w terenie
Współczynnik aerodynamicznej szorstkości terenu Z0 = 0,18 [m]
Dane wyjściowe (Output data)
Stała obliczeniowa m = 0,27 [ ]
Stała obliczeniowa a = 0,818 [ ]
Stała obliczeniowa b = 0,822 [ ]
Prędkość wiatru na wysokości wylotu z emitera Uh = 18,1 [m/s]
Emisja ciepła Q = 4655 [kJ/s]
Zastosowana formuła Hollanda
Wyniesienie gazów odlotowych Dh = 2,6 [m]
Efektywna wysokość emitera Hp = 46,6 [m]
Średnia prędkość wiatru Usr = 18,3 [m/s]
(w warstwie od z = h do z = HP)
Współczynnik obliczeniowy A = 0,381 [ ]
Współczynnik obliczeniowy B = 0,218 [ ]
Współczynnik dyspersji poziomej Sy = 61,47 [m]
(w warstwie od z = 0 do z = HP)
Współczynnik dyspersji pionowej Sz = 36,06 [m]
(w warstwie od z = 0 do z = HP)
StÄ™\
enie 1-godz. dwutlenku azotu S = 2,5 E-16 [µg/m3]
VI Ogólnopolska Konferencja Naukowa 679
Jacek Piekarski, Marcin Lubierski
Wyniki przeprowadzonych obliczeń stę\
enia 1-godzinnego S [µg/m3]
dwutlenku azotu na powierzchni terenu dla przedstawionych danych wejścio-
wych i w odległościach x " <0 m; 10000 m> oraz y " <-1000 m; 1000 m>
(wartość ujemna składowej y oznacza odległość odmierzaną po lewej stronie
osi wiatru  do obliczeń podstawiana została wartość dodatnia) przy wyko-
rzystaniu programu  Atmo zostały przedstawione w formie graficznej na
rysunku 9.
S [ug/m3]
0,05
0,04
0,03
0,02
0,01
0
-1000
-500
10000
0
8000
6000
Y [m]
500 4000
X [m]
2000
Rys. 9. StÄ™\
enie NO2 obliczone przy pomocy programu  Atmo
Rys. 9. Concentration NO2 calculated by programme  Atmo
4. Podsumowanie
Istnieje wiele programów, przy pomocy których mo\
emy symulować
rozprzestrzenianie się zanieczyszczeń w atmosferze. Programy te umo\
liwiajÄ…
modelowanie na du\
ych obszarach i z uwzględnieniem ró\
norodnych przemian,
którym podlega smuga zanieczyszczeń. Jednak są one kosztowne
680 VI Ogólnopolska Konferencja Naukowa
Komputerowa symulacja rozprzestrzeniania zanieczyszczeń w atmosferze
i skomplikowane w obsłudze. Przedstawiony program komputerowy  Atmo
jest aplikacją darmową (freeware), którą mo\
na pobrać ze strony internetowej
http://www.wbiis.tu.koszalin.pl/jacek. SÅ‚u\
y do obliczania stÄ™\
eń dla niewiel-
kich odległości (do 10 km). Pomaga wyeliminować czasochłonne obliczenia
przy wyznaczaniu stÄ™\
enia zanieczyszczeń gazowych wg zaleceń Rozporządze-
nia Ministra Åšrodowiska i jako taki mo\
e on stanowić pomocne narzędzie.
Literatura
1. Lubierski M.: Numeryczne metody symulacji wybranych modeli rozprzestrzeniania
zanieczyszczeń w powietrzu atmosferycznym. Praca magisterska  promotor dr in\
.
J. Piekarski, Politechnika Koszalińska, Koszalin 2002.
2. Madany A., Bartochowska M.: PrzeglÄ…d polskich modeli rozprzestrzeniania zanie-
czyszczeń atmosferycznych. Prace Naukowe PW, In\
. Åšrodowiska 1995, Z. Nr 19
3. Markiewicz M.T.: Przegląd modeli rozprzestrzeniania zanieczyszczeń w powietrzu
atmosferycznym. Prace Naukowe PW, In\
. Åšrodowiska 1996, Z. Nr 21
4. Nowicki M.: Parametry empiryczne w modelach dyfuzji zanieczyszczeń
w atmosferze. Zeszyt problemowy Nr X/84-85 Ochrona atmosfery z serii in\
ynier-
skie metody badania i obliczania stanu zanieczyszczenia powietrza atmosferyczne-
go, PZITS, Warszawa, 1984-85.
5. Petterssen S.: Zarys meteorologii. PWN Warszawa 1966.
6. Piecuch T. i inni: Spalanie i piroliza odpadów oraz ochrona powietrza przed szko-
dliwymi składnikami spalin. Wydawnictwo Politechniki Koszalińskiej, Koszalin
2002.
7. Piecuch T. i inni: Ochrona powietrza. skrypt uczelniany Politechnika Częstochow-
ska, Częstochowa 1981.
8. Rozporządzenie Ministra Środowiska z dnia 5 grudnia 2002 r. w sprawie wartości
odniesienia dla niektórych substancji w powietrzu, Dz. U. z 2003 r. Nr 1, poz. 12.
9. Szklarczyk M.: Ochrona atmosfery. Wydawnictwo Uniwersytetu Warmińsko-
Mazurskiego, Olsztyn 2001.
10. Woś A.: Meteorologia dla geografów. Wydawnictwo Naukowe PWN, Warszawa
1997.
11. http://www.stronameteo.prv.pl/
12. http://prognoza.icm.edu.pl/
13. http://www.meteo.ids.pl/meu
14. http://www.wunderground.com/global/PL.html
VI Ogólnopolska Konferencja Naukowa 681
Jacek Piekarski, Marcin Lubierski
Streszczenie
Zanieczyszczenia ulatujÄ…ce ze z
ródeł emisji ulegają rozproszeniu
w atmosferze. Ich pomiar jest czasochłonny i kosztowny. Dlatego opracowano
matematyczne modele, które pozwalają z zadowalającą dokładnością określić
stÄ™\
enia zanieczyszczeń. Jednym z takich modeli jest model gaussowski zale-
cany przez RozporzÄ…dzenie Ministra Åšrodowiska z dnia 5 grudnia 2002 r. Auto-
rzy artykułu wykonali i przedstawili charakterystykę programu komputerowe-
go, opracowanego do obliczania stÄ™\
enia zanieczyszczeń, wg wytycznych. Opi-
sali obsługę programu, zakres jego stosowalności oraz przykład. Przedstawiony
program jest prosty w obsłudze, zawiera szereg pomocniczych funkcji
i procedur ułatwiających pracę.
Computer Simulation of Pollutants Propagation
in the Atmosphere
Abstract
The pollutants outflowing from the emission sources undergo propaga-
tion in the atmosphere. Their measurement is time-consuming and expensive.
Therefore the mathematical models which allow to calculate concentration of
pollutants with satisfactory accuracy were worked out. One of such models is
Gaussian (recommended by Directive of Minister of the Environment). The
Authors of the paper created and introduced characteristics of computer pro-
gramme, worked out for calculation of pollutants concentrations, according to
guidelines. They described of programme maintenance, range of its applicabili-
ty as well as example. Introduced programme is easy in operation, contains
several auxiliary functions and procedures facilitating operation.
682 VI Ogólnopolska Konferencja Naukowa