Metoda Elementów Skończonych w Modelowaniu Układów Mechatronicznych

(1)

_________________________________________________________________

Metoda Elementów Skończonych w

Modelowaniu Układów Mechatronicznych

_________________________________________________________________

Układy prętowe (Scilab)

(2)

I. MES 1D – układy prętowe. Podstawowe informacje

Istotą metody elementów skończonych jest sposób aproksymacji cząstkowych równań różniczkowych (które określają zachowanie funkcji niewiadomych wewnątrz rozpatrywanego obszaru dla danego zjawiska fizycznego), polegający na podziale obszaru obliczeniowego na małe podobszary o prostych kształtach, zwane elementami skończonymi oraz specjalny sposób konstruowania funkcji aproksymujących (tzw. funkcji kształtu) opierający się na funkcjach zdefiniowanych w elementach skończonych. MES jako jedna z niewielu metod potrafi modelować zjawiska w skomplikowanych obszarach obliczeniowych, co stanowi jedną z jej podstawowych zalet w zastosowaniach praktycznych. ^(*)

ALGORYTM POSTĘPOWANIA w MES

________________________________________________________________________________

(*) cyt. oprac. na podstawie – K.Banaś „Wprowadzenie do MES” http://www.metal.agh.edu.pl/~banas/wprowadzenie_do_MES.pdf

Odczyt rezultatów

(przemieszczenia w węzłach), wyznaczenie reakcji w podporach, obliczenie naprężeń w prętach konstrukcji (lub ich odpowiedników dla innych zjawisk fizycznych)

Rozwiązanie układu równań MES

KQ=F (statyka, liniowa spręzystość)

Uwzględnienie warunków brzegowych

(np. metoda kary)

Agregacja

utworzenie globalnej macierzy sztywności, oraz wektora obciążeń sztywności (lub jej odpowiedników dla innych zjawisk fizycznych)

Wyznaczenie lokalnych macierzy sztywności

(lub jej odpowiedników dla innych zjawisk fizycznych, np. macierzy przewodności cieplnej)

Dyskretyzacja

(podział na elementy skończone)

Przyjęcie modelu fizycznego rozważanego obiektu

Tworzenie/import geometrii, określenie typu stosowanych elementów skończonych

(3)

II. Przykład modelowania układu prętowego za pomocą MES

1. Przyjęcie modelu fizycznego rozważanego obiektu,

Pręt stalowy o zmiennym przekroju, utwierdzony jednym końcem, obciążony dwoma siłami.

2. Dyskretyzacja - podział na elementy skończone

3. Wyznaczenie lokalnych macierzy sztywności

Lokalna macierz sztywności dla elementu prętowego:

gdzie: Ae – pole przekroju elementu skończonego, Ee – moduł Younga, Le – długość elementu skończonego

W przykładzie należy wyznaczyć 4 lokalne macierze sztywności, ponieważ wyróżniono 4 elementy skończone

Ponieważ cały pręt jest wykonany z jednego materiału to E1 = E2 = E3 = E4 = E

4. Agregacja

Utworzenie globalnej macierzy sztywności K

oraz wektora obciążeń

w ogólności wektor obciążeń F może zawiera również obciążenia masowe (f^e) oraz równomiernie rozłożone T^e (tzw. trakcje)

1 1

e e e

e

A E L

  

     

k

e e

 

K k



^e ^e



^e

e

   

F f T F

1 1 4 4

1 1 11 12 4 4 11 12

1 1 4 4

1 21 22 4 21 22

1 1 1 1

...

1 1 1 1

k k k k

A E A E

L k k L k k

     

   

                     

k k

(4)

Zagregowana macierz sztywności dla przykładu ma postać:

Uwagi:

 Zauważ, że „łączenia” lokalnych macierzy sztywności występują w węzłach które są wspólne dla sąsiadujących elementów. Tworzy się tzw. pasmo, zaznaczone powyżej liniami

przerywanymi.

 Tam gdzie nie występuje „połączenie” poszczególnych węzłów elementem, tam w macierzy sztywności występują wartości zerowe (np. 1-3, 1-4, 1-5, 2-4, itd.). Macierz posiadająca wiele elementów zerowych nazywa się macierzą rzadką.

 Poza tym, że macierz sztywności jest macierzą rzadką i pasmową jest ponadto macierzą symetryczną, tj.

Zagregowany wektor obciążeń (macierz kolumnowa) dla przykładu ma postać:

5. Uwzględnienie warunków brzegowych – metoda kary

Metoda kary (ang. Penalty approach) polega na zadaniu dużej sztywności w miejscach określenia warunku przemieszczeniowego. W ogólności, modyfikacji podlegają zarówno, globalna macierz sztywności, jak i globalny wektor obciążeń wg następującego schematu:

1. Modyfikacja globalnej macierzy sztywności K

dodajemy stałą „C” do elementów leżących na przekątnej macierzy odpowiadającym odebraniu danego stopnia swobody

2. Modyfikacja macierzy kolumnowej obciążeń F

dodajemy iloczyn „C u1” do elementu odpowiadającemu odebraniu przemieszczenia

stała C dobierana jest zwykle jako iloczyn modułu maksymalnej wartości macierzy sztywności oraz pewnego mnożnika np. 10⁴ (duża sztywność) , tj.

max|Kij|*10⁴

1 1

11 12

1 1 2 2

21 22 11 12

2 2 3 3

21 22 11 12

3 3 4 4

21 22 11 12

4 4

21 22

0 0 0

0 0

0 0 0

k k

k k k k

k k

 

  

 

 

 

  

 

 

K

ij



ji

K K

0 0 0 2

P

 

  

 

  

 

 

F

(5)

W rozważanym przykładzie utwierdzenie występuje w węźle nr 1, więc modyfikacja układu równań MES będzie następująca:

Uwaga:

W ogólności zadawane przemieszczenia nie muszą mieć zerowych wartości (np. luz montażowy).

Układ równań uwzględniający modyfikację wektora obciążeń ma wtedy postać:

6. Rozwiązanie układu równań MES

Zastosowanie metody kary umożliwia potraktowanie wektora przemieszczeń jako niewiadomej w równaniu MES:

Odpowiada to przekształceniu równania do postaci

Uwaga:

Zagregowana macierz sztywności bez uwzględniania warunków brzegowych jest macierzą osobliwą.

Uwzględnienie warunków brzegowych (tutaj metodą kary) powoduje, że macierz ta przestaje być osobliwa i możliwe jest rozwiązanie powyższego układu równań.

7. Odczyt rezultatów - przemieszczenia w węzłach, obliczenie naprężeń w poszczególnych prętach

Rozwiązanie układu równań MES wyznacza przemieszczenia węzłowe (wektor Q) Naprężenie w poszczególnych prętach wyznacza się z zależności

gdzie:

q^e jest dwuelementowym wektorem przemieszczeń dla elementu, B^e jest macierzą geometryczną elementu określoną zależnością:

11 12 1 1 1 1

21 22 2 2 2

1 2

C ... C

...

N N

N N NN N N

K K K Q F u

K K K Q F

  

     

     

      

     

     

1 1

11 12 1

1 1 2 2

21 22 11 12 2

2 2 3 3

21 22 11 12 3

3 3 4 4

21 22 11 12 4

4 4

21 22 5

C 0 0 0 0

0 0

0 0 0

0 0 0 2

k k Q

Q P

k k k k

k k k k Q

Q P

k k

      

       

     

       

      

     

         

 

KQ = F

1

Q = K F

E

e

e e e

σ = B q

 

1 1 1

Le

 

Be

(6)

III. Kod do programu SciLab dla powyższego przykładu

Dane do przykładu:

Moduł Younga E= 200 000MPa Siła P=1000N

Długości poszczególnych elementów: L1=100mm; L2=50mm; L3=200mm; L4=100mm Pola przekroju poszczególnych elementów: A1=50mm²; A2=40mm²; A3=20mm²; A4=5mm²

//określenie wartości obciążenia P=1000;

//określenie długości poszczególnych elementów L1=100;

L2=50;

L3=200;

L4=100;

//określenie własności materiałowych E=2e5;

//określenie pól przekroju dla poszczególnych elementów A1=50

A2=40;

A3=20;

A4=5;

//wyznaczenie macierzy sztywności poszczególnych elementów K1=E*A1/L1*[1, -1; -1, 1];

K2=E*A2/L2*[1, -1; -1, 1];

K3=E*A3/L3*[1, -1; -1, 1];

K4=E*A4/L4*[1, -1; -1, 1];

//Zainicjowanie globalnej macierzy sztywności KG=zeros(5,5);

//Agregacja globalnej macierzy sztywności KG(1:2,1:2)=KG(1:2,1:2)+K1;

KG(2:3,2:3)=KG(2:3,2:3)+K2;

KG(3:4,3:4)=KG(3:4,3:4)+K3;

KG(4:5,4:5)=KG(4:5,4:5)+K4;

//Określenie wektora sił węzłowych F=[0; -P; 0; 0; 2*P];

//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;

//Rozwiązanie układu równań Q=inv(KG)*F;

//określenie wektora przemieszczeń elementu q1=Q(1:2);

q2=Q(2:3);

q3=Q(3:4);

q4=Q(4:5);

//określenie macierzy geometrycznych dla poszczególnych elementów B1=1/L1*[-1,1];

B2=1/L2*[-1,1];

B3=1/L3*[-1,1];

B4=1/L4*[-1,1];

//obliczenie naprężeń w poszczególnych elementach sig1=E*B1*q1;

sig2=E*B2*q2;

sig3=E*B3*q3;

sig4=E*B4*q4;

//wyświetlenie wektora przemieszczeń Q disp(Q);

//wyświetlenie naprężeń w elementach disp(sig1, "Sigma 1=");

disp(sig2, "Sigma 2=");

(7)

IV. Zadania do samodzielnego wykonania

1. Zapoznaj się dokładnie ze sposobem modelowania układów prętowych MES (pkt I i II) 2. Uruchom kod z pkt. III w programie SciLab oraz przeanalizuj jego działanie.

3. Oblicz analitycznie wydłużenia oraz naprężenia dla poszczególnych części pręta wykorzystując dane z pkt. III.

Rezultaty obliczeń analitycznych porównaj z wynikami otrzymanymi w programie Scilab.

4. Napisz własny skrypt do programu SciLab wyznaczający przemieszczenia węzłowe oraz naprężenia w prętach dla następującego przykładu:

Dane do przykładu:

Moduł Younga E¹= 200 000MPa (stal); E²= 70 000MPa (aluminium);

Siła P=2000N

Długości poszczególnych elementów: L1=100mm; L2=70mm;

Pola przekroju elementów: A1=60mm²; A2=30mm²;

5. Oblicz analitycznie przemieszczenie przekroju w którym przyłożona jest siła P oraz naprężenia dla obydwu części pręta.

Rezultaty obliczeń analitycznych porównaj z wynikami otrzymanymi w programie Scilab.