5. PODSTAWOWE SFORMUA
OWANIA METODY ELEMENTÓW
SKOC
CZONYCH W NAWI
ZANIU DO RÓWNAC

MECHANIKI KONTINUUM
5.1. Podstawowe równania liniowej sprężystości
Stosowany w niniejszym opracowaniu opis omawianych zagadnień odnosi się do
prostokątnego układu współrzędnych kartezjańskich [x, y, z] lub, gdy stosujemy zapis
wskaz
nikowy, do układu [x1, x2, x3,]. Chcąc zdefiniować podstawowy układ równań opisujący stan
naprężenia, stan odkształcenia i pole przemieszczeń układu, będziemy się posługiwać następującymi
oznaczeniami.
Niech stan naprężenia w nieskończenie małej objętości ciała poddanego działaniu
obciążenia będzie opisany w układzie współrzędnych za pomocą składowych tensora,
uporządkowanych w macierzy à w postaci:
ij
Ã11 Ã12 Ã13
îÅ‚
ïÅ‚Ã
à = à à , (5.1)
ij 21 22 23
ïÅ‚
ïÅ‚Ã Ã Ã
ðÅ‚ 31 32 33
gdzie skÅ‚adowe Ã11 , à , à sÄ… naprężeniami normalnymi, natomiast Ã12 , Ã13 , Ã
22 33 23
opisują naprężenia styczne. Tensor stanu naprężenia (5.1) jest symetryczny, to
znaczy, że zachodzą następujące równości:
Ã12 = Ã , Ã13 = Ã , Ã = Ã . (5.2)
21 31 23 32
Stosując konsekwentnie zapis macierzowy, wygodnie jest niekiedy opisać stan naprężenia za
pomocÄ… wektora naprężenia à o nastÄ™pujÄ…cych skÅ‚adowych:
T
à = [à , à , à , à , à , à ] . (5.3)
xx yy zz xy xz yz
Tutaj, tak jak i poprzednio, za pomocą identycznych indeksów oznaczono składowe normalne, indeksy
zaÅ› różne, opisujÄ…ce skÅ‚adowe macierzy à informujÄ… o skÅ‚adowych stycznych stanu naprężenia.
Stan odkształcenia, podobnie jak poprzednio nawiązujący do opisu tensorowego,
reprezentuje macierz skÅ‚adowych µij w postaci:
42 Podstawowe sformułowania metody elementów skończonych...
µ11` µ12 µ13
îÅ‚
ïÅ‚µ
µij = µ22 µ23 . (5.4)
21
ïÅ‚
ïÅ‚µ31 µ32 µ33
ðÅ‚
W zapisie macierzowym posÅ‚ugiwać siÄ™ bÄ™dziemy wektorem odksztaÅ‚cenia µ , którego skÅ‚adowe
odniesione do układu [x, y, z], są następujące:
T
µ = [µ , µ , µ , Å‚ , Å‚ , Å‚ ] (5.6)
xx yy zz xy xz yz
Zwracamy w tym miejscu uwagę, że we wzorze (5.5) posługujemy się tzw. inżynierskimi definicjami
odkształceń stycznych, związanymi z odpowiednimi składowymi tensora odkształceń za pomocą
związków:
Å‚ = 2 Å"µ , Å‚ = 2 Å"µ , Å‚ = 2 Å"µ . (5.6)
xy xy yz yz xz xz
Przyjęcie w zapisie macierzowym miar inżynierskich odkształceń (ł - kąt odkształcenia postaciowego)
podyktowane jest dwoma faktami. Pierwszy to ich powszechne używanie w klasycznych zagadnieniach
liniowej sprężystości. Drugi zaś wynika z konieczności spójnego potraktowania miar naprężeń i odkształceń,
by w prosty sposób można było zapisać wyrażenie na pracę w obu zapisach - wskaz
nikowym i
macierzowym:
T
à Å"µij = à Å"µ . (5.7)
ij
Oprócz pól naprężeń i odkształceń do zapisania podstawowego układu równań konieczne jest jeszcze
pole przemieszczeń, którego składowe w punkcie opisane są w zapisie wskaz
nikowym:
T
ui = [u1, u2, u3] , (5.8)
lub w macierzowym:
T
u = [u, v, w] . (5.9)
5.1.1. Podstawowe równania w zapisie wskaz
nikowym
Dla przypomnienia zapiszmy podstawowy układ równań liniowej teorii sprężystości. Typowa
analiza ciaÅ‚a odksztaÅ‚calnego wymaga znalezienia funkcji naprężeÅ„ à lub przemieszczeÅ„ u
spełniających następujące równania: - trzy równania różniczkowe cząstkowe równowagi (równania Naviera)
Rozdział 5 43
à + bi = 0 i, j = 1,2,3, (5.10)
ij
gdzie w zapisie wskaz
nikowym zastosowano umowę sumacyjną, co znaczy, że powtarzający się w
jednomianie wskaz
nik informuje o konieczności dokonania sumowania po wszystkich możliwych jego
wartościach, a występujący między wskaz
nikami znak przecinka jest symbolem różniczkowania względem
"Ã
21
odpowiedniej zmiennej przestrzennej; na przykład à = ;
21,1
"x1
- sześć równań różniczkowych cząstkowych geometrycznych (równania Cauchy'ego)
1
µij = Å"(ui, j + u ), (5.11)
j,i
2
- sześć równań algebraicznych fizycznych (równania Hooke'a)
à = Eijkl Å"µkl , (5.12)
ij
Z powyższego zapisu nie wynika deklarowana uprzednio liczba równań, ale biorąc pod uwagę
założenia o izotropii, układ (5.12) redukuje się tylko do sześciu niezależnych równań i
występujących w nich tylko dwóch stałych materiałowych. Ponadto poszukiwane rozwiązania muszą
dodatkowo spełniać:
- równania nierozdzielności geometrycznej
µij,kl + µkl ,ij - µik , jl - µ = 0 (5.13)
jl,ik
w każdym punkcie obszaru, oraz
- naprężeniowe i przemieszczeniowe warunki brzegowe
à Å" n = pi* na brzegu Sà , (5.14)
ij j
ui = ui* na brzegu Su , (5.15)
przy czym oba deklarowane brzegi Sà i Su są rozłączne, tworząc w sumie cały brzeg
rozpatrywanego obszaru, tzn. SÃ )" Su = 0 oraz SÃ *" Su = S .
Domyślamy się, że ze względu na złożoność wymagań nakładanych na rozwiązania problemów,
określenie funkcji analitycznych, spełniających równania (5.10) - (5.15), nie jest łatwe. W szczególności
zadanie staje się niemożliwe do rozwiązania, jeśli skomplikuje się warunki brzegowe problemu. Zauważmy
jeszcze, że rozwiązanie problemu mechanicznego, opisanego za pomocą tak skonstruowanego modelu
matematycznego, prowadzi do zadania analizy matematycznej.
Trudności, na które natrafia się przy takim sformułowaniu, skłaniają do poszukiwania innych rozwiązań, tym
razem już nie analitycznych lecz rozwiązań przybliżonych.
44 Podstawowe sformułowania metody elementów skończonych...
5.1.2. Podstawowe równania w zapisie m acierzowym
Rozpocznijmy tym razem od równań geometrycznych. Odpowiednie składowe wektora
odkształceń można zapisać w następującej postaci:
"u "v "w "u "v
µ = , µ = , µz = , Å‚ = + ,
x y xy
"x "y "z "y "x
(5.16)
"v "w "u "w
Å‚ = + , Å‚ = + .
yz xz
"z "z "z "x
Używając poprzednio wprowadzonych oznaczeń, zapiszemy powyższy układ zależności w postaci:
µ = L Å" u , (5.17)
gdzie macierz operatorów różniczkowych L ma wymiar (6x3) a jej składowe można przedstawić
jako
" "x 0 0
îÅ‚
ïÅ‚
0 " "y 0
ïÅ‚
ïÅ‚ 0 0 " "z
L = . (5.18)
ïÅ‚" "y " "x 0
ïÅ‚
ïÅ‚
0 " "z " "y
ïÅ‚
ðÅ‚" "z 0 " "x
Równania równowagi można teraz zapisać krótko:
LT Å"Ã + b = 0 , (5.19)
gdzie b jest wektorem sił masowych. Zwróćmy uwagę na fakt, że macierz operatorów różniczkowych
równań równowagi (5.19) jest transponowana do odpowiedniej macierzy związków geometrycznych.
Równania fizyczne (konstytutywne), jako zależności między składowymi wektorów naprężeń i
odkształceń, określone są następująco:
à -½ Å"à -½ Å"à Ã
x y z xy
µ = , Å‚ = ,
x xy
E G
à -½ Å"à -½ Å"à Ã
y x z yz
µ = , Å‚ = , (5.20)
y yz
E G
à -½ Å"à -½ Å"Ã
Ã
z x y
zx
µz = , Å‚ = ,
zx
E G
Rozdział 5 45
gdzie przez E oznaczono moduÅ‚ odksztaÅ‚calnoÅ›ci podÅ‚użnej (moduÅ‚ Younga), zaÅ› G = E 2 Å" (1 +½ ) jest
moduÅ‚em odksztaÅ‚calnoÅ›ci postaciowej (moduÅ‚ Kirchhoffa), ½ jest liczbÄ… Poissona. W postaci równania
macierzowego powyższą zależność konstytutywną można wyrazić jako
µ = C Å"Ã , (5.21)
gdzie
1
îÅ‚ -½ -½ 0 0 0
ïÅ‚-½ 1 -½ 0
0 0
ïÅ‚
ïÅ‚-½ -½ 1 0 0 0
1
C = Å" . (5.22)
ïÅ‚
E 0 0 0 2 Å" (1 +½ ) 0 0
ïÅ‚
ïÅ‚
0 0 0 0 2 Å" (1 +½ ) 0
ïÅ‚
0 0 0 0 0 2 Å" (1 +½ )
ðÅ‚
Zależność (5.21) jest jednoznaczna, a kwadratowa macierz konstytutywna jest nieosobliwa, istnieje
więc odwzorowanie odwrotne
à = D Å"µ , (5.23)
gdzie macierz D = C-1 i jej reprezentacja przedstawia się następująco:
1
îÅ‚ -½ ½ ½ 0 0 0
ïÅ‚
½ 1 -½ ½ 0 0 0
ïÅ‚
ïÅ‚ ½ ½ 1 -½ 0 0 0
E
D = Å" . (5.24)
ïÅ‚
(1 +½ ) Å" (1 - 2 Å"½ ) 0 0 0 (1 - 2 Å"½ ) / 2 0 0
ïÅ‚
ïÅ‚
0 0 0 0 (1 - 2 Å"½ ) / 2 0
ïÅ‚
0 0 0 0 0 (1 - 2 Å"½ ) / 2
ðÅ‚
Podsumowując łatwo zauważyć i docenić zwięzłość stosowanego zapisu macierzowego, który pozwala
widzieć podstawowy układ równań w następującej postaci:
µ = L Å" u
LT Å"Ã + b = 0
(5.25)
µ = C Å"à lub à = D Å"µ
46 Podstawowe sformułowania metody elementów skończonych...
5.2. Analiza przybliżona problem u brzegowego
Jak wspomniano wyżej, możliwość znalezienia rozwiązań problemów brzegowych w postaci zamkniętych
formuł analitycznych ogranicza się, niestety, do wąskiej klasy zadań. W większości przypadków ważnych z
inżynierskiego punktu widzenia, to znaczy dla przypadków znajdujących zastosowania praktyczne, skom-
plikowane warunki podparcia układów, nietypowe obciążenia czy inne nieregularności uniemożliwiają
otrzymanie rozwiązań analitycznych. Chęć otrzymania wartościowych jakościowo i ilościowo wyników
opisujących stan układów zmusza do szukania odpowiedzi na drodze dyskretyzacji. Zamiast więc szukać
odpowiedzi układu w postaci pól naprężeń, odkształceń i przemieszczeń, poszukuje się wartości tych pól w
skończonej liczbie punktów należących do obszaru i jego brzegu. Z punktu widzenia zastosowań aparatu
matematycznego w przypadku stosowania dyskretyzacji uwalniamy siÄ™ od rozwiÄ…zywania problemu
różniczkowego, zastępując go zadaniem algebraicznym. Nie chcemy w tym miejscu dyskutować o różnych
możliwościach stosowania dyskretyzacji, a co za tym idzie, o różnych metodach rozwiązywania problemów
brzegowych. Podkreślamy tylko, że omawiana tutaj MES zakłada analizę przybliżoną, polegającą na
podziale całego układu na mniejsze części (elementy), posiadające charakterystyczne punkty zwane
węzłami, w których to punktach skoncentrowana jest niejako pełna informacja o zachowaniu się tych
elementów i ich własnościach. W spomniane przybliżenie polega - w najbardziej podstawowej wersji - na
przyjęciu pola przemieszczeń opisującego przemieszczenie dowolnego punktu elementu, jako funkcji
przemieszczeń węzłów i położenia danego punktu (jest to tzw. wersja przemieszczeniowa MES). Niewiado-
me są więc przemieszczenia węzłów. Musimy być świadomi, że przyjmowane funkcje określające pole
przemieszczeń elementu zwykle nie odpowiadają w pełni funkcjom analitycznym rozwiązującym problem
różniczkowy. Innymi słowy, popełniamy na tym etapie błędy, które, jak można to udowodnić, maleją w miarę
jak rośnie liczba elementów, na które podzielono cały układ. Musimy być także świadomi, że przyjmując pole
przemieszczeń w postaci określonych funkcji, deklarujemy tym samym przez związki geometryczne pole
odkształceń i dalej przez zależności konstytutywne - pole naprężeń. Jeśli w określonych przypadkach
szczególnie zależy nam na w miarę jak najlepszym odwzorowaniu pola odkształceń bądz
naprężeń, istnieją
inne możliwości przyjęcia funkcji aproksymacyjnych, zakładających wprost te właśnie pola. Takie
sformułowania MES nie będą jednak przedmiotem niniejszego opracowania.
Rozważana wersja przemieszczeniowa MES w celu przeanalizowania problemu brzegowego
wymaga podjęcia następujących kroków:
- dokonania podziału układu (konstrukcji, kontinuum) na skończoną liczbę pod-obszarów o
prostej geometrii,
- wybrania punktów węzłowych (węzłów), w których zostaną zapewnione warunki równowagi i
zgodności przemieszczeń,
- założenia funkcji przemieszczeń w obszarach każdego elementu, takiego że
przemieszczenia wszystkich punktów zależą od przemieszczeń węzłów,
- speÅ‚nienia w elemencie zależnoÅ›ci µ = L Å" u oraz à = D Å"µ ,
- wyznaczenia sztywności elementów i równoważnych sił węzłowych,
- zbudowania układu równań równowagi dla węzłów zdyskretyzowanego kontinuum,
- rozwiązania układu równań równowagi dla przemieszczeń w ęzłów,
- obliczenia przemieszczeń, odkształceń i napręże ń w wybranych punktach elementów,
Rozdział 5 47
- obliczenia reakcji podpór.
5.3. Podstawy MES wynikające z równania pracy wirtualnej
Załóżmy, że trójwymiarowy element skończony jest zdefiniowany w kartezjańskim układzie współrzędnych
[x, y, z,]. Niech wektor u , opisujący przemieszczenie dowolnego punktu elementu, jest wyrażony za
T
pomocą składowych: u = [u, v, w] , gdzie u, v, w są - odpowiednio - przemieszczeniami w kierunku osi
T
x, y, z . Siły masowe oznaczymy za pomocą wektora b = [bx , by , bz] , gdzie składowe oznaczają
siły przypadające na jednostkę objętości, powierzchni lub długości. Przez d oznaczymy wektor
przemieszczeń węzłowych elementu.. W ymiar tego wektora jest równy liczbie węzłów elementu pomnożonej
przez liczbę przyjętych stopni swobody węzła. Jeśli założymy, że przemieszczenia węzła opisują składowe
przesunięć w kierunku osi x, y, z oraz jeśli n jest liczbą węzłów w elemencie, to
d = [di], i = 1, 2,...,nen ,
(5.26)
gdzie
di = [dxi , d , dzi].
yi
(5.27)
Zauważmy tylko, że inne typy przemieszczeń, takie jak obroty czy krzywizny, mogą również być
i będą dalej traktowane jako składowe wektora przemieszczeń. Podobnie przyjmijmy siły węzłowe p
jako składowe sił we wszystkich węzłach elementu w kierunkach osi x, y i z :
p = [pi], i = 1,2,...,nen , (5.28)
gdzie
pi = [pxi , pyi , pzi]. (5.29)
Po tych definicjach wstępnych załóżmy pole przemieszczeń w elemencie jako funkcję przemieszczeń
węzłów elementu w postaci:
u = N Å" d . (5.30)
Ponieważ wektor u ma wymiary (3x1) , zaś wektor przemieszczeń węzłów d wymiary liczby stopni
swobody elementu nedf = nen x 3, więc macierz funkcji próbnych, inaczej zwanych funkcjami kształtu, jest
macierzą prostokątną o wymiarach 3 x nedf . Każda ze składowych macierzy N jest funkcją i określa wpływ
48 Podstawowe sformułowania metody elementów skończonych...
danej składowej wektora przemieszczeń d na przemieszczenie dowolnego punktu elementu o
współrzędnych x, y, z .
Zależność µ(u) otrzymuje siÄ™ przez różniczkowanie stosownych wyrażeÅ„ na przemieszczenia,
µ = L Å" u µ = L Å" N Å" d µ = B Å" d . (5.31)
Macierz B opisuje więc odkształcenia w każdym punkcie elementu, spowodowane jednostkowym
przemieszczeniem kolejnych stopni swobody węzłów. Z prawa fizycznego łatwo więc wyprowadzić, że
à = D Å"µ à = D Å" B Å" d , (5.32)
gdzie iloczyn macierzy D Å" B opisuje - podobnie jak poprzednio - zmiany naprężeÅ„ jako funkcje przemieszczeÅ„
węzłów.
Zasada prac wirtualnych głosi, że jeśli układ znajdujący się w równowadze poddany jest wirtualnym
przemieszczeniom (kinematycznie zgodnym stanom deformacji), wówczas praca wirtualna zewnętrznych
obciążeń jest równa wirtualnej energii odkształcenia naprężeń wewnętrznych:
´ Ue = ´ We , (5.33)
gdzie U jest wewnÄ™trznÄ… energiÄ… odksztaÅ‚cenia, W - pracÄ… siÅ‚ zewnÄ™trznych, zaÅ› ´ oznacza
wariację (stan wirtualny - pomyślany, zgodny z więzami).
W prowadz
my wiÄ™c wirtualny stan przemieszczeÅ„ wÄ™zÅ‚owych i oznaczmy go przez´d = [´di]
(i = 1, 2,...,nen ) . W irtualne przemieszczenia i odkształcenia można wówczas wyrazić jako
´ u = N Å"´ d oraz ´ µ = B Å"´ d . (5.34)
W irtualna energia układu i praca wirtualna sił zewnętrznych wyrażają się teraz w postaci wzorów:
T
´ Ue = ´µ Å"Ã Å" dV i ´ We = ´ pT Å" p + ´uT Å" b Å" dV (5.35)
V V
kolejno podstawiajÄ…c otrzymujemy z (5.33)
T
´µ Å"Ã Å" dV = ´ pT Å" p + ´uT Å" b Å" dV . (5.36)
V V
Rozdział 5 49
Uwzględniając (5.32) i (5.34), otrzymujemy
T T T T
´ d BT Å" D Å"µ Å" dV = ´ d Å" p +´ d N Å" b Å" dV , (5.37)
V V
T
i raz jeszcze wykorzystujÄ…c (5.34) i upraszczajÄ…c przez ´ d , otrzymujemy
ëÅ‚
T
ìÅ‚
BT Å" D Å" B Å" dV Å" d = p + N Å" b Å" dV (5.38)
ìÅ‚
íÅ‚V V
lub ostatecznie
K Å" d = p + pb , (5.39)
gdzie K jest tzw. macierzą sztywności elementu, której składowe mogą być interpretowane jako
fikcyjne siły w węzłach, spowodowane jednostkowymi ich przemieszczeniami, pb zawiera
równoważne siły węzłowe spowodowane masą ciała.
µ0
Obecność początkowego stanu odkształceń można uwzględnić w następujący sposób:
- założyć superpozycję stanów odkształceń
µ = µ0 + C Å"Ã , (5.40)
- stąd naprężenie
à = D Å" (µ - µ0 ) (5.41)
i po podobnych podstawieniach i przekształceniach, jak to uczyniono powyżej, otrzymujemy:
K Å" d = p + pb + p0 , (5.42)
gdzie p0 = BT Å" D Å"µ0 Å" dV jest wektorem równoważnym obciążeniom wÄ™złów od poczÄ…tkowego stanu
V
odkształceń (np. wpływ temperatury). Czytelnik mógłby dla nabrania umiejętności sprawdzić poprawność
wyprowadzonego wzoru (5.42).
50 Podstawowe sformułowania metody elementów skończonych...
Prześledz
my na prostym przykładzie pręta (rys. 5.1) postacie opisywanych macierzy oraz sposób
dojścia do sformułowania macierzy sztywności prostego elementu. Podkreślmy jednak w tym miejscu, że jak
dotąd próbujemy wyłącznie zdefiniować składowe stosownych macierzy i wektorów, odniesione do lokalnego
układu współrzędnych i tylko do jednego elementu.
Rys. 5.1. Funkcje kształtu dla dwuwęzłowego elementu kratownicy
W ektor przemieszczenia upraszcza się tutaj do jednej tylko składowej u = [ux], podobnie zresztą jak wektor
sil masowych b = [bx]. Element jest dwu węzłowy i ma po jednym stopniu swobody w każdym węz
le, tak
więc globalny wektor przemieszczeń jest tylko dwuelementowy d = [d1, d2]= [u1, u2]. Podobnie rzecz się
ma z obciążeniami p = [p1, p2]= [px1, px2]. Przyjmijmy funkcję przemieszczeń w postaci liniowej:
u = c1 + c2 Å" x , (5.43)
gdzie stałe ci wyznaczymy z warunków brzegowych
dla x = 0 u = d1 c1 = d1,
(5.44)
dla x = L u = d2 c2 = (d2 - d1)/ L .
Przemieszczenie dowolnego punktu wyraża się zatem wzorem
d1
x x îÅ‚
îÅ‚1
u = - , Å" = N Å" d , (5.45)
ïÅ‚
L L ïÅ‚d2
ðÅ‚
ðÅ‚
Rozdział 5 51
gdzie macierz funkcji kształtu N składa się z dwóch funkcji liniowych. Przebieg tych funkcji
zilustrowano na rysunku 5.1.
Odkształcenia dla tego prostego przypadku opisano tylko jedną składową
du dN
µ = [µ ]= L Å" u = = Å" d = B Å" d (5.46)
x
dx dx
więc
1
B = Nx = Å"[-1, 1]. (5.47)
L
Stan naprężenia również sprowadza się do jednej tylko składowej:
à = à = D Å"µ = E Å"µ = E Å" B Å" d , (5.48)
x x
gdzie operator konstytutywny uprościł się do jednej tylko stałej.
Rys. 5.2. Obciążenia pręta kratownicy: osiowe ciągłe i liniowo zm ienne
Macierz elementu otrzymujemy teraz z podstawienia:
L
E Å" A îÅ‚ -1
E îÅ‚-1 1
K = BT Å" D Å" B Å" dV = Å" Å"[-1 1] dA dx = (5.49)
ïÅ‚-1 1
1
L2 ïÅ‚ L
ðÅ‚ ðÅ‚
V 0 A
gdzie przyjęto, że pole powierzchni przekroju pręta jest stale na całej jego długości.
Przyjmijmy dodatkowo, że pręt obciążony jest siłą masową, zmieniającą się liniowo, tak jak to
pokazano na rysunku 5.2 według funkcji:
52 Podstawowe sformułowania metody elementów skończonych...
b2 - b1
bx = b1 + Å" x ; (5.50)
L
wówczas wektor sił masowych działających w węzłach wynosi:
L
2 Å" b1 + b2
1 îÅ‚
T
pb = N Å" bx Å" dx = Å" . (5.51)
ïÅ‚b + 2 Å" b2
6
ðÅ‚ 1
0
Gdy element poddany jest działaniu temperatury "T , mamy do czynienia z początkowymi odkształceniami
µ0 = µT = Ä… Å" ("T ) , gdzie przez Ä… oznaczono współczynnik rozszerzalnoÅ›ci cieplnej materiaÅ‚u. SiÅ‚y
przykładane w węzłach elementu, spowodowane początkowymi odkształceniami, wynoszą:
L
îÅ‚-1
p0 = pT = BT Å" D Å"Ä… Å" ("T )dA dx = E Å" AÅ"Ä… Å" ("T ) Å" . (5.52)
ïÅ‚
1
ðÅ‚
0 A
Czytelnik mógłby zadać sobie trud sprawdzenia poprawności wyników wzorów (5.51) i (5.52).
5.4. Podstawy MES wyprowadzone z twierdzenia o m inim um całkowitej
energii potencjalnej
Otrzymane w poprzednim rozdziale równania MES uzyskuje się również przez zastosowanie twierdzenia o
minimum całkowitej energii potencjalnej. Twierdzenie to głosi, że spośród wszystkich kinematycznie
dopuszczalnych pól przemieszczeń spełnia się to, które całkowitej energii potencjalnej zapewnia minimum.
Całkowita energia potencjalna układu wyraża się jako:
 = U -W , (5.53)
gdzie U oznacza energię sprężystą ciała, a W jest pracą sił zewnętrznych. A
atwo wykazać, że  dla ciała
liniowo-sprężystego jest funkcjonałem kwadratowym i ma jedno globalne minimum. Rozwiązanie jest więc
jednoznaczne. Energię  zapiszemy więc w postaci:
1
T
 = Ã Å"µ Å" dV - uT Å" b Å" dV - uT Å" p* Å" dS , (5.54)
2
V V S
gdzie p* jest danym obciążeniem brzegu S . Przyjmując interpolację dla elementu w znanym już nam
kształcie:
T
u = N Å" d, µ = B Å" d, Ã = D Å"µ,
możemy powyższe twierdzenie ograniczyć do obszaru elementu i zapisać:
Rozdział 5 53
1
T T T T T
 = d Å" BT Å" D Å" B Å" d Å" dV - d Å" N Å" b Å" dV - d Å" N Å" p* Å" dS ; (5.55)
2
Ve Ve SÃ
ze stacjonarności tego wyrażenia wynika równowaga elementu:
" e
= 0 BT Å" D Å" B Å" dV Å" d - pe = K Å" d - pe , (5.56)gdzie pe opisuje siÅ‚y
"d
Ve
węzłowe danego elementu jako efekt obciążeń masowych b i powierzchniowych p* :
T T
pe = N Å" b Å" dV + N Å" p* Å" dS . (5.57)
Ve SÃ
Sprawę modyfikacji wyprowadzonych wzorów dla przypadków uwzględniających udział odkształceń
wstępnych pozostawia się Czytelnikowi.
Energia sprężysta pojedynczego elementu belkowego bez uwzględnienia wpływu ścinania wynosi:
1 1 E
2
T T
Ue = Ã Å"µ Å" dV = µ Å" D Å"µ Å" dV = µ Å" dV . (5.58)
x
2 2 2
Ve Ve Ve
Jeśli przyjmiemy klasyczne założenie belki Bernoulli'ego, że odkształcenie jest funkcją
przemieszczenia (ugięcia),
2
d v
µ = - y Å" , (5.59)
x
dx2
wówczas energia wewnętrzna elementu zginanego wynosi:
2 2 2 2
1 1
2 2 2 2
E ëÅ‚ E ëÅ‚ ëÅ‚ E Å" I d v
d v d v E d v ëÅ‚
ìÅ‚ ìÅ‚
Ue = ìÅ‚- y Å" Å" dV = y2 Å"ìÅ‚ 2 Å" dV =
2 2
ìÅ‚ ìÅ‚dx dA dx = 2 ìÅ‚ dx . (5.60)
ìÅ‚dx
2 dx2 2 2
íÅ‚ íÅ‚dx íÅ‚ íÅ‚
Ve Ve 0 A 0
Równanie powyższe interpretuje energię wewnętrzną elementu belkowego jako funkcję przemieszczeń,
v = f (x)
.
54 Podstawowe sformułowania metody elementów skończonych...
Rys. 5.3. Postacie funkcji kształtu dla elem entu belkowego
Rozpatrzmy element belkowy, płaski, dwuwęzłowy zginany w płaszczyz
nie x0 y , jak na rysunku 5.3. W ektor
przemieszczeń węzłowych przyjmijmy w postaci d = [d1, d2, d3, d4]= [v1, Ć1, v2, Ć2], gdzie przez v
oznaczono przemieszczenie prostopadłe do osi pręta, zaś Ć jest kątem obrotu przekroju. Indeksy 1,2
odnoszą się do numeracji węzłów. Zgodnie z założeniami klasycznej teorii belek kąty obrotu są pochodnymi
przemieszczeń:
Rozdział 5 55
dv1 dv2
Ć1 = , Ć2 = . (5.61)
dx dx
Odpowiedni wektor sił węzłowych p = [p1, m1, p2, m2] zawiera siły skupione działające w kierunku
przemieszczeń oraz momenty zginające zgodne z kątami obrotów przekrojów.
Załóżmy funkcję przemieszczeń w postaci kompletnego wielomianu trzeciego stopnia:
v(x) = c1 + c2 Å" x + c3 Å" x2 + c4 Å" x3 . (5.62)
W funkcji tej współczynniki ci wyznaczymy z warunków brzegowych, które definiują wielkości
przemieszczeń i kątów obrotów na końcach elementów jako równe składowym wektora d. Zapiszmy te
warunki:
dv(0)
dla x = 0 v(0) = v1 oraz = Ć1
dx
(5.63)
dv(l)
dla x = l v(l) = v2 oraz = Ć2
dx
Po wyznaczeniu stałych ci zapiszemy macierz funkcji kształtu:
1
N = Å"[2 Å" x3 - 3Å" l Å" x2 + l3, l Å" x3 - 2 Å" l2 Å" x2 + x Å" l3, - 2 Å" x3 + 3Å" l Å" x2, l Å" x3 - l2 Å" x2]. (5.64)
l3
Funkcje kształtu, które przedstawiono na rysunku 5.3, opisują zmianę przemieszczenia v(x) , spowodowaną
jednostkowymi przemieszczeniami węzłów. Jeśli założymy dalej prawdziwość hipotezy płaskich przekrojów,
wówczas przemieszczenie podłużne wyniesie:
dv
u(x) = - y Å" , (5.65)
dx
skąd odkształcenie
2 2
du d v d v
µ = = - y Å" = - y Å"º , gdzie º = . (5.66)
x
dx dx2 dx2
W idzimy zatem, że operator różniczkowy L , transformujący przemieszczenie v(x) w odkształcenie
µ , ma postać:
x
56 Podstawowe sformułowania metody elementów skończonych...
2
d
L = - y Å" , (5.67)
dx2
skÄ…d macierz B = L Å" N otrzymujemy w postaci:
y
B = L Å" N = - Å"[12 Å" x - 61, 61Å" x - 412, -12 Å" x + 61, 61Å" x - 212]. (5.68)
l3
PamiÄ™tajÄ…c, że dla tego prostego przypadku zwiÄ…zek fizyczny ma postać à = E Å"µ (czyli operator
x x
D = E ), otrzymujemy macierz sztywności K dla elementu belkowego:
16 6 Å" l -12 6 Å" l
îÅ‚
E Å" Iz ïÅ‚ 6 Å" l 4 Å" l2 - 6 Å" l 2 Å" l2
ïÅ‚
Ke = Å" , (5.69)
l3 ïÅ‚-12 - 6 Å" l 12 - 6 Å" l
ïÅ‚
6 Å" l 2 Å" l2 - 6 Å" l 4 Å" l2
ðÅ‚
gdzie Iz = y2dA . Drobne przekształcenia, które należało wykonać, by w końcu otrzymać jawną postać
A
macierzy Ke , pozostawiamy Czytelnikowi.
Równoważne obciążenia węzłowe, wynikające z przyjęcia ciężaru równomiernie rozłożonego by
bÄ…dz
liniowo zmieniajÄ…cego siÄ™ by (x) , jak na rysunku 5.4 wynoszÄ… odpowiednio:
l
by Å" l
T
pb = N Å" by Å" dx = Å"[6, l, 6, - l] - obciążenie staÅ‚e,
12
0
by Å" l
pb = Å"[9, 2 Å" l, 21, - 3Å" l] - obciążenie liniowo zmienne (5.70)
60
Rys. 5.4 Obciążenia elem entu belkowego: obciążenie równomierne i liniowo zm ienne
Rozdział 5 57
Chcąc uwzględnić również wpływ odkształceń początkowych załóżmy, że element poddany jest liniowej
zmianie temperatury od "T1 na części dolnej do "T2 na górnej. Jeśli "T1 > "T2 i wysokość elementu jest
równa h , to zmiana temperatury w każdym punkcie wynosi:
1 y
"T = Å"("T1 + "T2)- Å"("T1 - "T2). (5.71)
2 h
Pierwszy człon opisuje efekt równomiernego ogrzania, a ponieważ nie wywołuje zginania, zostanie w
dalszych rozważaniach pominięty. Drugi człon powoduje odkształcenia od zginania:
y
µ = -Ä… Å" Å"("T1 - "T2 ),
(5.72)
xT
h
znajdujemy więc
pT = BT Å"Ã Å" dV =
T
V
(5.73)
l
y Ä… Å" E Å" Iz
= - BT Å" E Å"Ä… Å" Å"("T1 - "T2)Å" dA dx = Å"("T1 - "T2)Å"[0, - l, 0, l]
h h
0 A
Jeżeli znane sÄ… wyrażenia okreÅ›lajÄ…ce krzywizny poczÄ…tkowe º0 elementu, to poczÄ…tkowe
odkształcenia wyrażają się zależnością:
µ = - y Å"º0 , (5.74)
x0
siły węzłowe wyznaczymy zgodnie z (5.42) jako:
pT = BT Å" D Å"µ0 Å" dV . (5.75)
V
5.5. Podsum owanie
Spróbujmy na koniec tego rozdziału uświadomić sobie, w jaki czysto formalny sposób możemy
zbudować macierze sztywności elementów oraz wektory obciążeń, wynikające bądz
z działania sił
masowych, bÄ…dz
z wstępnych odkształceń. Zapamiętajmy następujący tok postępowania:
1. Rozpoczynamy od aproksymacji pola przemieszczeń, którą można wyrazić następująco :
u = g Å" c , (5.76)
58 Podstawowe sformułowania metody elementów skończonych...
gdzie przez g oznaczyliśmy tzw. macierz geometryczną, która najczęściej gromadzi odpowiednie potęgi
stosowanych wielomianów interpolacyjnych, zaś c jest macierzą stałych. Stałe te wyznaczymy z warunków
brzegowych, ( przemieszczenia w węzłach muszą być zgodne z wartościami przemieszczeń, wynikającymi z
przyjętych funkcji )
2. W arunki brzegowe wyrażą się w postaci:
d = h Å" c , gdzie h = [gi] dla i = 1,2,...,nedf . (5.77)
Macierz h jest macierzą kwadratową i nieosobliwą, tak więc z układu równań (5.77) można wyznaczyć stałe
wielomianów interpolacyjnych jako funkcji przemieszczeń węzłów,
c = h-1 Å" d . (5.78)
3. Funkcje kształtu otrzymamy teraz automatycznie i formalnie:
u = g Å" h-1 Å" d = N Å" d , (5.79)
wiÄ™c N = g Å" h-1 .
4. ZnajÄ…c postać operatora różniczkowego L , z Å‚atwoÅ›ciÄ… wyznaczymy macierz B = L Å" N .
5. Teraz zupeÅ‚nie formalnie przy ustalonym prawie konstytutywnym à = D Å"µ otrzymujemy
macierz sztywności K oraz pozostałe wektory pb, p0 lub pT .
Zaproponowany sposób postępowania spróbujemy wykorzystać w dalszych rozważaniach. Należy
jednak zaznaczyć, że o ile dla elementów o niewielkiej liczbie stopni swobody taki formalny sposób podejścia
jest wygodny, o tyle dla elementów bardziej skomplikowanych może okazać się nieskuteczny. W takich
przypadkach wygodniej będzie od razu próbować zdefiniować postacie funkcji kształtu N, a nie uzyskiwać
ich w sposób formalny. Ma to miejsce głównie w sytuacjach, gdy unika się budowania jawnej postaci
macierzy sztywności elementu, a otrzymuje się ją w wyniku zabiegów numerycznych.
Na koniec rozważań na temat formułowania elementów skończonych podejmijmy próbę odpowiedzi
na pytanie, kiedy rozwiązanie równania różniczkowego, opisującego dane zagadnienie brzegowo-
początkowe otrzymane za pomocą MES będzie zbiegać się z rozwiązaniem analitycznym (dokładnym). Czy
w miarę zwiększania liczby elementów skończonych rozwiązanie to będzie zbieżne z rozwiązaniem do-
kładnym? Zaznaczmy przed rozpatrzeniem tego problemu, że rozwiązania otrzymywane metodą elementów
skończonych są obarczone kilkoma typami błędów, wynikającymi z: błędów zaokrągleń obliczeń
komputerowych, błędów wynikających z aproksymacji praw konstytutywnych, błędów powstałych z
Rozdział 5 59
całkowania macierzy i błędów metod rozwiązywania równań (sposobu całkowania równań ruchu). Poniżej
rozpatrzymy tylko błędy wynikające z dyskretyzacji, czyli idealizacji konstrukcji czy kontinuum materialnego
elementami skończonymi.
Można wykazać, że w celu zapewnienia monotonicznej zbieżności rozwiązań elementy skończone
muszą spełniać dwa zasadnicze kryteria: zupełności i zgodności. Jeżeli są one spełnione, to
dokładność wyników rośnie w miarę zagęszczania siatki podziału na elementy.
W arunek zupełności wymaga, by funkcje przemieszczeń elementu mogły reprezentować jego ruch sztywny
(beznaprężeniowy) oraz stan stałych odkształceń. Na przykład dla elementu płaskiego wymagane jest by
funkcje przemieszczeń mogły przedstawić 3 postacie ruchu sztywnego (dwa translacyjne i jeden sztywny ob-
rót) oraz stan stałego odkształcenia. Stan stałego odkształcenia można zinterpretować wymaganiem, by w
miarę zagęszczania elementów przez ich pomniejszanie, w elementach takich odkształcenia powinny być
stałe, by móc skolei reprezentować dowolne zmienny stan odkształcenia całego układu.
Drugie kryterium, tzw. kryterium zgodności elementu, oznacza, że przemieszczenia wewnątrz elementu jak i
na jego brzegach powinny być ciągłe. Chodzi o to, by nie pojawiały się nieciągłości pola przemieszczeń
pomiędzy elementami w sytuacji, gdy układ elementów zostanie poddany obciążeniu. W przypadku gdy ma-
my do czynienia tylko z translacyjnymi stopniami swobody wymaganie to sprowadza siÄ™ do sprawdzenia
tylko ciągłości przemieszczeń u, v i w . W przypadku występowania rotacyjnych stopni swobody
zdefiniowanych jako pochodne przemieszczeń (w elementach belkowych i płytowych), należy spełnić to
wymaganie również dla tych stopni swobody, czyli spełnić ciągłość ich pierwszych pochodnych. Ciągłość tę
zazwyczaj trudno jest spełnić dla elementów płytowych, w których kąty obrotu są otrzymywane przez
różniczkowanie przemieszczeń poprzecznych.
Elementy nie spełniające powyższych kryteriów nazywane są elementami niedostosowanymi a ich
stosowanie nie gwarantuje . monotonicznej zbieżności wyników.
To, czy element jest zgodny i zupełny, zależy od użytego sformułowania i każde sformułowanie należy
sprawdzić indywidualnie.
Zadania
1. Zapisz funkcjonały całkowitej energii potencjalnej w notacji wskaz
nikowej i macierzowej dla
ogólnego problemu kontinuum.
2. Podaj wyprowadzenie wzoru na macierz sztywności elementu belkowego z twierdzenia o minimum
całkowitej energii potencjalnej.
3. W yprowadz
wzory na postaci macierzy sztywności elementu prętowego (kratownicy) z równania
pracy wirtualnej i twierdzenia o minimum całkowitej energii potencjalnej. Jakie są wspólne cechy
tego wyprowadzenia?
4. Dana jest rama o geometrii przedstawionej poniżej na rysunku. Pozostałe dane o przekrojach
przyjmij z tablicy. Należy :
- ponumerować węzły i pręty,
- sformułować macierze połączeń węzłów,
60 Podstawowe sformułowania metody elementów skończonych...
- obliczyć macierze sztywności wybranych elementów w układzie globalnym,
- dokonać agregacji globalnej macierzy sztywności układu,
- zmodyfikować układ równań zgodnie z warunkami brzegowymi.
Rysunek geom etrii ram y płaskiej
a A A I I E P
[kN]
[m] [cm2 ] [cm2 ] [cm4 ] [cm4 ] [GPa]
1.2 40 22 3000 1000 200 10
2.0 18 28 570 1450 200 12
1.0 240 96 9000 1200 80 6
1.5 240 72 9000 900 100 3