Etapy dyskretyzacji przestrzennej MES

Dyskretyzację przestrzenną metodą elementów skończonych przeprowadza się najczęściej w dwóch etapach. W pierwszym uzyskuje się sformułowanie słabe na poziomie nieskończenie wymiarowym, wykorzystując zerowanie się iloczynu skalarnego residuum rozważanego równania z funkcjami testującymi⁵. Istotną rolę, zwłaszcza dla numerycznej analizy zbieżności i błędu, odgrywają prze-strzenie funkcyjne, w których znajdują się funkcje niewiadome i testujące [73, 143, 1].

W drugim etapie dokonuje się redukcji zagadnienia do problemu skończenie wymiarowego poprzez ograniczenie rozważanych funkcji niewiadomych i testu-jących do elementów przestrzeni V_h, będącej skończenie wymiarową przestrze-nią funkcji określonych dla zadanego podziału obszaru Ω na elementy skoń-czone. Ta ostatnia przestrzeń jest najczęściej, także dla wszystkich rozważań w niniejszej pracy, przestrzenią funkcji będących wielomianami dla pojedyn-czego elementu. Dla tak zdefiniowanej przestrzeni sformułowanie słabe na po-ziomie skończenie wymiarowym nazywane będzie sformułowaniem skończenie elementowym.

Ze względu na położenie nacisku na zagadnienia implementacji, przedsta-wione poniżej wyprowadzenia sformułowań MES operują od razu pojęciem

5Sformułowanie słabe jest rozumiane jako pojęcie ogólniejsze od sformułowania waria-cyjnego. Obejmuje także przypadki, gdy nie istnieje minimalizowany lub maksymalizowany funkcjonał związany ze sformułowaniem.

przestrzeni funkcyjnej skończenie wymiarowej. Warunki brzegowe naturalne wchodzą w skład sformułowań, natomiast szczegóły implementacji warunków zasadniczych (np. funkcja kary) pozostawione są do implementacji specyficz-nych problemów (w notacji pominięte są możliwe modyfikacje przestrzeni V_h wprowadzane w celu uwzględnienia warunków brzegowych Dirichleta).

2.4.3 Podział obszaru obliczeniowego na elementy skończone W celu zdefiniowania dyskretyzacji MES dla określonego problemu, obszar obliczeniowy Ω dzielony jest na elementy skończone Ω_e. Każdy element przed-stawia się jako obraz elementu odniesienia ˆΩ_e poprzez odpowiednią wzajemnie jednoznaczną, odwracalną transformację T_e:

T_e : ˆΩ_e → Ω_e (2.21)

Każdemu punktowi elementu odniesienia, parametryzowanego przez współ-rzędne ξ, przyporządkowuje się punkt elementu rzeczywistego:

Te: ξ → x = x(ξ) (2.22)

Najczęściej elementy siatki są obrazami kilku podstawowych elementów od-niesienia: trójkąta, prostokąta, czworościanu, sześcianu, pryzmy (graniasto-słupa o podstawie trójkątnej, rys. 2.6) [91, 42]. Elementy jednej siatki mogą być obrazami różnych elementów odniesienia, przy spełnieniu założeń niena-kładania się elementów na siebie i sumowania się elementów do całości obszaru obliczeniowego. Transformacje T_e, będąc w dużym stopniu arbitralne, mu-szą jednak pozwalać na efektywne całkowanie funkcji za pomocą kwadratur zdefiniowanych dla elementów odniesienia i różniczkowanie funkcji określonych dla elementów odniesienia (nakłada to pewne ograniczenia na kształt obszaru obliczeniowego dla sformułowania MES, który często jest tylko przybliżeniem rzeczywistego obszaru, w którym zachodzi modelowane zjawisko).

W nieadaptacyjnej MES, a także w pewnych odmianach adaptacyjnej MES, wykorzystywane są wyłącznie siatki regularne, zwane także siatkami zgodnymi (conforming). Dla takich siatek wierzchołki elementów mogą znajdować się tylko na końcach krawędzi innych elementów, a każdy bok elementu jest jed-nocześnie bokiem elementu sąsiadującego. Przykładem takiej siatka jest siatka elementów trójkątnych przedstawiona na rys. 3.5.

Innym typem siatek są siatki nieregularne, zwane także siatkami niezgod-nymi (non-conforming). Najczęściej powstają one w wyniku zastosowania, dla początkowej siatki regularnej albo uprzednio zaadaptowanej siatki regularnej

x y

z

1.0

0.0 1.0

−1.0

1.0

Rys. 2.6. Pryzmatyczny element odniesienia

lub nieregularnej, adaptacji typu h, dopuszczającej dzielenie elementów bez modyfikacji ich sąsiadów. Wariantem nieregularności przydatnym do obliczeń dopuszczających zagęszczanie i rozrzedzanie siatki jest jednonieregularność.

Polega ona na tym, że nie pozwala się na dalsze podziały elementu, dopóki wszystkie jego boki i krawędzie nie są jednocześnie bokami i krawędziami ele-mentów sąsiadujących. Przykładami siatek jednonieregularnych są siatki po-kazane na rys. 2.1 i 2.4.

2.4.4 Wyprowadzenie sformułowania słabego dla przypadków aproksymacji ciągłej i nieciągłej

Poniżej przedstawione są wyprowadzenia sformułowań skończenie elemento-wych dla dwóch różnych typów aproksymacji: klasycznej ciągłej dyskretyzacji MES i nieciągłej dyskretyzacji Galerkina w wersji zaproponowanej w [125].

Nieciągła aproksymacja Galerkina ma kilka cech istotnie odróżniających ją od podejścia klasycznego (por. zbiór artykułów w [56]). Najważniejszą zaletą opisywanej wersji aproksymacji nieciągłej jest jej zachowawczy charak-ter (dla zachowawczych procesów ciągłych), lokalnie na poziomie pojedynczego elementu. Sformułowania klasyczne potrafią odtworzyć zachowawczość lokalną tylko za pomocą odpowiednich procedur postprocessingu [92]. Zachowawczość dyskretyzacji jest szczególnie istotna w modelowaniu zagadnień przepływów, zwłaszcza przy uwzględnieniu reakcji chemicznych, kiedy oscylacje rozwiąza-nia prowadzą do zachowań niefizycznych. Cechami negatywnymi prezentowa-nej wersji dyskretyzacji nieciągłej są: brak stabilności dla stopni aproksymacji niższych od dwóch i nieoptymalna zbieżność dla zagadnień eliptycznych w nor-mie L₂ [152].

W celu uzyskania sformułowania słabego równanie (2.16) jest mnożone przez funkcję testującą w, a następnie całkowane indywidualnie po obszarze każdego elementu. Po zastosowaniu wzorów na całkowanie przez części dla obszarów wielowymiarowych, otrzymuje się: względem brzegu elementu Γ_e.

Dla klasycznej ciągłej dyskretyzacji wysumowanie całek elementowych w całym obszarze powoduje zniesienie się wszystkich wyrażeń na brzegu mię-dzy elementami. Na brzegach obszaru, dla których zadane zostały warunki Neumanna i Robina, funkcje z warunków brzegowych podstawia się w miejsce odpowiednich wyrażeń. Warunki Dirichleta uwzględnia się poprzez modyfi-kacje sformułowania lub modyfimodyfi-kacje przestrzeni V_h⁶ [184, 91, 42]. Po kilku

6Modyfikacje te nie są analizowane w niniejszej pracy (choć model implementacji

uwzględ-prostych transformacjach uzyskuje się ostateczne sformułowanie MES dla

Inaczej jest w przypadku rozwiązań, dla których nie wymaga się ciągłości na brzegu międzyelementowym. Dla każdej pary sąsiadujących elementów Ω_e, Ω_f bok Γ_ef jest wspólną częścią ich brzegów, Γ_ef = Γ_e∩ Γ_f. Całość tak zdefi-niowanego brzegu międzyelementowego oznacza się poprzez Γ_int, Γ_int=^SΓ_ef. Dla każdego brzegu Γ_ef wybiera się (na podstawie arbitralnej konwencji) jeden wektor normalny jednostkowy n. Na brzegu zewnętrznym definiuje się n = n_e. Wprowadza się operatory skoku [ ] i uśrednienia <> dla funkcji określonych na Γ_ef:

[v] = v^L− v^R = v|_Γe∩Γ_ef − v|_Γf∩Γef (2.25)

< v >= 0.5 ∗ (v^L+ v^R) (2.26) Wysumowanie całek elementowych i zastosowanie powyższych operatorów daje w wyniku: Dla wyrazów dyfuzji na brzegu Γ_int stosuje się podstawienie:

[ab] =< a > [b] + [a] < b > (2.28)

nia możliwość ich stosowania).

oraz dodaje człon, zerujący się dla ciągłych funkcji u:

< A^ijnⁱw,j> [u]+ < w > [A^ijnⁱu,j] (2.29) Warunki brzegowe Neumanna i Robina stosuje się w analogiczny sposób jak w przypadku ciągłym, natomiast warunki Dirichleta wymusza się w sposób słaby, poprzez dodanie do obu stron równania (2.27) wyrazów:

Z

Ze względu na nieciągłość funkcji u i w istnieje możliwość modyfikacji wyra-zów konwekcyjnych (pierwszego rzędu) nie zmieniającej rozwiązań ciągłych. W niniejszym wyprowadzeniu przyjmuje się najprostszy, gwarantujący stabilność schemat dyskretyzacji „pod prąd” (upwind). Oznacza ona, że na brzegu mię-dzyelementowym używa się w wyrazach konwekcji wartości ¯u obliczanych po stronie wpływu. Strumienie qⁱ, które w ogólnym przypadku mogą być zarówno nieliniowe, jak i nieciągłe, są również dyskretyzowane pod prąd. Zastosowanie wymienionych wyżej kroków prowadzi do sformułowania:

Z

2.4.5 Ostateczna postać sformułowania skończenie