Wstęp teoretyczny
Czynniki determinujące rozprzestrzeniania się zanieczyszczeń w atmosferze [1]:
1.Czynniki meteorologiczne: a)turbulencja atmosfery,
b)pionowy gradient temperatury,
c)poziomy ruch dużych mas powietrza - wiatr, d)grubość strefy mieszania,
e)opady atmosferyczne. 2.Czynniki topograficzne:
a)naturalne ukształtowanie podłoża – tzw. rzeźba terenu, b)czynniki naturalne np. lasy, zbiorniki wodne itp.,
c)czynniki antropogeniczne np. zabudowa miejska, wyrobiska pokopalne. 3.Czynniki technologiczne:
a)wysokość źródła emisji,
b)prędkość wylotu zanieczyszczeń z komina, c)temperatura emitowanych gazów lub pyłów.
Typy źródeł emisji zanieczyszczeń [2]: 1.Źródła punktowe np. komin,
Źródła o wymiarach znacznie mniejszych od rozpatrywanej odległości transportu zanieczyszczenia.
2.Źródła liniowe np. ulica,
Emisja zanieczyszczenia odbywa się wzdłuż krzywej o długości porównywalnej z rozpatrywaną odległością transportu zanieczyszczenia. 3.Źródła powierzchniowe np. rozlewisko niebezpiecznej substancji,
powstałe w wyniku awarii zbiornika, wyodrębnione obszary miast.
Wymiary płaszczyzny emisji zanieczyszczenia nie mogą być zaniedbane, ze względu na porównywalną wielkość w stosunku do rozpatrywanej odległości transportu zanieczyszczenia.
4.Źródła objętościowe.
Ponadto źródła mogą charakteryzować się stałymi lub zmiennymi w czasie parametrami emisji, a samo źródło może mieć charakter stacjonarny (np. komin), lub też może poruszać się np. samochód.
Podział modeli rozprzestrzenia zanieczyszczeń w powietrzu [2,3].
Modele fizyczne – modele w mniejszej skali, realizowane najczęściej w warunkach laboratoryjnych, w których próbuje się odtworzyć sytuację rzeczywistą mającą miejsce w powietrzu atmosferycznym – np. tunele wiatrowe. Mogą być one wykorzystywane do weryfikacji prognoz oraz kalibracji modeli matematycznych.
Modele matematyczne:
a)deterministyczne – wykorzystują matematyczny opis zjawisk zachodzących w atmosferze, np. równanie transportu zanieczyszczeń i ich depozycji, równania przemian chemicznych. W zależności od podejścia do rozwiązania zestawu równań matematycznych wyróżnić można modele analityczne (wykorzystujące rozwiązanie dokładne, zazwyczaj przy szeregu uproszczeń i założeń) oraz numeryczne (stosują przybliżone rozwiązanie równań różniczkowych, zazwyczaj przy zastosowaniu metod komputerowych – metoda
2
różnic skończonych, metoda elementu skończonego, metoda objętości skończonych).
b)stochastyczne – do generowania prognoz wykorzystują statystyczne zależności pomiędzy informacjami odnośnie źródeł emisji, a stanem zanieczyszczenia powietrza atmosferycznego w danym punkcie, który znany może być np. z pomiarów. W modelach stochastycznych brak bezpośredniego opisu zjawisk i procesów mających miejsce w atmosferze. Wykorzystywane są niejawne zależności pomiędzy emisją, warunkami meteorologicznymi, a imisją (stanem aerosanitarnym powietrza), dlatego też modele stochastyczne mogą być stosowane tylko dla tych przypadków, których dotyczyły dane wejściowe wykorzystywane do ich tworzenia.
Literatura:
[1] Orzechowski Z., Prywer J., Zarzycki R., Mechanika płynów w inżynierii środowiska, WNT Warszawa, 1997, str. 433-452.
[2] Markiewicz M. T., Podstawy modelowania rozprzestrzeniania się zanieczyszczeń w powietrzu atmosferycznym, Oficyna Wydawnicza Politechniki Warszawskiej, 2004.
[3] Holnicki-Szulc P., Modele propagacji zanieczyszczeń atmosferycznych w zastosowaniu do kontroli i sterowania jakością środowiska, Akademicka Oficyna Wydawnicza EXIT, Warszawa, 2006.
Równania z mechaniki płynów wykorzystywane przy opisie rozprzestrzeniania zanieczyszczeń w atmosferze w modelach matematycznych.
Równania Naviera-Stokesa(płyn lepki, nieściśliwy)
dla kierunku X:
dla kierunku Y:
dla kierunku Z:
vx, vy, vz – składowe prędkości w kierunku x, y, z [m/s] t – czas [s]
g – przyspieszenie w kierunku x, y, z [m2/s] µ - lepkość płynu [Pa·s]
ρ – gęstość płynu [kg/m3]
µt – lepkość turbulentna [Pa·s]
Równania ciągłości(płyn nieściśliwy)
+ ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ x p g z v v y v v x v v t v x x z x y x x x
ρ
ρ
(
)
(
)
(
)
∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ x v z v z x v y v y x v x z x t y x t x tµ
µ
µ
µ
µ
µ
2 + ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ y p g z v v y v v x v v t v y y z y y y x yρ
ρ
(
)
(
)
(
)
∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ y v z v z y v y y v x v x z y t y t x y tµ
µ
µ
µ
µ
µ
2 + ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ z p g z v v y v v x v v t v z z z z y z x zρ
ρ
(
)
(
)
(
)
∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ z v z z v y v y z v x v x z t y z t x z tµ
µ
µ
µ
µ
µ
0 = ∂ ∂ + ∂ ∂ + ∂ ∂ z v y v x vx y zModel turbulencji np. K-model (inaczej: równanie transportu energii
kinetycznej turbulencji)
K – energia kinetyczna turbulencji (prędkości fluktuacyjnych) σK, CK – parametry modelu
lµ - droga mieszania
γ
&– lokalna szybkość ścinania w płynieInne modele turbulencji to np. k-ε, k-ω. Poza tym próbuje się stosować metodę symulacji wielkich wirów (LES - Large Eddy Simlution) oraz bezpośrednią symulację turbulencji (DNS - Direct Numerical Simulation) – tym niemniej wzrost zapotrzebowania na moc obliczeniową dla dwóch ostatnich przypadków limituje możliwości ich praktycznego wykorzystania.
Problem zadowalającego matematycznego opisu turbulencji, pomimo wielu istniejących modeli, pozostaje nadal nierozwiązany !
Równanie transportu zanieczyszczeń – równanie adwekcji-dyfuzji
c – stężenie [kg/m3]
∆ – szybkość uwalniania zanieczyszczeń [kg/s] D – współczynniki dyfuzji molekularnej [m2/s] Dtx, D
t y , D
t
z – współczynniki dyspersji burzliwej [m 2
/s]
W powyższym równaniu obok źródła zanieczyszczenia ∆ istnieje również możliwość uwzględnienia zaniku zanieczyszczenia jako dodatkowego składnika. Grupa modeli wykorzystujących takie równanie nazywana jest tzw. modelami teorii K, w których Kx,y,z=D+D
t
x,y,z. K są anizotropowymi
współczynnikami dyfuzji uwzględniającymi wpływ dyspersji burzliwej. Trudności w zastosowaniu tych modeli spowodowane są koniecznością opracowania równań opisujących zmienność parametru K.
Modele gaussowskie jako przykład modeli z dokładnym rozwiązaniem analitycznym.
Modele gaussowskie stosujące tzw. równanie Pasquilla są grupą modeli najczęściej wykorzystywaną w praktyce ze względu na ich prostotę obliczeniową i łatwość zastosowania. Podstawowym założeniem prowadzącym do otrzymania równań gaussowskich jest przyjęcie, że rozkład stężeń zanieczyszczenia w smudze jest zgodny z rozkładem Gaussa, a uzyskiwane pola stężeń ustalone w czasie.
+ − + − − − − = 2 2 2 0 2 1 exp 2 1 exp 2 1 exp 693 , 0 exp 2 ) , ( z e z e y z y H z H z y u Q r s c δ δ δ τ τ δ πδ = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ z K v y K v x K v t K z y x
ρ
µρ
γ
µ
σ
µ
µ
σ
µ
µ
σ
µ
µ
l K C z K z y K y x K x K t D t K t K t 2 3 2 − + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ = & µ µρ
µ
t = C K1/2l δ µ =0.4 l = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ z c v y c v x c v t c z y x(
)
(
)
(
)
+∆ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ z c D D z y c D D y x c D D x t z t y t x4
u - średnia prędkość wiatru
He - efektywna wysokość emisji (suma wysokości emitora i wyniesienia smugi)
Q - natężenie emisji zanieczyszczenia
c(s,r)- stężenie zanieczyszczenia w punkcie r wynikające z emisji w punkcie s
δy, δz - współczynniki poziomej i pionowej dyspersji, ustalane jako funkcje pionowych stanów turbulencji i odległości receptora od źródła emisji
τ0 - czas połowicznego zaniku danej substancji
τ - rozpatrywany czas rozprzestrzeniania się zanieczyszczenia
Modele gaussowskie nie mogą być stosowane w przypadku małych prędkości wiatru i w warunkach bezwietrznych, ponadto nie uwzględniają występowania przeszkód terenowych i wymagają oszacowania parametrów δy, δz w oparciu o stany równowagi atmosfery.
Rozkład stężeń zanieczyszczenia w smudze zanieczyszczenia zgodny z rozkładem Gaussa [3]
Zadania laboratoryjne do wykonania
ALOHA - Areal Locations of Hazardous AtmospheresZastosowanie programu ALOHA, jako przykładowego modelu matematycznego wykorzystującego dokładne rozwiązanie analityczne, do oszacowania stref zagrożenia dla ludności w przypadku awaryjnych uwolnień substancji niebezpiecznych, dla dwóch scenariuszy zdarzenia:
• scenariusz 1 – wyciek benzenu ze zbiornika w zakładzie przemysłowym,
• scenariusz 2 – emisja chloru w wyniku kolizji na bocznicy kolejowej. Scenariusz zdarzenia 1 [3]
Na terenie zakładów przemysłowych zlokalizowanych na obrzeżach Baton Rouge w Luizjanie (USA) 20 wrześnie 2006 roku strażnik ochrony o godzinie 22.30 zauważa wyciek ze zbiornika, w którym przechowywany jest benzen. Wyciek następuje przez okrągły otwór o średnicy 15 cm, który powstał w płaszczu zbiornika, na wysokości 25 cm od dna. Wypływający z uszkodzonego zbiornika ciekły benzen, gromadzi się na wybetonowanym placu tworząc rozlewisko. Na terenie gdzie nastąpiło zdarzenie znajduje się kilka
niewysokich budynków, natomiast obszar na północ od zakładu pokryty jest niską roślinnością trawiastą.
Warunki meteorologiczne w trakcie zdarzenia:
• temperatura powietrza 27°C,
• wilgotność powietrza 75%,
• wiatr południowo-zachodni, 3 m/s na wysokości 10 m (dane z wieży meteorologicznej),
• zachmurzenie duże - 7/10 nieba pokryte przez chmury,
• z południowego-zachodu nadciąga burza,
• nie ma inwersji temperatury na niskim pułapie.
Wg informacji pozyskanych od strażnika zbiornik o pojemności 2 m3 i średnicy 1,2 m napełniany był całkowicie w godzinach popołudniowych w dniu wystąpienia zdarzenia.
1) Określić zasięg toksycznej chmury, analizując następujące poziomy stężeń zanieczyszczenia: 1000 ppm, 150 ppm, 50 ppm. Porównać wyniki dla klasycznego modelu gaussowskiego i modelu gazu ciężkiego.
2) Określić zagrożenie promieniowaniem cieplnym w przypadku gdyby w wyniku burzy nastąpił zapłon rozlanego benzenu tworząc płonące rozlewisko.
6
Scenariusz zdarzenia 2 (do samodzielnego rozwiązania)
Na bocznicy kolejowej, w odległości około 550 m na zachód od głównego dworca kolejowego w Bytomiu, stoi wagon transportujący ciekły chlor (Cl2). Na skutek błędu maszynisty lokomotywy manewrowej dochodzi do kolizji, w wyniku której następuje całkowite zniszczenie trzech pojemników zawierających po 50 kg chloru każdy. Zaistniała sytuacja ma miejsce w
ciepłym okresie późnowiosennym, we wczesnych godzinach
popołudniowych - hipotetyczny czas wystąpienia zdarzenia: 4 czerwca 2010, godzina 1500. Zdarzenie ma miejsce na wysokości 273 m n.p.m. Należy określić zasięg stref zagrożenia dla ludności Bytomia wykorzystując klasyczny model gaussowski oraz model traktujący uwolnioną substancję jako gaz cięższy od powietrza.
Na potrzeby symulacji zaproponowano następujące warunki meteorologiczne:
• kierunek wiatru: zachodni (W),
• prędkość wiatru: 2 m/s (pomiar na wysokości 3 m od podłoża),
• temperatura powietrza: 22°C,
• wilgotność względna powietrza: 80%,
• zachmurzenie: 3/10 powierzchni nieba pokryte jest chmurami,
• położenie wagonu transportowego z chlorem: E 18 54ۥ 23ۥۥ; N 50 20ۥ 36ۥۥ, (zlokalizować za pomocą planu miasta prezentowanego w Internetowym Serwisie Bytomskiej Infrastruktury Informacji Przestrzennych
http://sitplan.um.bytom.pl/ibidp/).
Analiza okoliczności zaistniałego zdarzenia
Analiza ortofotomapy (z 2008 roku) wykazuje, że bezpośrednie otoczenie wagonu z chlorem stanowi głównie teren otwarty. W kierunku wschodnim, najbliższe budynki oddalone są o około 550 m od miejsca zdarzenia. W kierunku zachodnim, w odległości około 100 m występują tylko pojedyncze, małe budynki, a dalej znajduje się wolna przestrzeń (tory kolejowe). Taki układ na osi wchód-zachód pozwala wnioskować, że zmierzona na potrzeby symulacji prędkość wiatru (2 m/s) nie powinna zmieniać się istotnie w najbliższym sąsiedztwie miejsca wystąpienia uwolnienia awaryjnego. W oparciu o te informacje zaproponować współczynnik szorstkości terenu z0.
Kolejność postępowania
• podać informacje o lokalizacji miejsca zdarzenia (SiteData Location Add);
• zdefiniować czas wystąpienia zdarzenia (SiteData Data&Time);
• wybrać rodzaj uwolnionej substancji (SetUp Chemicals Chlorine);
• zdefiniować warunki meteorologiczne (SetUp Atmospheric User Input);
• określić parametry źródła emisji (SetUp Source Direct)
• wybrać:
a)klasyczny model gaussowski (SetUp Calculation Options Use Gaussian dispersion only),
b)model gazu ciężkiego ((SetUp Calculation Options Use Heavy Gas dispersion only).
• przeanalizować uzyskane wyniki.
Uzyskany zasięg stref zagrożenia zapisać do pliku programu MS Word (w programie ALOHA: Edit Copy, a następnie wkleić rysunki do programu MS
Word jako mapy bitowe: Edycja Wklej specjalnie Mapa bitowa) i sformułować wnioski odnośnie wykonanego ćwiczenia.
Literatura:
[3] The CAMEO® Software System - ALOHA®, User’s Manual, U.S. Environmental Protection Agency & National Oceanic and Atmospheric Administration.
http://www.epa.gov/emergencies/docs/cameo/ALOHAManual.pdf
ANSYS – moduł FLOTRAN CFD
Zastosowanie modułu FLOTRAN programu ANSYS, jako przykładowego modelu matematycznego wykorzystującego numeryczne rozwiązanie równań opisujących rozprzestrzenianie zanieczyszczeń metodą elementu skończonego, do symulacji w czasie ewolucji pola stężeń substancji emitowanej z komina.
Rozprzestrzeniania zanieczyszczeń emitowanych z komina
• schemat problemu do rozwiązania
vwiatr= 3 m/s
d = 1 m H = 5 m vwylot = 1 m/s
• określenie dziedziny obliczeniowej oraz dyskusja warunków brzegowych i początkowych v = 3 v = 0 x y v = 0 v = 0 x y v = 0y v = 0 v = 1 a = 0 a = 1 x y SP01 SP02 a = 1 a = 0 SP01_start SP02_start 45 m 5 m 1 m 79 m 20 m Wykonanie ćwiczenia 1.Przygotowania zadania:
• skopiować z serwera zawartość katalogu CONC_START, do folderu
D:\Zan_pow_lab (w przypadku braku folderu należy go utworzyć);
• uruchomić program Mechanical APDL Product Launcher:
8
• ustawić katalog roboczy programu ANSYS:
Working Directory przycisk Browse i wskazać folder D:\Zan_pow_lab;
• wprowadzić nazwę dla zadania nad którym pracujemy: Job Name i podać nazwę komin_conc;
• uruchomić program ANSYS: przycisk Run
• wczytać siatkę dziedziny obliczeniowej: przycisk RESUME_DB;
• zapisać projekt – użyć przycisku Save Analysis z paska narzędzi standardowych;
2.Określenie przepływu mas powietrza
Zestaw plików pobranych z serwera zawiera zdefiniowaną dziedzinę obliczeniową oraz rozwiązanie pola przepływu mas powietrza. W celu uzyskania wektorów prędkości (w każdym węźle siatki dziedziny obliczeniowej) wykonane zostały następując etapy:
• nałożono warunki brzegowe dotyczące prędkości,
• określono właściwości fizykochemiczne powietrza w oparciu o bazy danych dostępne w programie ANSYS,
• zdefiniowano charakter modelowanego zjawiska (przepływ turbulentny, proces ustalony),
• przeprowadzono obliczenia numeryczne.
2.1. Postprocessing – wizualizacja uzyskanego pola prędkości
• wczytać zbiór wyników obliczeń:
General Postproc Read Results Last Set (dwukrotnie zamknąć okno pojawiającego się komunikatu przyciskiem Close)
General Postproc Plot Results Vector Plot Predefined i zatwierdzić przyciskiem OK proponowane ustawienia w oknie Predefined
Vector Plot of Predefined Vectors;
3.Obliczenia pola stężeń
W rozważanym przykładzie pole stężeń wyrażone będzie w postaci ułamków masowych dwóch składników: powietrza (SP01) oraz emitowanego zanieczyszczenia (SP02). Ułamek masowy powietrza, liczony jest poprzez algebraiczne dopełnienie do wartości 1 sumy ułamków masowych dla SP01 i SP02 w oparciu o rozwiązania równań transportu i dyfuzji dla zanieczyszczenia.
Zestaw plików pobranych z serwera zawiera zdefiniowane warunki brzegowe dotyczące pola stężeń. Określone zostały także parametry fizykochemiczne analizowanych składników. Jako gęstość powietrza przyjęto 1.252 kg/m3 (dla temp. 0°C) oraz założono dla uproszczenia, że emitowane zanieczyszczenie jest mieszaniną gazów o właściwościach takich samych jak otaczające powietrze. Ponadto zainicjowany został także moduł rozwiązania równań transportu, dzięki czemu możliwe jest uruchomienie etapu obliczeniowego, po uprzednim określeniu warunków początkowych.
• wcisnąć przycisk RESUME_DB;
• zdefiniować warunek początkowy do obliczeń w całej dziedzinie obliczeniowej; w chwili początkowej ułamek masowy powietrza wynosi 1, a zanieczyszczenia 0:
Preprocesor Loads Define Loads Apply Initial Condit’n Define, a następnie w oknie Define Initial Cond… wcisnąć przycisk Pick
10
Initial value of DOF podać 1 (OK); następnie powtórzyć czynności dla SP02 (zanieczyszczenie), podając Initial value of DOF równą 0;
3.1. Etap obliczeniowy
W etapie obliczeniowym w kolejnych krokach czasowych rozwiązywane jest pole stężeń w oparciu o ustalone wcześniej pole prędkości. Obliczenia stężeń wykonywane dla stanu nieustalonego, aby pokazać ewolucję smugi zanieczyszczenia w czasie. Symulowany jest czas rzeczywisty wynoszący 1 min. (60 sekund), a w trakcie obliczeń zapisywane są wyniki pośrednie co 10 sekund czasu rzeczywistego.
• uruchomić moduł obliczeniowy: Solution Run Flotran
(UWAGA! Czas obliczeń może wynosić kilka minut. Ze względu na duże obciążenie procesora zalecane jest aby nie uruchamiać w tym czasie innych aplikacji);
3.2. Postprocessing - wizualizacja uzyskanego pola stężeń
• wczytać zbiór wyników obliczeń:
General Postproc Read Results Last Set;
• wyświetlić kontury uzyskanego pola stężeń w postaci ułamków masowych: General Postproc Plot Results Contour Plot Nodal Solution i w oknie Contour Nodal Solution Data rozwinąć kliknięciem opcję DOF Solution i wybrać podopcję Mass fraction of species SP02, następnie zatwierdzić przyciskiem OK;
• Dostosować sposób wyświetlania konturów, tak aby jak najlepiej zobrazować pole stężeń uwalnianego zanieczyszczenia:
Menu PlotCtrls Style Contours Uniform Contours
Podpowiedź: zdefiniować kontury użytkownika (wybrać opcję User Specified), podając wartość minimalną, maksymalną i liczbę wyświetlanych konturów;
• zamknąć program ANSYS;
3.3. obliczenia pola stężeń dla różnych współczynników dyfuzji/dyspersji
• skopiować z serwera zawartość katalogu CONC_START do folderu
D:\Zan_pow_lab2 (w przypadku braku folderu należy go utworzyć);
• uruchomić program Mechanical APDL Product Launcher:
Start Programy ANSYS 12.0 Mechanical APDL Product Launcher;
• ustawić katalog roboczy programu ANSYS:
Working Directory przycisk Browse i wskazać folder D:\Zan_pow_lab2;
Job Name i podać nazwę komin_conc;
• uruchomić program ANSYS: przycisk Run
• wczytać siatkę dziedziny obliczeniowej: przycisk RESUME_DB;
• zapisać projekt – użyć przycisku Save Analysis z paska narzędzi standardowych;
• ustalić wartość współczynnika dyfuzji – w dotychczasowych obliczeniach stosowana była wartość domyślna wynosząca 1:
Preprocesor FLOTRAN Set Up Multiple species, potwierdzić OK i w oknie Multiple species wybrać Species #2, a jako Select data to be
defined/edited wybrać Properties, zatwierdzić OK, następnie w oknie Main
Species Property Selection wybrać Mass Diffusion (OK); dalej w
pojawiającym się oknie Species Mass Diffusion Coef jako Nominal Value wpisać wybraną wartość współczynnika dyfuzji (np. 1e-1), zatwierdzić OK, wybrać opcję EXIT Properties Panel, zatwierdzić OK, zamknąć okno ustawień przyciskiem Cancel;
• zdefiniować warunek początkowy do obliczeń w całej dziedzinie obliczeniowej; w chwili początkowej ułamek masowy powietrza wynosi 1, a zanieczyszczenia 0:
Preprocesor Loads Define Loads Apply Initial Condit’n Define, a następnie w oknie Define Initial Cond… wcisnąć przycisk Pick
all ; w oknie Define Initial Conditions wybrać SP01 (powietrze) i jako
Initial value of DOF podać 1 (OK); następnie powtórzyć czynności dla SP02 (zanieczyszczenie), podając Initial value of DOF równą 0;
• uruchomić obliczenia dla pola stężeń (wszystkie ustawienia, w tym krok czasowy i ilość iteracji są ustawione na takie same wartości jak w poprzedniej części ćwiczenia):
Solution Run Flotran
• wyświetlić wyniki jak w poprzednim etapie ćwiczenia.
D=1 D=1e-1
