Implementacja Metody Elementów Skończonych
w środowisku Matlab

Marcin Maślanka
student V roku, spec. Automatyka i Metrologia, IMiR, AGH

Streszczenie - Metoda Elementów Skończonych (MES) pomimo niekwestionowanej, ogromnej przydatności w zastosowaniu do licznych, złożonych zadań obliczeniowych jest bardzo uciążliwa w użyciu ze względu na konieczność stosowania specjalistycznego, drogiego oprogramowania. Celem opisywanej pracy była próba przeniesienia algorytmów MES do ogólnie dostępnego w kręgach inżynierskich środowiska Matlab. Pakiet wykorzystano w zastosowaniu do analizy wybranych, prostych zagadnień, w celu zbadania jego przydatności jako platformy do zaawansowanych obliczeń z wykorzystaniem MES. W opracowaniu przedstawiono stworzony przez autora zestaw procedur, uruchomionych w środowisku Matlab, które wraz z interfejsem użytkownika umożliwiają, przy zastosowaniu tytułowej metody, na przeprowadzanie obliczeń statycznych dowolnie określanych przez użytkownika geometrycznych struktur płaskich.

  SPIS TREŚCI

1. Wstęp
2. Prezentacja programu EasyMesLab v1.0
    2.1 Zestawienie dostępnych funkcji programu
    2.2 Sposób uruchomienia i opis działania programu
    2.3 Wybrany przykład zastosowania
3. Podsumowanie
Literatura

    Metoda Elementów Skończonych wykorzystywana jest do analizy różnego rodzaju ośrodków ciągłych, przy czym, owe ciągłe kontinuum dyskretyzuje się na elementy o prostych kształtach połączonych ze sobą w węzłach. Przemieszczenia węzłów opisywane są za pomocą układu równań różniczkowych (w zagadnieniach dynamicznych) lub algebraicznych (w przypadku warunków statycznych). Wybrane wielkości węzłowe opisuje się dla wnętrza elementów przy pomocy funkcji interpolujących, dobranych tak aby zachować ciągłość na granicach poszczególnych elementów.
Postępowanie numeryczne sprowadza się zatem do rozwiązania odpowiedniego układu równań utworzonych dla przyjętej sieci elementów skończonych, a następnie wykorzystaniu funkcji interpolujących dla uzyskania poszukiwanych parametrów wewnątrz elementów.

Pomimo niekwestionowanej, ogromnej przydatności Metody Elementów Skończonych w zastosowaniu do licznych, złożonych zadań obliczeniowych, jej stosowanie wiąże się nadal nierozłącznie z koniecznością stosowania specjalistycznego, bardzo drogiego oprogramowania.
Celem badań było przeniesienie algorytmów MES do ogólnie dostępnego w kręgach inżynierskich środowiska programu Matlab. Badając przydatność Matlaba do zaawansowanych obliczeń tą metodą opracowałem program EasyMesLab (prezentowany na XXXVII Sesji Studenckich Kół Naukowych AGH 11 maja 2000r.) wykorzystujący MES do analizy dowolnie definiowanych przez użytkownika struktur płaskich, budowanych z elementów skończonych prętowego i trójkątnego. Wyniki obliczeń jakie uzyskiwałem w opracowanym przeze mnie programie porównywałem z wynikami podawanymi w literaturze oraz z wynikami uzyskiwanymi w profesjonalnym programie MSC/Nastran (obiektem badań była poddana różnym obciążeniom płaska tarcza z otworem, podzielona na ponad 200 trójkątnych elementów skończonych).
Analiza wyników wskazuje na poprawność stosowanych algorytmów.

W opracowaniu tym pominięto opis teoretyczny założeń MES. Zainteresowanym poleca się publikację O. Zienkiewicza [1].

  W obecnej wersji programu EasyMesLab zaimplementowano zbiór opracowanych przez autora procedur narzędziowych umożliwiających edycję struktur geometrycznych płaskich złożonych z elementów skończonych prętowego i trójkątnego oraz statyczną analizę przemieszczeń, odkształceń i naprężeń przy założeniu obciążenia przykładanego w węzłach.

    Algorytmy obliczeniowe programu zostały gruntownie przetestowane na wielu przykładach, dla których wyniki obliczeń uzyskane były inną drogą. Dostępność kodu źródłowego umożliwia rozbudowę lub modyfikację programu a także podnosi jego wartość dydaktyczną.

2.2 Sposób uruchomienia i opis działania programu

   Opis działania programu dokonany zostanie na przykładzie. Utwórzmy więc nową strukturę poleceniem Plik|Nowy. Na ekranie pojawia się obszar edycyjny geometrii struktury w płaszczyźnie X-Y. Obowiązującą jednostką geometrii jest metr. Wszystkie elementy wystroju okna graficznego można odpowiednio dobrać korzystając z opcji menu Dostosuj. Przyciski z lewej strony okna umożliwiają edycję struktury. Po wciśnięciu przycisku Edycja nowego elementu wyświetla się okno wprowadzania danych jak na rys.1 (18kB). Dla każdego z tworzonych elementów skończonych (składowych budowanej struktury) można w tym oknie utworzyć typ elementu określający rodzaj materiału i geometrii elementu. Po utworzeniu odpowiedniego typu (załóżmy w tym przykładzie, że tworzoną strukturą jest płaska kratownica, a wszystkie elementy są jednego typu - prętowe o przekroju 12 [mm²]) podaje się w odpowiednich polach współrzędne węzłowe kolejnych elementów skończonych podając dla elementu prętowego dwa węzły, dla elementu trójkątnego trzy. Numeracji węzłów dokonuje sam użytkownik programu - numeracja elementów jest automatyczna i zgodna z kolejnością wprowadzania elementów. Gdy kolejny wprowadzany element ma mieć węzeł wspólny z węzłem elementu już istniejącego należy w odpowiednim miejscu podać jedynie numer już istniejącego węzła (nie podajemy dwa razy współrzędnych jednego węzła). Przykładowo wpiszmy: w odpowiednich polach:

nr_węzła:1   x [m]= 0.1    y [m]=0.1
nr_węzła:2   x [m]= 0.1    y [m]=0.1

Po wciśnięciu przycisku Dodaj utworzony został element prętowy o numerze 1 i węzłach w punktach o podanych współrzędnych. Łączymy go z kolejny elementem:

nr_węzła:3   x [m]= 1.1    y [m]=0.3
nr_węzła:2   (nie podajemy kolejny raz współrzędnych węzła 2)

a ten z kolei z elementem:

nr_węzła:4   x [m]= 1.4    y [m]=0.1
nr_węzła:3   (uwaga j.w)

Wygląd obszaru edycyjnego po dokonaniu tych operacji pokazano na rys.2 (67kB).

Po dokonaniu edycji geometrii przechodzimy do edycji utwierdzeń i obciążeń. Opcje Edycja utwierdzeń oraz Edycja obciążeń pozwalają zarówno na przyłożenie obciążenia lub utwierdzenia do wybranego węzła jak i na uwolnienie tego węzła od przyłożonej wcześniej siły lub utwierdzenia. Odpowiednio utwierdzony i obciążony model poddaje się analizie - przycisk Analiza statyczna struktury.

W celu przedstawienia przykładowych wyników analizy otworzony został przykładowy plik z katalogu Example/ Wszebor/wszebor.wk1. Wynik analizy przedstawiono na poniższych rysunkach:

Rys.3 Oznaczone przemieszczenia węzłów dla rozpatrywanej struktury

Rys.4 Sposób prezentacji wybranych wartości
przemieszczeń, odkształceń i naprężeń

    Podając numer węzła w poszczególnych polach wprowadzania danych otrzymuje się wartości przemieszczeń, odkształceń i naprężeń. Skala wektorów przemieszczeń wyznacza ich wielkość na obszarze edycyjnym. Przemieszczenia są zazwyczaj znikomo małe w stosunku do geometrii modelu zatem aby móc je przedstawić w obszarze geometrii struktury wprowadzenie odpowiedniego współczynnika skali jest konieczne. Skalowanie przebiega względem największego wektora przemieszczeń - w tym przypadku maksymalny wektor przemieszczenia przyjmuje wartość 0.1 wg skali obszaru edycji, a pozostałe wektory przemieszczeń są proporcjonalnie do tej wartości przeliczone.

    Dla struktur zbudowanych z elementów trójkątnych wybiera się rodzaj analizy poprzez wybranie jednej z opcji spośród opcji wyszczególnionych żółtą czcionką. Możliwe jest ponadto pokrycie elementów kolorem odpowiadającym wartości obliczonych wielkości dla danego elementu wg przyjętej palety barw. Przykładowo dla struktury zapisanej w pliku Example/Tarcza/tarcza30.wk1 rozkład naprężeń zredukowanych (naprężenia zredukowane wyznaczane są w programie wg hipotezy Hubera) przedstawiono na rys.5 (115kB).

    Zmiany stosowanej palety barw dokonuje się w menu Dostosuj | Paleta kolorów wizualizacji. Po zmianie palety należy odświeżyć obszar edycyjny poprzez ponowne wybranie wizualizacji elementem popupmenu - Opcje wizualizacji wyników. Skalowanie barw odbywa się wg największej wartości wielkości, której rozkład się wizualizuje.

Na rys.6 (33kB) przedstawiono sposób w jaki zostają prezentowane wyniki obliczeń.

Uwaga: W przypadku nie importowania przez wykorzystywaną wersję programu Matlab użytych w programie EasyMesLab czcionek, może okazać się, że oznaczenia naprężeń, odkształceń lub symbole przyjęte przez autora programu na oznaczenie podpór i obciążenia przyłożonego do węzłów struktury będą inne.

2.3 Wybrany przykład zastosowania

    W pliku Example/Tarcza/tarcza205.wk1 zapisano geometrię i sposób utwierdzenia/obciążenia tarczy z otworem podzielonej na 205 elementów skończonych, dla której wyniki obliczeń uzyskane w programie EasyMesLab porównywane były z wynikami uzyskanymi w profesjonalnym programie Nastran.
Utwierdzenie: na lewostronnych narożach tarczy.
Obciążenie: Fx=-10 [N], Fy=-10[N] przyłożone do prawego, górnego naroża tarczy.

    Teoretyczny przykład obciążonej tarczy, który posłużył do sprawdzenia poprawności algorytmów EasyMesLab zaczerpnięto z pracy prof. Zienkiewicza [1] (str. 74), gdzie podano również literaturę przedstawiającą analitycznie uzyskane wyniki dla tego rodzaju przykładu.

    Wyniki obliczeń uzyskane w programie Nastran zapisano do pliku tarcza1 (119kB)
Numery węzłów i elementów przyjętych dla tego przykładu w Nastranie przedstawiają pliki rys.7 (77kB) (zbliżenie środka tarczy) oraz rys.8 (70kB)
Plik tarcza2 (119kB) prezentuje wyniki obliczeń programu Nastran dla tej samej struktury obciążonej symetrycznie czterema siłami przyłożonymi do krawędzi otworu. Wartości sił: F= 100 [N]. Wszystkie pliki przykładów udostępniane są przez autora wraz z plikami źródłowymi programu EasyMesLab.

    Po uwzględnieniu, że numeracja elementów i węzłów przyjęta w programie Nastran różni się od numeracji przyjętej w programie EasyMesLab analiza porównawcza wyników wykazuje ich dużą zgodność stanowiąc podstawę do stwierdzenia o poprawności stosowanych w programie EasyMesLab algorytmów. 
Na rys.9 zamieszczono przykładowy wynik wizualizacji rozkładu jednej z wyznaczonych wielkości na poszczególne elementy rozpatrywanej w przykładzie tarczy.

Rys.9 Przykładowy wynik wizualizacji wybranej wielkości dla badanej tarczy

    Dla prostego przykładu dwóch elementów trójkątnych wyniki całkowicie zgodne z wynikami otrzymanymi analitycznie przedstawiono w katalogu Example/Orlos364. Dla struktur złożonych z elementów prętowych zagadnienie znacznie się uprasza. Poprawność wyników otrzymywanych w programie potwierdzona została wielokrotnie. Wybrane, proste przykłady zamieszczono w katalogu Example/Orlos264 oraz Example/Wszebor.

    Zaawansowana implementacja Metody Elementów Skończonych w Matlabie byłaby z pewnością bardzo przydatna. Zwraca uwagę łatwość posługiwania się wbudowanym w to środowisko językiem programowania, dostępność zaawansowanych metod numerycznych, sposób funkcjonowania zmiennych w przestrzeni roboczej oraz łatwość tworzenia zgodnego ze standardem Windows 9x graficznego interfejsu użytkownika. Jest to zatem niezwykle sprzyjające środowisko do szybkiego tworzenia dowolnie zaawansowanych obliczeniowo aplikacji o strukturze otwartej. Po zaimplementowaniu szeregu podstawowych procedur można byłoby stosować pakiet do budowy aplikacji służących do rozwiązywania pewnych niestandardowych zagadnień przy wykorzystaniu MES.

Kierunek dalszych badań: kolejnym krokiem implementacji mogłoby być wprowadzenie przestrzennych elementów skończonych, a zatem analizy struktur trójwymiarowych, jak również przejęcie przez program automatycznego (z możliwością korekcji) dzielenia wprowadzonej do programu geometrii na elementy skończone. Możliwa jest również adaptacja programu do obliczeń Metodą Elementów Skończonych innych zagadnień. Wszystkich zainteresowanych rozwojem programu proszę o kontakt: masmar1pl@yahoo.com

[1] O. Zienkiewicz Metoda Elementów Skończonych. WNT Warszawa 1977