Metoda Elementów Skończonych
Strona – 1
METODA ELEMENTÓW SKOŃCZONYCH
– przykład rozwiązania kratownicy
materiały pomocnicze do zajęć z przedmiotu
Metody Komputerowe
Rzeszów 2007
Macierz kierunkowa
Transformacja
Metoda Elementów Skończonych
Strona – 3
Równanie metody elementów skończonych dla modelu w globalnym ukła-
dzie
odniesienia zapisać możemy jako:
K · a
= f ,
(1)
gdzie K jest globalną macierzą sztywności analizowanego układu, a wek-
torem parametrów węzłowych (przemieszczeń, obrotów - jeśli występują),
f
wektorem (węzłowych) obciążeń zewnętrznych.
Podobną zależność zapisać możemy również dla pojedynczego elementu
skończonego:
k
i
· V
i
= S
i
,
(2)
gdzie k
i
jest macierzą sztywności i-tego elementu w układzie globalnym,
V
i
wektorem parametrów węzłowych elementu, S
i
obciążeniem działa-
jącym w węzłach i-tego elementu.
Metoda Elementów Skończonych
Strona – 4
Przejście pomiędzy lokalnym układem współrzędnych (dla danego ele-
mentu) oraz globalnym układem współrzędnych (dla całej konstrukcji)
wymaga uwzględnienia różnicy kątów pomiędzy osiami obu układów.
Macierz kosinusów w przestrzeni 2D zapisać możemy jako
c
i
=
cos(x, ξ)
cos(x, η)
cos(y, ξ)
cos(y, η)
=
cos(α) − sin(α)
sin(α)
cos(α)
,
(3)
gdzie x, y są osiami globalnego układu
współrzędnych, ξ, η są osiami układu lo-
kalnego, zaś α wyraża kąt pomiędzy osia-
mi obu układów odniesienia.
Przyjmując dla uproszczenia założenie, że
c
= cos(α) oraz s = sin(α), równanie (3)
skraca się do postaci
c
i
=
c
−
s
s
c
.
Równanie MES
Macierz kierunkowa
Metoda Elementów Skończonych
Strona – 5
Aby transformować wektory obciążenia (pomiędzy układami odniesie-
nia), posłużyć należy się następującymi formułami podanymi w notacji
macierzowej
S
i
= C
i
· S
i,l
⇔
S
i,l
= C
T
i
· S
i
.
(4)
W identyczny sposób transformują się wektory przemieszczeń
V
i
= C
i
· V
i,l
⇔
V
i,l
= C
T
i
· V
i
.
(5)
Wymiar tzw. macierzy kierunkowej C
i
zależeć może np. od liczby stop-
ni swobody w węźle i liczby węzłów w danym elemencie. Przykładowo
dla elementu prętowego z dwoma stopniami swobody w każdym węźle,
macierz kierunkowa będzie mieć wymiar 4×4 i przyjmie postać
C
i
=
c
i
0
0
c
i
=
c
−
s
0
0
s
c
0
0
0
0
c
−
s
0
0
s
c
.
(6)
Równanie MES
Macierz kierunkowa
Metoda Elementów Skończonych
Strona – 5
Wówczas formułę transformacji wektora obciążeń (4) zapisać możemy
jako
h
S
1
S
2
i
1
h
S
3
S
4
i
2
|
{z
}
S
i
=
"
c
−s
0
0
s
c
0
0
0
0
c
−s
0
0
s
c
#
|
{z
}
C
i
·
h
N
ξ
T
η
i
1
h
N
ξ
T
η
i
2
|
{z
}
S
i,l
Zapisując równianie MES dla elementu w układzie lokalnym (analogicz-
nie do (2)) jako
S
i,l
= k
i,l
V
i,l
możemy podstawić je do równania (4), a uwzględniając zależność (5)
dostaniemy
S
i
= C
i
k
i,l
V
i,l
= C
i
k
i,l
C
T
i
|
{z
}
k
i
V
i
= k
i
V
i
Równanie MES
Macierz kierunkowa
Metoda Elementów Skończonych
Strona – 5
Przez analogię do równania MES spostrzec można, że równanie trans-
formacji lokalnej macierzy sztywności elementu w macierz globalną ele-
mentu zdefiniować możemy zależnością
k
i
= C
i
· k
i,l
· C
T
i
(4)
Zatem dla wspomnianego elementu prętowego z dwoma stopniami swo-
body w każdym węźle, przekształcenie macierzy lokalnej w globalną
przyjmie postać:
k
i
=
"
c
−s
0
0
s
c
0
0
0
0
c
−s
0
0
s
c
#
·
EA
l
"
1
0
−1
0
0
0
0
0
1
0
−1
0
0
0
0
0
# "
c
−s
0
0
s
c
0
0
0
0
c
−s
0
0
s
c
#
k
i
=
EA
l
c
2
cs
−c
2
−cs
cs
s
2
−cs
−s
2
−c
2
−cs
c
2
cs
−cs
−s
2
cs
s
2
=
EA
l
h
k
(i)
11
k
(i)
12
k
(i)
21
k
(i)
22
i
(5)
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Agregacja
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 7
Dla przykładu rozwiążmy za pomocą MES prosty układ kratownicowy,
którego schemat pokazano na poniższym rysunku.
Rysunek 1: Schemat kratownicy.
Analiza elementu
Macierze sztywności
elementu
Agregacja
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 8
Pierwszym etapem algorytmu MES jest dyskretyzacja (idealizacja)
rozważanego modelu. Polega to przede wszystkim na podziale konstrukcji
na elementy skończone, a następnie ponumerowaniu zdefiniowanych
w ten sposób węzłów i elementów.
Rysunek 2: Dyskretyzacja układu kratownicy.
Macierze sztywności
elementu
Agregacja
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 9
Następnie w każdym z węzłów
przyjmujemy odpowiednią do da-
nego typu elementu liczbę stopni
swobody. Dla elementu kratowni-
cowego (T=0) sytuacja ta poka-
zana została na rysunku obok.
Lokalną macierz sztywności i–
tego elementu kratownicowego
zapisać możemy jako:
Rysunek 3: Dodatnie kierunki sił
dla układu lokalnego i globalnego.
k
i,l
=
EA
l
1
0
−
1
0
0
0
0
0
1
0
−
1
0
0
0
0
0
(6)
Schemat
Dyskretyzacja
Analiza elementu
Agregacja
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 10
Korzystając z równania (4) otrzymujemy globalną macierz sztywności k
i
dla i–tego elementu w postaci (5):
ki =
EA
l
"
c2
cs
−c
2
−cs
cs
s2
−cs
−s
2
−c
2
−cs
c2
cs
−cs
−s
2
cs
s2
#
Zatem dla kolejnych elementów naszego modelu, nachylonych do układu
globalnego pod kątem α, dostaniemy
(α = 0◦ , c = 1, s = 0)
k2 = k4 = k6 =
EA
l
1
0
−1
0
0
0
0
0
−1
0
1
0
0
0
0
0
(α = 90◦ , c = 0, s = 1)
k1 = k5 =
EA
l
0
0
0
0
0
1
0
−1
0
0
0
0
0
−1
0
1
(α = 45◦ , c = s =
√
2
2
)
k3 = k7 =
EA
l
"
1
2
1
2
−
1
2
−
1
2
1
2
1
2
−
1
2
−
1
2
−
1
2
−
1
2
1
2
1
2
−
1
2
−
1
2
1
2
1
2
#
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 11
Kolejnym etapem algorytmu MES złożenie macierzy sztywności każdego
z elementów w globalną macierz całego układu wg zależności
K
= K
1
+ K
2
+ . . . + K
7
(7)
Macierze
sztywności
K
i
muszą
uwzględniać
zarówno
informa-
cję o początkowym i końcowym węźle danego elementu, jak
również położenie danego elementu w analizowanym układzie.
Konstruując zatem globalną macierz sztywności K
1
(dla elemen-
tu nr 1), wiemy że element ten ma swój początek w węźle nr 2
(wg globalnej numeracji węzłów) oraz koniec w węźle nr 1. Jeśli
przyjmiemy dla układu lokalnego, że początek jest punktem nr
1
,
a koniec
2
, wówczas:
✔
przypadkowi węzła 1-1 odpowiada lokalna numeracja
2-2
, co
wskazuje na element k
(1)
22
globalnej macierzy sztywności k
1
✔
połączeniu węzłów 1-2 odpowiada lokalna numeracja
2-1
, co
wskazuje na element k
(1)
21
globalnej macierzy sztywności k
1
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 11
Przechodząc do węzła nr 2 napisać możemy, że:
✔
przypadkowi węzła 2-2 odpowiada lokalna numeracja
1-1
, co wska-
zuje na element k
(1)
11
globalnej macierzy sztywności k
1
✔
dla przypadku połączenia 2-1 odpowiada lokalna numeracja
1-2
, co
wskazuje na element k
(1)
12
globalnej macierzy sztywności k
1
Na tej podstawie zbudować możemy teraz globalną macierz K
1
, umiesz-
czając podmacierze k
(1)
ij
na miejscach odpowiadającym indeksom wyzna-
czonym przez numery węzłów, jak następuje
K
1
=
EA
l
k
(1)
22
k
(1)
21
0
0
0
k
(1)
12
k
(1)
11
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
gdzie wymiar macierzy K
1
odpowiada liczbie węzłów (5 × 5).
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 11
Dla elementu nr 2 (poziomy), kolejnym indeksom odpowiadają macierze:
✔
1-1 ⇒ k
(2)
11
✔
1-3 ⇒ k
(2)
12
✔
3-3 ⇒ k
(2)
22
✔
3-1 ⇒ k
(2)
21
Zatem globalna macierz K
2
dla elementu nr 2 przyjmuje postać
K
2
=
EA
l
k
(2)
11
0
k
(2)
12
0
0
0
0
0
0
0
k
(2)
21
0
k
(2)
22
0
0
0
0
0
0
0
0
0
0
0
0
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 11
Dla elementu nr 3 (skośny), kolejnym indeksom odpowiadają macierze:
✔
2-2 ⇒ k
(3)
11
✔
2-3 ⇒ k
(3)
12
✔
3-3 ⇒ k
(3)
22
✔
3-2 ⇒ k
(3)
21
Zatem globalna macierz K
3
dla elementu nr 3 przyjmuje postać
K
3
=
EA
l
0
0
0
0
0
0
k
(3)
11
k
(3)
12
0
0
0
k
(3)
12
k
(3)
22
0
0
0
0
0
0
0
0
0
0
0
0
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 11
Dla elementu nr 4 (poziomy), kolejnym indeksom odpowiadają macierze:
✔
2-2 ⇒ k
(4)
11
✔
2-4 ⇒ k
(4)
12
✔
4-4 ⇒ k
(4)
22
✔
4-2 ⇒ k
(4)
21
Zatem globalna macierz K
4
dla elementu nr 4 przyjmuje postać
K
4
=
EA
l
0
0
0
0
0
0
k
(4)
11
0
k
(4)
12
0
0
0
0
0
0
0
k
(4)
21
0
k
(4)
22
0
0
0
0
0
0
Analogicznie tworzymy macierze K
i
dla pozostałych elementów (5 do 7),
co zaleca się wykonać samodzielnie.
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 11
Ponieważ w przedstawionym na poprzednich slajdach zapisie, macierza-
mi są zarówno składniki k
(e)
ij
(por. równ. 5) jak również (pogrubione)
zera, stąd po uwzględnieniu wartości macierzy k
4
, macierz sztywności
K
4
ostatecznie przyjmuje postać:
K
4
=
EA
l
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
−
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
−
1
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 11
Mając zdefiniowane macierze sztywności K
i
dla wszystkich elementów
(i = 1 ÷ 7) i dokonując ich sumowania K =
P
7
i
=1
K
i
(por. równ. 7),
otrzymujemy globalną macierz sztywności K całego układu:
K =
EA
l
1
0
0
0
−
1
0
0
0
0
0
0
1
0
−
1
0
0
0
0
0
0
0
0
3
2
1
2
−
1
2
−
1
2
−
1
0
0
0
0
1
3
2
1
2
−
1
2
−
1
2
0
0
0
0
−
1
0
−
1
2
−
1
2
5
2
1
2
0
0
−
1
0
0
0
−
1
2
−
1
2
1
2
1
2
0
−
1
0
0
0
0
−
1
0
0
0
3
2
1
2
−
1
2
−
1
2
0
0
0
0
0
−
1
1
2
1
2
−
1
2
−
1
2
0
0
0
0
−
1
0
−
1
2
−
1
2
3
2
1
2
0
0
0
0
0
0
−
1
2
−
1
2
1
2
1
2
Uwaga: Otrzymana macierz K jest macierzą osobliwą (det K = 0)! Nie
jest możliwe zatem wyznaczenie (na tym etapie) macierzy odwrotnej K
−
1
=
1
det K
K
D
. Do tej pory, przy budowie macierzy K, uwzględnione zostały warun-
ki równowagi oraz zgodności przemieszczeń, natomiast nie narzucono żadnych
ograniczeń kinematycznych (warunków brzegowych).
Wektor przemieszczeń i obciążenia
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Agregacja
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 12
Przyjmijmy następujące założenia:
✔
niech wektor przemieszczeń a zawiera składowe przemieszczenia
w układzie globalnym) każdego z węzłów analizowanego układu,
✔
oraz na wektor sił f niech składają (znane) siły potencjalnie wystę-
pujących w węzłach.
Powyższe sformułowania zapisać możemy jako
a
=
a
1x
a
1y
a
2x
a
2y
a
3x
a
3y
a
4x
a
4y
a
5x
a
5y
,
f
=
f
1x
f
1y
f
2x
f
2y
f
3x
f
3y
f
4x
f
4y
f
5x
f
5y
.
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Agregacja
Wektor przemieszczeń
i obciążenia
Rozwiązanie równania
MES
Wyznaczenie sił
przekrojowych
Metoda Elementów Skończonych
Strona – 13
Z uwagi na warunki podparcia naszego układu zapisać możemy następu-
jące warunki brzegowe:
✔
a
2x
= 0,
✔
a
2y
= 0,
✔
a
4y
= 0.
Wartości tych jednak nie podstawiamy wprost do wektora przemiesz-
czeń a, lecz wpisujemy je w odpowiednie miejsca wektora sił f . W konse-
kwencji pozwala to zmodyfikować wartości macierzy sztywności K w ta-
ki sposób, że w miejscach przecięcia się wybranych wierszy i kolumn
wstawiamy wartość 1, natomiast pozostałe wartości w tych wierszach
i kolumnach zerujemy.
Dzięki temu macierz K staje się macierzą nieosobliwą i możliwe jest
wyznaczenie K
−1
. Ostateczny wygląd macierzy sztywności K oraz wek-
tora obciążenia f pokazany został na kolejnym slajdzie. Zmodyfikowane
składniki wyróżniono kolorem
czerwonym
.
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Agregacja
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Metoda Elementów Skończonych
Strona – 14
Po uwzględnieniu warunków brzegowych macierz sztywności K oraz wek-
tor obciążenia f przyjmą wartość:
K =
EA
l
1
0
0
0
−1
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
−1
0
0
0
5
2
1
2
0
0
−1
0
0
0
0
0
1
2
1
2
0
0
0
0
0
0
0
0
0
0
3
2
0
−
1
2
−
1
2
0
0
0
0
0
0
0
1
0
0
0
0
0
0
−1
0
−
1
2
0
3
2
1
2
0
0
0
0
0
0
−
1
2
0
1
2
1
2
,
f =
f1x
f1y
f2x
f2y
f3x
f3y
f4x
f4y
f5x
f5y
=
0
0
0
0
0
0
0
0
0
−P
.
W ten sposób pozbyliśmy się osobliwości macierzy K. Przekształcając
następnie równianie (1) do postaci
a
= K
−1
f
,
wyznaczyć możemy nieznane przemieszczenia węzłów rozpatrywanego
układu. Dalej, znając przemieszczenia a, określić możemy siły wewnętrz-
ne każdego z elementów.
Schemat
Dyskretyzacja
Analiza elementu
Macierze sztywności
elementu
Agregacja
Wektor przemieszczeń
i obciążenia
Warunki brzegowe
Rozwiązanie równania
MES
Metoda Elementów Skończonych
Strona – 15
Weźmy dla przykładu elementu nr 2, dla którego zapisać możemy rów-
nanie S
2
= k
2
V
2
. Wynika z niego, że:
S
1
S
2
S
3
S
4
= k
2
·
a
1x
a
1y
a
3x
a
3y
Korzystając z równań transformacji (4) lub (5) łatwo możemy przejść do
układu lokalnego i wyznaczyć wartości sił przekrojowych Q oraz N:
N
1
Q
1
N
2
Q
2
= C
T
i
·
S
1
S
2
S
3
S
4
lub
N
1
Q
1
N
2
Q
2
= k
2,l
· C
T
i
·
a
1x
a
1y
a
3x
a
3y
Tym sposobem dobrnęliśmy do końca przykładu. Z uwagi na ludzką niedoskonałość prezentacja
może zawierać pewne błędy numeryczne. Dlatego w trakcie samodzielnej pracy nigdy nie może
zabraknąć czujności.
