Metoda Objętości Skończonych

1. Metoda Objętości Skończonych (ang. Finite Volume Method, FVM)

2. Rozpoczniemy od dwóch podstawowych przykładów wprowadzających metodę skończonych objętości. Przykład 1. (Równanie transportu). Rozważmy równanie gdzie: v(x,t) jest danym polem prędkości;  jest pewnym obszarem płaszczyzny, w którym rozważamy równanie; u0(x) jest daną wartością poczatkową szukanej gestości u(x,t). W przypadku powyższego równania strumien dany jest przez a samo równanie jest wyrazem prawa zachowania (bilansu) wielkości opisanej gęstością u(x,t).

3. Niech oznacza podział obszaru  na elementy w postaci wielokątów wypukłych (np. trójkąty, prostokąty) w ten sposób, że każde dwa elementy albo są rozłączne, albo datykają się tylko wierzchołkami (v), albo mają tylko wspólny bok (e): Elementy zbioru  nazywamy „objętościami kontrolnymi”. Dla dowolnego K   po scałkowaniu równania bilansu otrzymujemy: co po skorzystaniu ze wzoru Gaussa (tw. o dywergencji) daje: gdzie (x) oznacza wektor normalny do brzegu K w punkcie xK.

4. Teraz musimy tą całkową postać prawa zachowania (poprzedni slajd) przekształcić do jakiejś formy zdyskretyzowanej, czyli możliwej do obliczeń numerycznych na komputerze. Niech k > 0 oznacza stały krok czasowy oraz tn=nk, nN. Zapisując równanie całkowe dla tn i dykretyzując pochodną czasową wg schematu jawnego Eulera otrzymujemy Następnie do każdej objętości skończonej (kontrolnej) przypisujemy nieznaną wartość przybliżenia rozwiązania. Najczęściej rozważa się tzw. centralne objętości skończone (ang. cell-centered finite volume method), w których każda nieznana wartość rozwiązania jest przypisana do całej komórki. Wartości te oznaczamy następująco: Ponadto niech K oznacza zbiór wszystkich krawędzi (boków) komórki  .

5. Całkę po brzegu możemy teraz zapisać w postaci sumy: K  (x) v(x,t) gdzie oznacza wektor normalny do  skierowany na zewnątrz komórki K. Wprowadzamy dodatkowe oznaczenia Są to wielkości w zasadzie znane, gdyż nie występuje tutaj niewiadoma funkcja. Całka zależy tylko od  (brzegu komórki elementarnej) oraz od zadanego pola prędkości v(x,tn).

6. Teraz jest kluczowy moment w metodzie objętości skończonych: dyskretyzacja strumieni na brzegu, czyli całek Niech L i K oznaczaja dwie komórki kontrolne, których wspólnym brzegiem jest krawędź : Strumień dyskretyzujemy następująco: Zgodnie z tą definicją mamy: Oznacza to, że nasz schemat przybliżania jest zachowawczy: strumienie dyskretne na brzegu bilansują się. Okazuje się, że jest to kluczowa własność. Dzięki niej można udowodnić matematycznie zbieżność schematu i jego numeryczną stabilność.

7. Całkę aproksymujemy natomiast w oczywisty sposób: ponieważ przybliżona wartość związana z komórką K jest oznaczona przez unK, więc gdzie |K| oznacza pole powierzchni komórki K (w ogólnym przypadki d-wymiarową objętość, gdy rozważamy równanie w Rd). Zatem Uwaga: Czasami miarę |K| zbioru oznacza się inaczej, np. m(K) lub (K).

8. Ostatecznie wyprowadziliśmy następujący schemat dyskretyzacji w metodzie objętości skończonych: dla danego podziału obszaru  szukamy przybliżonego rozwiązania, określonego przez wartości które spełniają następujący układ równań gdzie: Powyżej EKoznacza zbiór krawędzi (w ogólności – ścian) wielokąta (w ogólności – d-wymiarowego wielościanu) K.

9. Przykład 2. (Stacjonarny stan dyfuzji – równanie Poissona). Rozważmy równanie: W równaniu tym u(x) może oznaczać np. stężenie substancji, temperaturę lub potencjał elektryczny (w stanie stacjonarnym). Niech tym razem  oznacza siatkę zbudowaną z prostokątów. Podstawowe założenie dotyczące regularności podziału jest jak zwykle następujące: gdzie v oznacza wierzchołek, a e oznacza krawędź prostokąta.

10. Całkowanie równania Po ustalonym elemencie kontrolnym K i skorzystanie z twierdzenia o całkowaniu przez części lub twierdzenia Gaussa daje Dla każdej komórki kontrolnej K oznaczamy pzrzez xK jej środek, natomiast  oznacza wspólna krawędź dwóch przylegających komórek:  = K  L. Jednym ze sposobów przybliżenia strumienia jest użycie różnicy centralnej: gdzie (uK)K są szukanymi przybliżeniami rozwiązania u(x). Tzn. z każdą komórką K wiążemy wartość uK, którą możemy interpretować jako przybliżenie rozwiązania w środku komórki xK K. Ponadto d oznacza odlęgłość pomiędzy xK a xL.

11. W przypadku, gdy objętość kontrolna K ma krawędź wspólną z fragmentem brzegu  dziedziny , to aproksymacja strumienia jest podobna ale tym razem d oznacza odległość pomiędzy środkiem xK komórki K a brzegiem , natomiast g oznacza wartość funkcji warunku brzegowego w środku odcinka   . W obu przypadkach |  | oznacza długość (ogólnie: (d-1)-wymiarową objętość) odcinka . Stąd otrzymujemy dyskretyzację metodą skończonych objętości dla naszego stacjonarnego problemu: szukamy (uK)K, które spełniają następujący układ równań gdzie

12. Ogólny opis metody objętości skończonych dla równań praw zachowania Metoda objętości skończonych jest używana do dyskretyzacji praw zachowania (równań bilansu). Lokalna (różniczkowa) postać prawa zachowania w d-wymiarowej przestrzeni Rd: W równaniu tym q reprezentuje gęstość przestrzenną jakiejś wielkości (energia, masa, liczba moli jakiegoś składnika), pole wektorowe J(x,t) jest strumieniem tej wielkosci [m-2s-1], a funkcja S(x,t) wyraża możliwe źródła. Symbol div oznacza dywergencja przestrzenną, która w układzie kartezjańskim (prostokątnym) może być wyrażona poprzez pochodne cząstkowe następujaco: gdzie

13. Dyskretyzacja czasu w równaniu bilansu dokanujemy wprowadzając rosnący ciąg Dalej dla prostoty będziemy zakładać stały krok czasowy: tn=nk, dla ustalonego k > 0. Można używać różnych schematów dla pochodnej czasowej. Najprostrzym jest różnica w przód, czyli schemat Eulera: Aby dokonać dyskretyzacji przestrzennej wprowadzamy siatkę  w dziedzinie   Rd składającą się z objętości kontrolnych K. Ogólnie siatka (zwana też czasem triangulacją) składa się z wielościanów wypukłych d-wymiarowaych. Wymagamy, aby: gdzie s jest wspólną ścianą wymiaru  d-1.

14. Z każdym elementem K i dyskretnym czasem tn wiążemy przybliżoną wartość unK szukanego rozwiązania. Pierwszym krokiem klaycznej metody objętości skończonych jest jest całkowanie równania bilansu po każdym elemencie kontrolnym: gdzie vK(x) jest wektorem normalnym zewnętrznym do Kw punkcie x. Drugim krokiem metody jest aproksymacja strumieniaJ(qn,x,tn)·vK(x) przez brzeg K każdej objętości kontrolnej wyrażona przy pomocy {unK: K}. Gdy pochodną czasową przybliżamy różnicą wstecz (niejawna metoda Eulera), to aproksymacja strumienia będzie zawierała niewiadome {un+1K: K}.

15. Wyrażając drugi krok bardziej precyzyjnie: dokonujemy aproksymacji członu strumieniowego przy pomocy pewnej wielkości zależnej od dla metody jawnej jednokrokowa lub od dla metody niejawnej jednokrokowej. Podkreślmy najważniejszą cechę klasycznej metody objętości skończonych: lokalna zachowawczość dla każdej komórki kontrolnej (objętości skończonej) K, tzn. co oznacza, że strumień zdyskretyzowany przez ścianę  = K  L liczony od K do L jest taki sam jak od L do K ze znakiem przeciwnym.

16. Szczegółowy opis jednowymiarowego problemu eliptycznego Niech będzie dana funkcja f: (0, 1)  R. Rozważmy następujące równanie różniczkowe Zauważmy, że równanie –u’’(x)=f(x) może być zapisane w formie zachowawczej: Aby skonstruować procedurę do numerycznego rozwiązania tego problemu musimy najpierw zdefiniować specjalną siatkę  w zbiorze  = (0, 1) składającą się z N komórek (objętości kontrolnych), które będziemy oznaczać przez Ki.

17. Definicja. (Dopuszczalnej siatki w 1D) Dopuszczalna siatka na odcinku  = (a, b), oznaczana przez , dana jest przez rodzinę odcinków (Ki)i taką, że Ki = (xi-1/2, xi+1/2), oraz punkty (xi)i takie, że Ponadto przyjmujemy oznaczenia: Dalej będziemy dla uproszczenia zapisu przyjmowali odcinek jednostkowy:  = (a, b)=(0, 1)

18. Rozwiązaniem problemu jest oczywiście pewna funkcja u: [0, 1]  R, ciągła w zbiorze domnkinętym [0, 1], różniczkowalna w (0, 1) i spełniąjąca odpowiednie równości. Przybliżone wartości rozwiązania w punktach xi, czyli u(xi), oznaczamy przez ui. Równanie całkujemy po elemencie Ki = (xi-1/2, xi+1/2), co daje Teraz musimy dokonać aproksymacji strumienia –u’(x) na brzegu {xi+1/2} = Ki Ki+1 (oczywiście w przypadku 1D brzeg obszaru to jeden punkt i moglibyśmy po prostu powiedzieć, że aproksymujemy pochodną w punkcie). Naturalny wybór to pochodna centralna, zatem Ponieważ dotyczy to węzłów wewnętrznych, więc

19. Zauważmy, że ta aproksymacja jest zgodna w tym sensie, że dla uC2([0, 1], R) zachodzi Wracając do równości mamy Całki aproksymujemy na ogół numerycznie. Wygodnie jest jednak wprowadzić oznaczenie Teraz

20. Uwzgledniając warunki brzegowe mamy ponadto Można wszystkie wyrażenia Fi+1/2 zapisać jednolicie następujaco: Schemat numeryczny ma więc postać układu równań liniowych gdzie

21. Macierz A jest zdefiniowana tak, że Podsumujmy: rozwiązujemy układ równań

22. Metoda zostanie opisana na przykładzie ewolucyjnego zagadnienia ciepła z uwzględnieniem konwekcji. W przybliżeniu liniowym rozchodzenie się ciepła (ściślej: energii wewnętrznej) w ośrodku izotropowym jest proporcjonalne do gradientu temperatury: gdzie q oznacza wektor strumienia ciepła [J·m-2·s-1], a k jest współczynnikiem przewdodnictwa cieplnego w rozważanym ośrodku [J·m-1·s-1·K-1]. W przypadku, gdy na układzie nie jest wykonywana praca, energii wewnętrzna jednostkowej objętości materiału jest równa gdzie cp jest ciepłem własciwym, a   gęstością ośrodka.

23. Różniczkowa forma bilansu energii może być zapisana następująco: a zatem mamy Jeżeli teraz przyjmiemy, że ciepło właściwe i gęstość są stałe, to otrzymujemu gdzie cieplnym wpółczynnikiem dyfuzji (może zależeć od położenia!).

24. Jeżeli dopuścimy jeszcze możliwość konwekcyjnego transportu ciepła z prędkością v, to poprzednie równanie należy zastąpić następującym: gdzie kropka (·) oznacza standardowy iloczyn skalarny w R3, tzn. Ponieważ dalej będziemy rozważać przypadek bez źródeł, więc pole prędkości v nie może być dowolne. Brak źródeł lub ujść (ang. sinks) oznacza, że v musi być bezwirowe, a zatem:

25. Dzięki warunkowi div v =0 możemy jeszcze przekształcić nieco główne równanie, gdyż Ostatecznie będziemy rozwiązywać poniższe równanie: Równanie powyższe obowiązuje w jakimś obszarze   R3, który odpowiada naszemu ośrodkowi. Analogiczne równanie otrzymujamy dla np. transportu substancji zamiast energii („ciepła”). Wtedy zamiast T będziemy mieli c stężenie (lub gęstość) substancji.

26. Całkowa postać prawa zachowania (bilansu) energii Równanie całkujemu po dowolnie wybranym podobszarze V   ośrodka: Skorzystamy dale z tw. Gaussa: oraz z możliwości wyciągnięcia pochodnej przed znak całki (bo całka jest względem położenia, a nie czasu).

27. Całkowa postać prawa zachowania (bilansu) energii Stosując wspomniane zależności otrzymujemy: gdzie SV oznacza powierzchnię ograniczajacę elemnt V, a  jest polem wektorów normalnych do powierzchni SV. Wektor normalny do zamknietej powierzchni w jakimś jej punkcie oznacza: wktor jest prostopadły do powierzchni; ma długość jeden; jest skierowany na zewnątrz. Podkreślmy jeszcze raz: na tym etapie objętość kontrolna V („objętość skończona”  stąd nazwa metody) jest dowlnym objętościowym fragmentem dziedziny   R3, w której spełnione są równania. Oczywiście dalej przy konstrukcji schematu obliczeniowego wybiera się te elementy w sposób nieprzypadkowy.

28. Równanie transportu (w R2) Rozważmy liniowe równanie transportu