_________________________________________________________________
Metoda Elementów Skończonych w
Modelowaniu Układów Mechatronicznych
_________________________________________________________________
Układy belkowe (Scilab)
str.1
I. Podstawowe informacje
Dla elementów belkowych oprócz rozciągania (ściskania) występuje również zginanie oraz ścinanie.
W węzłach elementów belkowych wystąpią dwa stopnie swobody * (dla elementów prętowych był to tylko jeden stopień swobody – przemieszczenie wzdłuż osi x).
W każdym z węzłów wyróżniamy stopnie swobody przemieszczeniowe (indeksy nieparzyste) oraz przemieszczenia kątowe (indeksy parzyste).
Macierz sztywności płaskiego elementu belkowego ma postać:
Dla wektora obciążeń F indeksy nieparzyste odpowiadają siłom węzłowym, natomiast nieparzyste momentom skupionym w odpowiednich węzłach.
Agregacja, czyli łączenie lokalnych macierzy sztywności w globalny układ równań uwzględnia wszystkie stopnie swobody w węzłach.
W metodzie kary w zależności od typu podpory modyfikujemy odpowiednie elementy macierzy uwzględniające zarówno przemieszczenia jak i obroty w węzłach.
Po rozwiązaniu macierzowego równania MES
otrzymujemy wektor węzłowych przemieszczeń Q
Pozostałe wielkości, jak reakcje podporowe, odkształcenia oraz naprężenia w elementach wyznaczane są z odpowiednich formuł. W ramach tego laboratorium obliczenia MES dla belek dotyczyć będą wyłącznie wyznaczania przemieszczeń.
________________________________________________________________________________
(*) tak jest tylko w przypadku rozważanych płaskich elementów belkowych. W ogólnym przypadku przestrzennym, w którym tego typu elementami można modelować również ramy, w każdym z węzłów wystąpi 6 stopni swobody (3 przemieszczenia oraz 3 obroty)
3 2 3 2
2
3 2
12 6 12 6
4 6 2
12 6
4
e z
L L L L
L L L
EI
L L
sym L
−
−
=
−
k
KQ = F
str.2
II. Przykład modelowania układu belkowego
Dana jest belka przedstawiona na rysunku poniżej. Belka o przekroju prostokątnym podparta jest na końcach za pomocą podpór przegubowych, w środku obciążona siłą skupioną P
Model dyskretny składa się z 4 elementów skończonych (5 węzłów), tak więc układ w sumie posiada 10 stopni swobody. Na rysunku zaznaczono również warunki brzegowe dla zadania. Odebranie przemieszczeń liniowych dla węzła 1 oraz 5 (indeksy (1 i 9) oraz siłę w węźle nr 3 (indeks 5)
Na poniższym rysunku zaznaczono poglądowo macierzowy układ równań dla przykładu. Jasnoszare pola w globalnej macierzy sztywności odpowiadają kolejnym lokalnym macierzom sztywności elementu belkowego. Ciemnoszarym kolorem zaznaczono miejsce gdzie wartości ulegają zsumowaniu. Jak wspomniano wcześniej dotyczy to obydwu stopni swobody w węźle. Czarne pola w macierzy sztywności to miejsca w których modyfikacji podlega (zgodnie z zasadami metody kary) element macierzy w którym odebrano przemieszczeniowy stopień swobody. Czarne pole w globalnym wektorze obciążeń odpowiada zadanej sile P.
Uwaga:
Dla utwierdzeń, oprócz przemieszczeń liniowych odbierane są również przemieszczenia kątowe.
W takim przypadku należy to odpowiednio uwzględnić (indeksy parzyste). W wektorze obciążeń parzyste indeksy odpowiadają momentom skupionym w odpowiednich węzłach
str.3
III. Kod do programu SciLab dla powyższego przykładu
Dane do przykładu:
Moduł Younga E = 200 000MPa Siła P = 200N
Długości poszczególnych elementów: L = 100mm
Parametry przekroju poprzecznego belki b = 5mm, h = 10mm
clear
clc //określenie wartości obciążenia P=200;
//określenie długości elementu (podział jednorodny) L=100;
//określenie własności materiałowych E=2e5;
//obliczenie momentu bezwładności elementów (taki sam dal całej belki) b=5;
h=10;
Iz=b*h^3/12;
//wyznaczenie macierzy sztywności poszczególnych elementów K=E*Iz*[12/L^3, 6/L^2, -12/L^3, 6/L^2;
6/L^2, 4/L, -6/L^2, 2/L;
-12/L^3, -6/L^2, 12/L^3, -6/L^2;
6/L^2, 2/L, -6/L^2, 4/L;];
//Zainicjowanie globalnej macierzy sztywności KG=zeros(10,10);
//Agregacja globalnej macierzy sztywności KG(1:4,1:4)=KG(1:4,1:4)+K;
KG(3:6,3:6)=KG(3:6,3:6)+K;
KG(5:8,5:8)=KG(5:8,5:8)+K;
KG(7:10,7:10)=KG(7:10,7:10)+K;
//Określenie wektora sił węzłowych F=[0; 0; 0; 0; P; 0; 0; 0; 0; 0];
//Określenie przemieszczeniowych warunków brzegowych - metoda kary C=max(KG)*10^4; //Określenie stałej C modelującej sprężynę o dużej sztywności KG(1,1)=KG(1,1)+C;
KG(9,9)=KG(9,9)+C;
//Rozwiązanie układu równań Q=inv(KG)*F;
disp(Q)