Fluent w zastosowaniu do problemu transportu skażeń w otwartej atmosferze – próby wstępne

Piotr Korczyk
ZMiFP, 2008

1Wstęp i założenia

Praca wykonana w ramach grantu PBZ-MNiSW-DBO-03/I/2007. Celem pracy jest opracowanie metody szybkiej analizy rozprzestrzeniania się skażeń emitowanych do atmosfery. Dla warunków początkowych określonych na podstawie analizy chwilowej sytuacji należy znaleźć rozwiązanie problemu rozpraszania się zanieczyszczenia w otoczeniu.

Na obecnym etapie pracy przyjęto następujące założenia:

Warunki brzegowe

Domena obliczeniowa ograniczona jest od dołu powierzchnią ziemi i znajdującymi się na niej elementami (np. zabudowania). Warunkiem brzegowym nałożonym na te powierzchnie jest brak poślizgu. Reszta ścianek to powierzchnie wirtualne utworzone oby oddzielić domenę obliczeniową od reszty atmosfery. Na ściankach tych należy narzucić prędkość powietrza zgodnie z aktualnym kierunkiem i prędkością wiatru. Należy zadbać również o ustalenie odpowiedniego profilu prędkości wzdłuż składowej pionowej. W skalach czasowych rzędu kilkudziesięciu minut można przyjąć stałość wiatru związanego z wielkoskalowymi ruchami atmosfery. To założenie zostało wykorzystane w niniejszej pracy.

2Siatka i warunki brzegowe

W prezentowanych tutaj przykładowych obliczeniach użyto siatki prostokątnej stworzonej oprogramowaniem GAMBIT. Domena obliczeniowa ograniczona jest do prostopadłościanu o wymiarach:

długość -200m,

szerokość – 50m,

wysokość – 50m.

Domena została tak zorientowana, żeby kierunek wiatru był równoległy do najdłuższego boku(patrz rys 1). Szerokość powinna być dostosowana do dyfuzji skażenia w płaszczyźnie poziomej w kierunku prostopadłym do prędkości wiatru. Wysokość domeny powinna być nie większa od planetarnej warstwy grnicznej

W domenę obliczeniową wstawione zostały prostopadłościenne przeszkody symulujące zabudowania i prostopadłościenny komin jako emiter skażenia. Ważne jest żeby objętości wszystkich dodatkowych przedmiotów wyrzucić z domeny obliczeniowej (narzędzie „subtract”, patrz rys. 3) tak, żeby ostatecznie stworzyć jedną objętość a budynki stanowiły jedynie powierzchnie ograniczające domeną obliczeniową.

Na powierzchni budynków i dolnej powierzchni domeny obliczeniowej (ziemi) ustalono warunki brzegowe typu „wall” bez poślizgu.

W dalszych próbach geometria siatki została dodatkowo skomplikowana przez pofalowanie powierzchni ziemi i rotacje budynków.

Prędkość na pozostałych ściankach domeny obliczeniowej ustalono przyjmując założenie, że wielkoskalowa prędkość wiatru nie zmienia się w trakcie rozchodzenia się zanieczyszczenia.

Na ściankach reprezentujących wlot powietrza ustawiony został profil prędkości, która rośnie wraz z wysokością. Ścianka te są rodzaju „inlet velocity”. Podobnie na ściance górnej ustawiono warunek brzegowy typu „inlet velocity”. Wartości prędkości na tych ściankach są ustawione przez wykorzystanie następującej funkcji UDF:

#include "udf.h"  /* must be at the beginning of every UDF you write */

DEFINE_PROFILE(x_velocity,thread,index)
{
  real x[ND_ND];                /* this will hold the position vector */
  real y;
  face_t f;

  begin_f_loop(f,thread)  /* loops over all faces in the thread passed 
                              in the DEFINE macro argument  */
    {
      F_CENTROID(x,f,thread);
      y = x[1];
      F_PROFILE(f,thread,index) = log(y+1);
    }
  end_f_loop(f,thread)
}

W tym przypadku ustalony został profil logarytmiczny. Funkcja ta może być w prosty sposób zmodyfikowana tak, aby uwzględniać inny profil prędkości.

Na ściankach bocznych ustawiono warunek brzegowy typu „wall” z poślizgiem.

Komin emitujący zanieczyszczenia to prostopadłościan, na którego bocznych ścianach ustawiony jest warunek brzegowy typu „wall” bez poślizgu. Górna ściana symulująca wylot z komina jest typu „inlet-velocity”. Dzięki temu można na niej zadać prędkość wypływu skażenia i jego koncentrację.

3Obliczenia przepływu

Wykonano wstępne obliczenia na wyżej opisanej siatce. Zanieczyszczenie symulowane było przez parę wodną. Para wodna wprowadzany była do przepływu przez komin na którego górnej ściance zadano warunek brzegowy typu „velocity-inlet” i

Wykonano obliczenia dla modeli stacjonarnych i zmiennych w czasie.

4Wizualizacja przepływu

Dobrym sposobem na wizualizację trójwymiarowej struktury chmury zanieczyszczeń jest zobrazowanie powierzchni stałej koncentracji zanieczyszczenia (rys.4 i 5).

Okno tworzenia izopowierzchni: Surface → Iso-Surface. Po utworzeniu odpowiedniej powierzchni powinna pojawić się ona w spisie powierzchni w oknie rysowania konturów (Display → Contours). Możliwe jest wyświetlenie na izopowierzchni innych wielkości np. modułu prędkości (rys. 5).

Inną metodą wizualizacji wyników jest przedstawienie izolinii koncentracji skażenia na powierzchni ziemi (rys. 6). Dzięki temu widać, które miejsca narażone są na największe skażenia.

Wyświetlenie koncentracji skażenie na powierzchni ziemi i ścianach budynków pozwala na analizę skali potencjalnego skażenia w miejscach przebywania ludzi.

5Problemy do rozwiązania

• Generacja siatki – opracować metodę wprowadzania siatki z nierównością powierzchni ziemi.

• Dołączyć do modelu strumienie konwekcyjne zależne od rodzaju podłoża, pory roku i pory dnia.

• Dołączyć chropowatość podłoża w zależności od rodzaju podłoża.

• Włączyć reakcje symulujące odparowywanie kropelek skażenia.

• Zastanowić się nad modelem turbulencji.

• Opracować system analizy zależności rozwiązania od niedokładności warunków początkowych i brzegowych.

6Dodatki

6.1Tworzenie siatki

To jest instrukcja do tworzenia siatki obliczeniowej w programie GAMBIT. Przeznaczona jest dla osób, które są już z tym programem zaznajomione.

Materiały do tej części znajdują się tutaj.

Opisana zostanie procedura tworzenia siatki o dowolnie pofalowanej powierzchni ziemi.

Na początek potrzebne są punkty opisujące topografię, czyli pofalowanie powierzchni ziemi. Zbiór takich punktów mozna stworzyć na podstawie map topograficznych sczytując wysokość niektórych punktów a następnie interpolując te dane na regularną siatkę (na przykład w przy użyciu polecenia griddata w MATLABIE). Regularna siatka jest potrzebna do wczytania do programu GAMBIT. Przykładowy zbiór punktów na regularnej siatce można wygenerować na przy użyciu takiego skryptu w MATLABIE:

function []=tworz_pow



startX=0; endX=200; stepX=1;
startY=0; endY=50; stepY=1;

X=[startX:stepX:endX];
Y=[startY:stepY:endY];

lX=length(X);
lY=length(Y);

X=repmat(X,lY,1);
Y=repmat(Y,lX,1)';

Z=(.03*(Y-25).^2+.1*Y+.01*X+.003*X.*Y+  10*exp(((X-100).^2+(Y-40).^2)/-500)  -15*exp(((X-150).^2+(Y-20).^2)/-400) +20*exp(((X-50).^2+(Y-8).^2)/-1000))/4
Z=Z-min(Z(:))+.1;
M=[X(:),Z(:),Y(:)];

surf(X,Y,Z)
axis equal
shading interp

save gorki.txt M -ASCII

W tym przypadku zostanie stworzona powierzchnia opisywana odpowiednią funkcją na siatce o wymiarach 201 na 51 punktów. Dane zostaną zapisane do pliku gorki.txt.

Tak utworzony zbiór punktów należy wczytać do programu GAMBIT.

Po wczytaniu plików z uprzednio utworzonego pliku należy je następnie połączyć, żeby GAMBIT rozumiał je jako powierzchnię.

Po utworzeniu powierzchni można wyrzucić wszystkie punkty.

Teraz trzeba utworzyć prostopadłościan, z którego utworzymy domenę obliczeniową.

Teraz należy przeciąć objętość prostopadłościanu wcześniej utworzoną powierzchnią ziemi.

Ta operacja spowoduje, że objętość prostopadłościanu zostanie podzielona na dwie objętości. Jedną z nich – tę która znajduje się pod powierzchnią ziemi należy wyrzucić.

Teraz należy ewentualnie dodać przeszkody takie jak zabudowania. Dodajmy zwykły prostokąt, który będzie reprezentował pojedyncze zabudowanie.

Teraz obracamy prostopadłościan o dowolny kąt (w tym przypadku 40^o).

Po obróceniu możemy przemieścić nasze zabudowanie.

Kolejnym krokiem jest wyrzucenie objętości zabudowania. Przepływ będzie liczony na zewnątrz budynków i innych przeszkód, więc ich wnętrze musi zostać usunięte z domeny obliczeniowej. Idea dodawania kolejnych elementów polega więc na wyrzucaniu ich objętości z domeny obliczeniowej. Ta operacja modyfikuje powierzchnię, która ogranicza domenę. W ten sposób otrzymujemy tylko jedną objętość w której utworzona zostanie siatka obliczeniowa.