2. Medyczne i techniczne podstawy pracy
2.4. Metoda elementów skończonych
Analiza obiektów o złożonej budowie, własnościach materiałowych czy działającym układzie sił sprowadza się do stworzenia modelu matematycznego opisanego równaniami różniczkowymi. Uzyskanie analitycznego rozwiązania takiego problemu jest nierzadko karkołomnym lub niewykonalnym zadaniem. W tychże przypadkach stosuje się metody numeryczne w szczególności metodę elementów skończonych (ang.
Finite Element Method, FEM), której główną ideą [105][106] jest aproksymacja skomplikowanego opisu matematycznego zestawem równań algebraicznych. W zagadnieniach inżynierskich takie przybliżenie rozwiązania jest wystarczające.
Najprostszym przykładem działania metody elementów skończonych jest znalezienie długości okręgu o nieznanym promieniu R. Na takim okręgu (rys. 8a) można opisać (lub wpisać weń) wielokąt foremny, w którym każda ze ścian jest w istocie elementem skończonym (rys. 8b). W zależności od przyjętej metody tworzenia wspomnianego wielokąta uzyskujemy sumę długości jego boków jako wartość nieco mniejszą lub większą od rzeczywistej długości okręgu. Należy tu zauważyć, iż zwiększanie liczby elementów skończonych będzie zbliżać aproksymowane rozwiązanie do właściwego (rys. 8c), jednak czas potrzebny na obliczenia może się wydłużać.
36
Rysunek 8. Schemat ideowy MES: a) okręg o nieznanym promieniu R, b) sześciokąt foremny wpisany w okręg, c) szesnastokąt foremny wpisany w okręg.
Historia metody elementów skończonych sięga lat czterdziestych XXw. , kiedy to Hrennikoff [107] i McHenry [108] zastosowali układ elementów liniowych jako reprezentację obiektu ciągłego. Jednak dopiero na przełomie lat pięćdziesiątych i sześćdziesiątych Argyris i Kelsey [109] oraz Turner [110] przedstawili pełne pryncypia owej metody. Samo sformułowanie „element skończony” zostało po raz pierwszy użyte przez Clough w 1960 roku [111]. Martin [112] a także Melosh [113] byli pionierami stosowania MES w zagadnieniach układów trójwymiarowych. Większość prac z tamtego okresu podejmowała tematykę obiektów o elastycznych właściwościach materiałowych, poddanych niewielkim, statycznym obciążeniom i przemieszczeniom.
Na lata sześćdziesiąte XX wieku przypada okres intensywnego rozwoju metody elementów skończonych. I tak, analiza zjawisk termalnych została wprowadzona przez Turner [114], materiałów nieliniowych przez Gallagher [115], własności lepkosprężystych oraz zagadnień polowych przez Zienkiewicza [107][108], mechaniki płynów przez Martin [118] i ostatecznie przewodnictwa cieplnego przez Wilson i Nickel [119]. Współcześnie, nowym polem zastosowania MES jest bioinżynieria, która wciąż stanowi wyzwanie dla metod numerycznych ze względu na złożoną geometrię i silnie nieliniowe właściwości materiałowe badanych obiektów.
W ogólności, metoda elementów skończonych składa się z kilku etapów [120][121].
Kroki te zostaną po krótce przedstawione poniżej.
1. Dyskretyzacja
Jest to proces, podczas którego obiekt ciągły zostaje podzielony na mniejsze regiony tj. elementy skończone (ang. finite element). Zbiór wszystkich elementów zwany jest siatką elementów skończonych (ang. finite element mesh).
37 Siatka ta może być jednorodna, wówczas wszystkie jej składniki są tych samych rozmiarów, lub niejednorodna, gdy elementy są różne. Elementy są ze sobą połączone w punktach zwanych węzłami (ang. nodes).
Istnieją różne typy elementów skończonych, z których każdy optymalny jest dla specyficznego zastosowania. Jeśli dany obiekt ma jeden wymiar znacznie większy od pozostałych, to kształt obiektu można przybliżyć elementem liniowym. W najprostszej wersji jest to odcinek, na którego obu końcach znajdują się węzły (rys. 9a). Ma on zastosowanie w takich zagadnieniach jak rozkład temperatury w pręcie czy deformacja belki pod wpływem osiowego rozciągania, a także, jako złożenie wielu liniowych elementów skończonych, w modelowaniu mostów, wież i budynków.
Podstawowym kształtem elementu dla zagadnień dwuwymiarowych jest trójkąt, choć niekiedy przekształca się go do odpowiedniego złożenia w kwadrat, prostokąt czy inny równoległobok (rys. 9b). W tym przypadku węzły znajdują się w narożnikach elementu. Ten typ elementu stosowany jest jako aproksymacja obiektów powierzchniowych, gdzie grubość jest znikoma w stosunku do dwóch pozostałych wymiarów, np. płachty metalowe, poszycie skrzydeł samolotu.
Problemy rozpatrywane w trzech wymiarach wymagają objętościowych elementów takich jak tetrahedrony czy heksahedrony (rys. 9c). Są to elementy najczęściej stosowane w analizie złożonych geometrii. Ponieważ nie upraszczają aproksymowanego obiektu, dają najwierniejsze rezultaty, lecz wymagają przy tym bardzo drobnej siatki elementów skończonych.
Dodatkowo, by poprawić dopasowanie elementów do obiektów o krawędziach w formie łuków, stosuje się elementy o krawędziach wyższych rzędów (rys. 10).
Typ ten charakteryzuje się dodatkowymi węzłami w krawędziach elementu skończonego.
38
Rysunek 9. Podstawowe kształty elementów skończonych pierwszego rzędu: a) element 1D, b) elementy 2D, c) elementy 3D [122]
Rysunek 10. Przykładowe elementu skończone drugiego rzędu: a) element 1D, b) elementy 2D, c) elementy 3D [122]
Powyżej przedstawione zostały bazowe typy elementów. Niemniej, istnieje cały szereg elementów przeznaczonych do modelowania specyficznych obiektów, jak choćby elementy osiowosymetryczne. Charakterystyka wszystkich ich rodzajów wykracza poza tematykę niniejszej pracy.
2. Równania elementu
Należy przypomnieć, że do opisu zachowania dowolnego obiektu, także pojedynczego elementu skończonego, stosuje się model matematyczny wiążący zarówno ogólne prawa fizyczne oraz równania konstytutywne. Te pierwsze, będące zazwyczaj równaniami różniczkowymi, opisują oddziaływania
39 niezależne od substancji rozpatrywanego obiektu. Drugie, najczęściej równania algebraiczne, wyrażają odpowiedź konkretnego materiału na zewnętrzne wymuszenia. Na podstawie powyższych związków tworzone jest równanie elementu (ang. element equation) opisujące odpowiedź w węzłach rozpatrywanego fragmentu obiektu na zadane obciążenie. Dla zagadnień mechanicznych równanie to związane jest z przemieszczeniem węzłów elementu. Dla problemów niestrukturalnych rozwiązaniem równania może być, przykładowo, temperatura w przewodnictwie cieplnym czy prędkość przepływu w dynamice płynów. Te węzłowe wartości zwane są pierwotnymi lub kinematycznymi niewiadomymi (ang. primary unknowns), natomiast inne, wyliczone na ich podstawie, takie jak naprężenia czy deformacje noszą nazwę wtórnych lub statycznych niewiadomych (ang. secondary unknowns). Równanie elementu przyjmuje postać:
(2)
gdzie [U] jest wektorem sił działających w węzłach, [k] oznacza macierz sztywności elementu, [D] jest wektorem odpowiedzi w węzłach.
Macierz sztywności elementu [k] uzyskiwana jest, dla prostych układów, z równań równowagi sił, natomiast w ogólności stosuje się takie metody jak:
metoda Rayleigh-Ritz, metoda residuów ważonych (metoda Galerkina, metoda najmniejszych kwadratów).
Należy pamiętać, iż wartości węzłowe nie są obliczane na tym etapie, lecz dopiero po stworzeniu globalnego równania opisującego cały obiekt. Ma to miejsce w punkcie czwartym.
Informacje węzłowe służą następnie do obliczenia reakcji na obciążenie w dowolnym punkcie elementu skończonego. Aby tego dokonać, wykorzystuje się uprzednio założoną funkcję interpolacyjną, zwaną także funkcją kształtu (ang.
interpolation function, shape function), która wiąże wartości węzłowe z wartością w danym punkcie w przestrzeni elementu. Funkcja kształtu ma najczęściej postać wielomianu, choć szeregi trygonometryczne bywają również stosowane. Poniżej pokazany zostanie związek funkcji interpolacyjnej z
40 rzędowością elementu skończonego na przykładzie analizy mechanicznej obiektu dwuwymiarowego.
Rysunek 11. Element skończony 2D pierwszego rzędu o trzech węzłach przed transformacją (linia ciągła) i po transformacji (linia przerywana). D1, D2,D3 oznaczają przesunięcia węzłów.
Na rysunku 11 przedstawiono element dwuwymiarowy pierwszego rzędu z liczbą węzłów równą trzy. Obiekt podlega niesztywnej transformacji, toteż każdy jego punkt doznaje przemieszczenia D wzdłuż osi X oraz Y zależnego od obu jego aktualnych współrzędnych. Zadaniem funkcji interpolacyjnej jest estymacja przemieszczenia dowolnego punktu wewnątrz elementu na podstawie wartości przemieszczeń skończonej liczby punktów opisujących ten element, które to przemieszczenia wyznaczono jako reakcje w węzłach. Do kompletnego opisu translacji dowolnego punktu w przestrzeni 2D rozpatrywanego elementu są zatem potrzebne dwa równania Dx(x,y) i Dy(x,y). Są to równania liniowe:
.
(3) Powyższy układ zawiera sześć nieznanych dotąd współczynników an, a zatem
potrzebnych jest sześć równań do ich ustalenia. Wymagany układ równań pozyskany jest na podstawie znanego przemieszczenia trzech wierzchołków elementu skończonego, a zatem:
(4)
41
.
Powyższe układy równań można zapisać w formie macierzowej:
Aby obliczyć wartości współczynników równań, należy dokonać następującego
przekształcenia:
Po obliczeniu współczynników an należy uzupełnić o nie równanie (3) i zapisać korzystając z notacji macierzowej:
42 ,
(9) gdzie Dx i Dy są składowymi wzdłuż osi X i Y przemieszczenia dowolnego
punktu wewnątrz elementu skończonego, Dn w macierzy [D] stanowią przemieszczenia w węzłach rozpatrywanego elementu, Nn w macierzy [N] są funkcjami kształtu.
Powyższe rozważania dotyczyły analizy elementu dwuwymiarowego o trzech wierzchołkach. Poniżej zostanie przedstawione analogiczne rozumowanie dla elementu dwuwymiarowego drugiego rzędu (rys.12).
Rysunek 12. Element skończony 2D drugiego rzędu o sześciu węzłach przed transformacją (linia ciągła) i po transformacji (linia przerywana). D1, D2,D3,D4,D5,D6 oznaczają przesunięcia węzłów.
W przypadku gdy opis elementu z użyciem funkcji liniowej jest niewystarczający, stosuje się wielomiany wyższych rzędów, np. funkcję kwadratową. Powstaje wówczas potrzeba znajomości przesunięć sześciu punktów, aby możliwe było sformułowanie funkcji kształtu. Podstawowe równanie opisujące przemieszczenie dowolnego punktu wewnątrz elementu ma zatem postać:
.
(10)
43 Po stworzeniu układu równań dla sześciu znanych węzłów i przekształceniu go do zapisu macierzowego otrzymujemy: Następnie wyznaczone zostają współczynniki an równań kwadratowych. Po
wprowadzeniu do równania (9) tak wyznaczonych współczynników równanie to przyjmuje formę: Pamiętając o związkach wyprowadzonych w poprzednim przykładzie, tj.
[A]=[H][D] i [N]=[R][H], równanie zawierające funkcje kształtu ma ostatecznie postać:
.
(13) Jak widać na przykładzie równania (8) oraz (12) funkcja kształtu wiąże
przemieszczenia węzłowe Dn z przesunięciami dowolnego punktu D w elemencie skończonym. Należy zauważyć, iż wybór typu funkcji opisującej
44 determinuje liczebność węzłów danego elementu. Innymi słowy, liczba stopni swobody posiadanych przez wszystkie węzły Dn musi być równa liczbie współczynników an równania opisującego przemieszczenie dowolnego punktu na podstawie wartości węzłowych. Warto też podkreślić, iż w praktyce wybór funkcji interpolacyjnej dokonuje się nie wprost na etapie wyboru typu elementów skończonych, na które rozpatrywany obiekt ma zostać podzielony.
3. Globalna macierz sztywności i warunki brzegowe
Gdy równania dla każdego elementu zostaną już ustanowione, tworzona jest globalna macierz sztywności [K] oraz globalne równanie w postaci:
, (14)
gdzie jest wektorem globalnych sił w węzłach, a oznacza wektor wszystkich stopni swobody w węzłach.
Zabieg agregacji poszczególnych równań elementu w jedność przeprowadzany jest z użyciem metody superpozycji zwanej tutaj metodą bezpośredniej sztywności (ang. direct stiffness method) i jest podyktowany potrzebą przedstawienia wszelkich wymuszeń i reakcji (np. przemieszczeń, orientacji itp.) w globalnym układzie współrzędnych oraz uzgodnieniem jednorodnej numeracji węzłów w całym obiekcie.
Na tym etapie można wykazać, iż macierz [K] jest osobliwa (jej wyznacznik równy jest zero), czego fizyczną interpretacją jest brak odebranych stopni swobody w którymkolwiek węźle w analizowanym obiekcie. Wprowadza się zatem warunki brzegowe, które rozwiązują problem osobliwości macierzy.
Dzieje się tak dlatego, iż warunki brzegowe oznaczają konkretne wartości rozpatrywanej zmiennej (podstawowe, kinematyczne warunki brzegowe) lub powiązanych z nią zmiennych, np. pochodnych (naturalne, statyczne warunki brzegowe), na granicy obszaru analizy np. struktury fizycznej.
4. Rozwiązanie globalnego równania
W tym kroku globalne równanie z punktu 3 zostaje rozwiązane metodami eliminacyjnymi (np. metoda Gaussa) lub iteracyjnymi (np. metoda Gaussa-Seidela). Wynik stanowią wartości w węzłach dla przemieszczeń, czyli
45 pierwotnych niewiadomych. Wtórne niewiadome, takie jak naprężenia, wylicza się następnie w oparciu o rozwiązanie równania z pkt 3.
MES znajduje zastosowanie w szeregu inżynierskich zagadnień takich jak mechanika strukturalna, przewodzenie ciepła, dynamika płynów i pole elektromagnetyczne. W ostatnich czasach implementacja metody elementów skończonych w zagadnieniach biologicznych także znajduje szerokie zastosowanie. W szczególności jest ona wykorzystywana do badania wytrzymałości mechanicznej naczyń krwionośnych takich jak tętniak aorty brzusznej, co jest dokładnie omówione w następnym rozdziale. W niniejszej pracy metoda elementów skończonych odgrywa rolę narzędzia w przeprowadzonych analizach, a jej rozwój nie jest przedmiotem rozprawy. Stosowana jest zatem zgodnie z aktualnie przyjętym standardem w odniesieniu do systemów biomechanicznych.
46