Elementy wyższych rzędów
Oprócz elementów liniowych (dwuwęzłowych), w których zakłada Numeryczne Modelowanie Układów Ciągłych
się liniową zmianę wartości wielkości analizowanej wzdłuż krawędzi Podstawy Metody Elementów Skończonych
elementu, można stosować elementy wyższych rzędów.
Część II
Postępując podobnie jak poprzednio można wyprowadzić równania elementu dla większej ilości węzłów w elemencie prętowym w lokalnym układzie odniesienia.
Ireneusz Czajka
AGH, WIMiR
u
u
u
u
u
1
2
1
2
3
1
2
1
3
2
marzec 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Elementy wyższych rzędów
Interpolacja
Rozkład przemieszczeń wzdłuż pręta niech będzie wyznaczony przez trzy węzły – czyli może być opisany parabolą.
Niech zostanie przyjęty lokalny układ współrzędnych ze zmienną Algorytm
ξ ∈< 0 , 1 > i węzłami na końcach elementu oraz na jego środku.
Węzły mają współrzędne ξ 0 = 0, ξ 1 = 0 . 5, ξ 2 = 1. W takiej 1
wyznaczyć (przyjąć) rozkład przemieszczeń wzdłuż pręta sytuacji można napisać następujące funkcje kształtu, wynikające 2
wyznaczyć macierz sztywności ze znanego wzoru
wprost z interpolacji Lagrange’a
3
wyznaczyć macierz bezwładności ze znanego wzoru
ξ
ξ
ξ
ξ
ξ
ξ
N
− ξ
− x
− ξ
− ξ
− ξ
− ξ
= [ l
1
2
0
2
0
1
0
l 1 l 2] =
=
ξ 0 − ξ 1 ξ 0 − ξ 2
ξ 1 − ξ 0 ξ 1 − ξ 2
ξ 2 − ξ 0 ξ 2 − ξ 1
= [1 − 3 ξ + 2 ξ 2 4( ξ − ξ 2) 2 ξ( ξ − 1)] (1) Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Macierz sztywności
Macierz bezwładności
Przy takiej postaci funkcji kształtu, można wyznaczyć macierz sztywności, pamiętając, iż
dN
Podobnie macierz bezwładności można wyliczyć jako
[ B] = ldξ
Z
4
2
− 1
ρAl
czyli
[ m] =
NT ρNdΩ =
2
16
2
30
(3)
1 h
i
Ω
− 1 2
4
[ B] = l − 3 + 4 ξ 4 − 8 ξ 4 ξ − 1
W ten sposób można wyznaczać macierze sztywności oraz co daje macierz sztywności w postaci
bezwładności dla elementów wyższych rzędów. Warto pamiętać, iż
Z
wyższy rząd elementu daje znaczne korzyści obliczeniowe.
EA
7
− 8
1
[ k] =
[ B] T E [ B] dΩ =
3 l − 8 16 − 8
(2)
Ω
1
− 8
7
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Element trzeciego stopnia
Macierz sztywności
Funkcje kształtu są następujące
Macierze sztywności i bezwładności przyjmjują następującą postać
− 9
2 ξ 3 + 9 ξ 2 − 11
2 ξ + 1
27
37
− 189
27
− 13
NT
10
40
20
40
=
2 ξ 3 − 45
2 ξ 2 + 9 ξ
(4)
− 189
54
− 297
27
− 27
AE
2 ξ 3 + 18 ξ 2 − 92 ξ
[ k] =
40
5
40
20
(5)
9
l 27
− 297
54
− 189
2 ξ 3 − 92 ξ 2 + ξ
20
40
5
40
− 13
27
− 189
37
40
20
40
10
Następnie różniczkując wyznaczamy macierz zgodności geometrycznej ( B) oraz całkując macierz bezwładności.
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Macierz bezwładności
Podstawowe zagadnienia dynamiki
Metoda elementów skończonych znajduje zastosowanie nie tylko do rozwiązywania problemów statyki. Metoda pozwala także na rozwiązywanie problemów dynamicznych. Równanie opisujące Jak poprzednio całkując funkcje kształtu po objętości całego dynamikę układu z tłumieniem proporcjonalnym do prędkości, elementu uzyskujemy macierz bezwładności.
zdyskretyzowanego elementami skończonymi przedstawia się
następująco
8
33
− 3
7
105
560
140
619
33
27
− 27
− 3
[ m]
[ m] = ρAl
{¨ u
560
70
560
140
(6)
i } + [ c] { ˙ ui } + [ k] {ui } = {Pi }, (7)
− 3
− 27
27
33
140
560
70
560
7
− 3
33
8
gdzie kropkami nad zmienną oznaczono pochodną względem 619
140
560
105
czasu, czyli
d 2
{¨ ui } = dt 2 {ui}, d
{ ˙ ui } = dt {ui}.
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Częstości własne
Macierze bezwładności i sztywności można wyznaczyć bez większych problemów, natomiast najwięcej problemów sprawia macierz tłumienia. Jej postać można znaleźć na podstawie zależności
Z
Jedno z najczęściej spotykanych
[ c] =
[ N] T µ[ N] dΩ
Zagadnienie własne, czyli wyznaczenie częstości drgań własnych Ω
układu.
Często przyjmuje się tzw. proporcjonalne tłumienie, czyli macierz Zazwyczaj nie uwzględnia się tłumienia, ponieważ tłumienia wyraża się przy pomocy macierzy sztywności i bezwładności
nieznane współczynniki tłumienia
[ c] = k[ m] + w [ k]
znacznie większa trudność rozwiązania problemu własnego gdzie współczynnik w oznacza współczynnik tłumienia niewielkie różnice częstości własnych z tłumieniem i bez wewnętrznego – materiałowego, zaś k współczynnik tłumienia konstrukcyjnego – zewnętrznego.
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Problem własny
Po pomnożeniu równania (8) lewostronnie przez [ m] − 1 można zapisać
[ m] {¨ ui } + [ k] {ui } = 0
[ m] − 1[ k] − ω 2I {ui } = 0
(9)
dla rozwiązania harmonicznego ui = Aeiωt daje można zastosować podstawienia λ = ω 2, A = [ m] − 1[ k], X = {ui}
{¨ ui } = −ω 2 {ui }
uzyskując następującą, dobrze znaną z matematyki, postać problemu znajdowania wartości i wektorów własnych czyli
[ k] − ω 2[ m] {ui } = 0 .
(8)
(A − λ I)X = 0
Rozwiązanie tego problemu można sprowadzić do znalezienia Po rozwiązaniu tak sformułowanego zagadnienia wartości własne pierwiastków równania charakterystycznego
λi stanowią kwadraty częstości własnych. Wektory własne X i przedstawiają postacie drgań odpowiadające odpowiednim
[ k] − ω 2[ m] = 0
wartościom własnym.
i wyznaczenia postaci drgań własnych {ui}.
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Drgania swobodne pręta
Uwaga
Wartości własne mają konkretne wartości, ale dla wektorów własnych czyli postaci drgań własnych istotne są nie absolutne wartości dla poszczególnych węzłów, lecz stosunki tych wartości.
Zatem można, a nawet należy dokonać normalizacji tych wektorów.
Zazwyczaj normalizuje się je do jedynki lub względem masy. Przy Znaleźć pierwszą częstość drgań własnych stalowego pręta normalizacji do jedynki postępuje się według następującej utwierdzonego lewostronnie, o długości l = 1 m, ρ = 7 , 65 kg/m 3, zależności
E = 2 , 1 · 105 MPa, A = 25 cm 2. Ocenić dokładność rozwiązania
{u
w zależności od liczby elementów prętowych przyjętych do
{u
i }
i }N =
(10)
max( u
dyskretyzacji.
i )
Wpraktyce często stosuje się normalizację względem masy. Przy normalizacji względem masy wykorzystywane jest równanie
{ui }Tj [ m] {ui}j = I
gdzie {ui}j jest j-tym wektorem własnym.
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Rozwiązanie
Pierwsza częstość drgań własnych pręta wynosi
s
s
Ważne
π
E
π
2 , 1 · 1011
rad
ω
Przy problemach drgań własnych niezwykle istotne jest nałożenie 1 =
=
= 8230
2 l
ρ
2 · 1 7 , 65 · 103
s
warunków brzegowych, ponieważ postacie drgań własnych zależą od warunków brzegowych.
Numeryczne wyznaczenie częstości drgań własnych Bardzo dobrze to widać przy belce jednostronnie utwierdzonej, gdzie pierwsza postać drgań własnych bez uwzględnienia warunków Aby wyznaczyć częstość drgań własnych należy
brzegowych jest zdecydowanie różna od tej z uwzględnieniem dokonać dyskretyzacji
utwierdzenia.
dokonać agregacji macierzy sztywności
dokonać agregacji macierzy bezłwadności
nałożyć warunki brzegowe Dirichleta (utwierdzenie) Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Rozwiązanie
Rozwiązanie
Agregacja macierzy sztywności przebiega dokładnie jak Wyniki w zależności od ilości elementów skończonych wyglądają w poprzednim przykładzie. Identycznie przebiega agregacja następująco
macierzy bezwładności.
12
Dla modelu z trzema elementami skończonym macierze przyjmują postać następującą:
10
1
0
0
0
8
0
3 . 15 · 109
− 1 . 575 · 109
0
[ k] =
0 − 1 . 575 · 109
3 . 15 · 109
− 1 . 575 · 109
6
Blad [%]
0
0
− 1 . 575 · 109
1 . 575 · 109
4
1
0
0
0
0
4 . 25
1 . 0625
0
[ m] =
2
0 1 . 0625
4 . 25
1 . 0625
0
0
1 . 0625
2 . 125
0
0
10
20
30
40
50
Elementow
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Elementy paraboliczne
Elementy sześcienne
Tok obliczeń jest dokładnie taki sam jak dla elementów liniowych.
Przy czym macierze sztywności i bezwładności mają postać jak w zależnościach poprzednio wyprowadzonych.
Po zdyskretyzowaniu pręta przy pomocy jednego elementu trzeciego stopnia, pierwsza częstość drgań własnych została Ilość elementów ω [rad/s] Błąd [%]
obliczona jako
1
8260,9
0,375
ω 1 = 8230 , 5
2
8232,1
0,026
3
8230,4
0,005
co daje błąd około 0 , 006 %.
4
8230,1
0,001
5
8230,0
–
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Elementy wyższych rzędów
Elementy wyższych rzędów
Wnioski
No to
Dla elementów liniowych podział na dwa elementy daje dokładność wyznaczenia pierwszej częstości drgań własnych około 3%. Podział na trzy elementy daje dokładność rzędu 1%.
Jednocześnie widać, iż dla elementów liniowych dopiero podział na 111 elementów skończonych dał dokładną wartość pierwszej częstości drgań własnych, podczas gdy podział na 5 elementów Może na dziś koniec
parabolicznych daje dokładną wartość pierwszej częstości drgań własnych.
Metody adaptacyjne
Adaptacyjny MES występuje zasadniczo w dwóch odmianach typu h – zagęszcza się siatkę
typu p – podnosi się stopień funkcji kształtu Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011
Ireneusz Czajka
Fakultet Numeryczne Modelowanie 2011