Symulacje komputerowe- dlaczego istnieją różne podejścia do wykonywania symulacji?

Pamiętam, jak zaczynałem wchodzić w świat symulacji komputerowych i zadawałem sobie pytanie jaka jest różnica między różnymi podejściami do wykonywania symulacji? Dlaczego jest coś takiego jak MES (Metoda Elementów Skończonych) oraz MOS (Metoda Objętości Skończonych)? Jaka jest różnica pomiędzy tymi dwoma podejściami? Przecież w obu przypadkach dzieli się swoją domenę obliczeniową na elementy siatki oraz wykonuję się obliczenia. Dlaczego w jednym przypadku ten pojedynczy sześcian nazywa się elementem skończonym, a w drugim przypadku nazywa się go objętością skończoną? Przecież one wyglądają tak samo. Dodatkowo za pomocą jednego podejścia i drugiego podejścia dostane wynik końcowy- dlaczego więc jest między nimi jakaś różnica? Jeżeli sam zadawałeś sobie takie pytania, to postaram się na nie odpowiedzieć w poniższym artykule.

Na samym początku należy zadać sobie pytanie: dlaczego w ogólne istnieją takie metody? Po co ktoś je wymyślił?

W symulacjach numerycznych, zawsze zależy użytkownikowi na rozwiązaniu pewnych równań fizyki, przykładami takich równań są równania Naviera- Stokesa, przenoszenia ciepła itp. Cały problem polega na tym, że tych równań nie da się rozwiązać dokładnie w postaci analitycznej dla rzeczywistych geometrii.

Dlatego też stosuje się metody dyskretyzacji, czyli dzieli się cały ciągły obszar i równania na małe fragmenty, które komputer potrafi już policzyć. Dwoma przykładami, które są bardzo często stosowane to właśnie wspomniane we wstępie metody: MES oraz MOS. Cała różnica pomiędzy tymi metodami polega na podejściu do problemu.

W MESie ideą jest spróbowanie opisania pola za pomocą „ładnych” funkcji i sprawdzenia jak dobrze one odwzorowują fizykę. Te „ładne funkcje” nazywają się funkcjami kształtu, które przybliżają wartości zmiennych (np. prędkość) za pomocą równań matematycznych. Filozofią MOSa jest podzielenie obszaru na małe kostki i policzenie przepływu masy/energii przez ich ściany tak, żeby nic się nie zgubiło. Czyli MOS bazuje na bilansie fizycznym, czyli zachowują dokładną równowagę masy, pędu i energii w każdej objętości siatki (małej kostce).

Głównymi równaniami MESu są równania różniczkowe, które są przekształcane w słabe formy (integralne). Dla MOSu głównymi równaniami są równania zachowania (konserwatywne), które są integrowane po objętości. Na tym etapie należy to zaznaczyć, że MES nie zawsze dokładnie zachowuje bilans masy/energii. Wszystko zależy od przybliżeń.

Czyli jak to wszystko wygląda w praktyce?

Proszę sobie wyobrazić, że należy zasymulować przepływ powietrza przez rurę. Sposób działania obu metod wygląda następująco:

MOS

Dzielisz rurę na wiele małych objętości (np. sześcianów).
Dla każdej z nich mówisz:

„Powietrze, które wpłynie przez jedną ścianę, musi wypłynąć przez inną – masa się nie może zgubić.”

Liczysz strumienie przez każdą ścianę (masy, pędu, energii).
Wynik: zachowujesz równowagę – wszystko się zgadza fizycznie.

MES

Dzielisz rurę na elementy (np. trójkąty albo czworościany).
W każdym elemencie przybliżasz prędkość i ciśnienie równaniami (np. liniowo, kwadratowo).
Łączysz elementy przez tzw. funkcje kształtu i rozwiązujesz globalny układ równań.

Czyli można powiedzieć, że MOS myśli fizycznie, czyli pilnuje, żeby bilanse się zamykały, natomiast MES myśli matematycznie, czyli dopasowuje pole zmiennych do równań różniczkowych.

Dobrze, sama idea działania obu metod została pokazana powyżej, ale… Przecież dla obu przypadków dzieli się domenę obliczeniową na elementy siatki. Na ekranie komputera jest „chmura małych bryłek” i dlaczego czasami mówi się o elementach skończonych, a czasami o objętościach kontrolnych?

Żeby to dobrze zrozumieć, trzeba zrozumieć filozofię tego, co ta siatka reprezentuje matematycznie i fizycznie.

Zarówno w MES, jak i MOS dzielimy analizowany obszar (np. rurę, profil, zbiornik) na małe fragmenty, żeby komputer mógł obliczyć, co się dzieje lokalnie.

Ale cel i znaczenie tych fragmentów są inne. W MESie siatka składa się z elementów geometrycznych (np. trójkątów, czworościanów, sześcianów), które służą do aproksymacji pola zmiennych (np. temperatury, przemieszczenia, prędkości), czyli w każdym elemencie zakłada się, że dana wielkość (np. temperatura T) zmienia się według znanej funkcji (funkcji kształtu). Funkcja te są ciągłe między elementami, czyli sąsiednie elementy „dogadują się” na swoich granicach. Równania fizyczne przekształca się w postać matematyczną (słabą formę) i rozwiązuje się globalny układ równań. Czyli element to fragment, w którym opisujemy, jak zachowuje się pole wielkości, niekoniecznie jak przepływa przez jego granice. Można powiedzieć, że element w MESie to miejsce, w którym przybliżamy zachowanie pola – jego matematyczny kształt w przestrzeni.

W MOSie mamy do czynienia z objętościami kontrolnymi, które reprezentują realne fragmenty przestrzeni, w których zachodzi bilans masy, energii i pędu. Czyli dla każdej objętości zakładamy, ile masy/energii wchodzi i wychodzi przez jej ścianki. Zapisywane jest równanie bilansowe:

wejście-wyjście+źródło=zmiana wewnątrz

Dzięki termu metoda dokładnie zachowuje bilanse fizyczne, czyli nic się nie „zgubi”. Objętość kontrolna ma więc znaczenie fizyczne – to jak mały „zbiorniczek”, w którym mierzymy co wpływa, co wypływa i ile zostało.

Żeby lepiej zrozumieć różnicę między tymi podejściami to należy wyobrazić sobie rzekę i zadać sobie pytanie- jak płynie woda? Dla MESu, dzieli się rzekę na odcinki i mówi się: „W tym kawałku prędkość zmienia się w przybliżeniu liniowo. Sprawdźmy, jak dobrze to pasuje do równań ruchu.” Czyli patrzy się na kształt profilu prędkości, niekoniecznie na to, ile wody przepływa przez przekrój rzeki. Natomiast dla MOSu, rzekę dzieli się na „zbiorniczki” i mówi się: „Ile wody wpłynęło do tego zbiorniczka? Ile wypłynęło? Ile przybyło w środku? Bilans musi się zgadzać.”.

Czyli patrzy się na przepływ przez granice, a nie na kształt rozkładu wewnątrz.

Podsumowując oba podejścia tworzą „siatkę”, ale:

MES myśli: „Jak najlepiej przybliżyć pole?”
MOS myśli: „Jak policzyć, co wpływa i wypływa, żeby zgadzał się bilans?”

Dobrze teraz znana jest różnica między podejściem MESem, a MOSem, ale jak to wpływa na wyniki symulacji? Oczywiście zrozumiałe jest już na tym etapie, że w MESie nie musi zostać spełniony bilans, ale przecież ta metoda też dostarcza wyników np. dla temperatury. Dlaczego to podejście jest gorsze/lepsze?

Zanim odpowiem na to pytanie, chce zaznaczyć, że to co znajduje się poniżej jest główną przyczyną, dlaczego CFD prawie zawsze używa MOS, a nie MES.

W obu metodach dostaje się ładny wykres pola, ale przez to, że w MESie nie musi zostać zachowany bilans, to nie musi zostać tam zachowana fizyka przepływu. Bilans masy lub energii może nie być dokładny lokalnie, bo jest spełniony tylko w „średnim” sensie (po całej domenie lub po zbiorze elementów). Oznacza to, że temperatura w jednym miejscu może być lekko zawyżona, a w innym zaniżona, ale średnio równanie przewodzenia się zgadza. Takie podejście może się mścić, bo nawet drobne naruszenie bilansu masy powoduje niefizyczne przyspieszenia, wiry albo „uciekającą” energię. Powoduje to, że odwzorowanie, np. turbulencji może być trudna do odwzorowania bez dodatkowych zabiegów. Dodatkowo stabilność obliczeń np. przy wymianie ciepłą jest cięższa do uzyskania, jeżeli nie kontroluje się bilansu energii pomiędzy poszczególnymi elementami.
 Przez to, że w MOSie stosuje się bilans energii, to bilans całkowity domeny jest dokładany. Żadna masa, energia czy pęd nie „znika” w obliczeniach. Jeżeli coś się grzeje, to tylko dlatego, że wpłynęła energia, albo występuje źródło ciepła. Jeżeli ciśnienie maleje, to znaczy, że rzeczywiście energia została rozproszona. Powoduje to, że obliczenia są stabilne, bo są oparte na strumieniach.

A dlaczego jest to tak istotne w symulacjach CFD, żeby masa, pęd oraz energia bilansowały się dokładnie? To dlatego, że równania Naviera-Stokesa są równaniami ciągłości. Nawet mikroskopijny brak bilansu może się „rozlać” po domenie — np. ciśnienie może rosnąć bez powodu albo pojawią się sztuczne źródła energii kinetycznej.

Teraz to co chciałem pokazać to przykład przewodzenia ciepła w pręcie, który będzie rozważany tylko w 1 wymiarze. Takie podejście pozwoli zobaczyć, gdzie są obliczane wartości w obu przypadkach.

Założenia przykładu:

Metalowy pręt ma długość 1 m.
Na jednym końcu panuje temperatura 100°C, a na drugim końcu panuje temperatura 0°C.
Przewodnictwo cieplne k = 10 ^W/_(m∙K)
Brak źródeł cieplnych, stan ustalony

Analityczne Rozwiązanie to prosta linia, opisana równaniem:

T(x)=100-100x

Czyli temperatura maleje liniowo od 100 °C do 0 °C.

Pierwszym krokiem do zastosowania obu metod jest posiatkowanie pręta. W tym przypadku pręt został podzielony na 4 fragmenty, czyli ma 5 węzłów:

Numer	Pozycja x[m]
1	0,0
2	0,25
3	0,5
4	0,75
5	1,0

MES

Pierwszym etapem jest stworzenie 4 elementów: (1–2), (2–3), (3–4), (4–5)
Zakładam, że temperatura zmienia się liniowo w każdym elemencie.
Węzły to punkty, w których rozwiązanie jest szukane
W środku elementu temperatura jest obliczona z funkcji kształtu, czyli interpolowana z węzłów
Dla równania przewodzenia

– ^d/_dx (k ^dT/_dx) = 0

MES tworzy tzw. „sztywną macierz”
W stanie ustalonym równania są liniowe i po rozwiązaniu otrzymuje się

Węzeł	MES-temperatura [0 °C]
1	100
2	75
3	50
4	25
5	0

Sprawdzenie bilansu energii po lewej stronie

q_L= -k ^{T₂ – T₁} / _∆x = -10 ^. ^{75 – 100} / _0,25 = 10000 W/m²

Sprawdzenie bilansu energii po prawej stronie

q_p = -k ^{T5 – T₄} / _∆x = -10 ^. ^{0 – 0,25} / _0,25 = 1000 W/m²

Sprawdzenie bilansu energii po prawej stronie
Bilans się zgadza, bo to trywialny przypadek liniowy.
Ale gdyby:

przewodność k się zmieniała,
elementy były różnej długości,
albo dodano źródła ciepła,

to lokalny bilans energii w każdym elemencie nie byłby zachowany dokładnie — byłby spełniony tylko w sensie „średnim”.

MOS

Pręt jest dzielony na 4 objętości kontrolne, każda wokół środka przedziału.
Każda objętość ma ścianki (lewa/prawa), przez które płynie ciepło.
Wartość temperatury przypisana jest do środka każdej objętości
Równania zapisujemy w postaci bilansu energii:

q_wejście – q_wyjście = 0

Dla i-tej objętości:

k ^{T_{i – 1} -T_i} / _∆x = k ^{T_i – T_i+1} / ∆

Rozwiązując układ równań (dla 4 komórek + warunków brzegowych):

Środek komórki x [m]	MOS-temperatura [0 °C]
0,125	87,5
0,375	62,6;
0,625	37,5
0,875	12,5

Wynik też liniowy, zgodny z fizyką, ale teraz uwaga:

Każda objętość ma lokalnie dokładny bilans:

strumień przez lewą ściankę=strumień przez prawą ściankę

To jest gwarantowane konstrukcją metody. Nic nie można zgubić

Na tym przykładzie jasno widać, że MES rozwiązuje równania w węzłach siatki, a wartości między węzłami są interpolowane z funkcji kształtu. Natomiast MOS rozwiązuje obliczenia w środku objętości kontrolnej, a strumienie mas i energii liczone są na ścianach między komórkami.

Po przeczytaniu tego wszystkiego można sobie teraz zadać pytanie „A to, dlaczego w symulacjach wytrzymałościowych stosuje się MES, a nie MOS?”. Na pierwszy rzut oka przecież „bilans energii” czy „równowaga sił” wydają się tak samo ważne w mechanice ciała stałego, jak bilans masy czy pędu w CFD. Dlaczego, więc w analizach wytrzymałościowych używa się MES, a nie MOS?

Kluczowa różnica: co opisuje równanie równowagi w mechanice ciała stałego?

W mechanice płynów podstawą jest lokalne prawo zachowania (bilans), czyli:

zmiana w objętości = strumienie przez granice + źródła

Wszystko dzieje się dynamicznie w przestrzeni i czasie, a mały lokalny błąd powoduje naruszenie zachowania masy lub energii.

W mechanice ciała stałego (dla małych odkształceń) podstawowe równanie ma inną postać:

∇∙σ+f=0

Czyli równowagę sił w stanie ustalonym.

To nie jest równanie zachowania w sensie „przepływu przez granice” (jak w CFD), lecz równanie równowagi — czyli warunek, że suma sił (wewnętrznych i zewnętrznych) jest równa zeru.

Nie ma tu fizycznego przepływu masy, energii ani pędu przez granice elementów. MES nie gwarantuje, że w każdej pojedynczej komórce równowaga sił jest lokalnie spełniona,
ale gwarantuje, że równowaga globalna całego układu jest dokładna — i to wystarcza.

Kolejnym pytaniem, które może się pojawić jest:” A co z energią odkształcenia?”.
W MES energia potencjalna odkształcenia jest automatycznie zachowana w sensie globalnym, bo cała metoda opiera się na zasadzie minimalizacji energii potencjalnej:

δΠ=0

Czyli system znajduje stan, w którym energia potencjalna całkowita (odkształcenia + praca sił zewnętrznych) jest stacjonarna. To daje poprawne rozwiązanie globalne, bez potrzeby zachowania energii w każdej komórce z osobna.

Teoretycznie można korzystać z metody MOS do symulacji wytrzymałościowych, ale:

MOS jest trudniejszy do zastosowania w nieliniowej sprężystości (bo wymaga uśredniania gradientów odkształcenia),
MES łatwiej radzi sobie z dużymi odkształceniami, kontaktami, elementami powłokowymi,
MES ma naturalną postać macierzową, pasującą do klasycznych solverów statycznych.

Po tym wszystkim może pojawić się kolejne pytanie „A co w symulacjach wytrzymałościowych, ale zmiennych w czasie? Czy MOS nie byłby odpowiednim podejściem np. w crash testach?”.

Tutaj nadal wykorzystuje się metodę MES. Dlaczego? Dlatego że, w MOS interesuje użytkownika ruch płynu, czyli pola prędkości, ciśnienia czy temperatury. Płyn zmienia kształt dowolnie- nie ma pojęcia „sztywnego elementu”. Dlatego przestrzeń jest dzielona na objętości kontrolne. Natomiast w analizie mechanicznej, użytkownik ma do czynienia z ciałem o określonym kształcie, które się odkształca, ale punkty materiału pozostają „połączone”. To co jest interesujące to naprężenia i odkształcenia.

Oczywiście bilans energii (kinetycznej+ odkształcenia+ tłumienia) oraz bilans pędu w crash testach jest też ważny, ale w trochę inny sposób. Tam wszystko wynika z praw ruchu bryły (drugiej zasady Newtona), a nie z przepływu przez objętości kontrolne. Oczywiście w crash testach ma się do czynienia z ekstremalnymi przemieszczeniami, zderzeniami elementów bądź z lokalnymi wgnieceniami. W takich przypadkach w MESie korzysta się z siatki związanej z materiałem (siatka Lagrange’owska), czyli takiej siatki, w której węzły poruszają się razem z materiałem. Takie podejście umożliwia dokładne śledzenie deformacji, naprężeń i kontaktu między częściami. W MOSie siatka jest związana z przestrzenią (Eulerwoska) i w takim przypadku materiał przepływa przez komórki. Oznacza to, że nie da się naturalnie opisać zjawiska, w którym materiał się zgina, zgniata i zachowuje spójność. Gdyby próbowało się użyć MOS w crash teście to konieczne byłoby opisać deformujące się ciało jako przepływającą ciecz. Dodatkowo należałoby śledzić granicę między powietrzem a materiałem i mieć nieliniowe modele lepko-sprężyste o bardzo dużej lepkości. To wszystko byłoby ekstremalnie skomplikowane i mało praktyczne.

Podsumowując, MOS zachowuje lokalny bilans masy, pędu i energii, dlatego jest podstawą symulacji przepływów w CFD, gdzie kluczowe jest ścisłe spełnienie równań zachowania w każdej objętości kontrolnej. Natomiast MES opiera się na tzw. słabej postaci równań i śledzi deformacje materiału w układzie Lagrange’a, dzięki czemu dokładnie opisuje naprężenia, kontakt i pękanie, ale nie gwarantuje lokalnego bilansu. W efekcie MOS dominuje w analizach płynów, a MES w analizach mechanicznych, ponieważ każda z metod lepiej odpowiada fizyce zjawisk, które opisuje.