SiM - z.258 - Pop PDF

Studia i Monografie z.
258
Lilianna Sadecka
Metoda różnic skończonych i metoda

elementów skończonych w zagadnieniach
mechaniki konstrukcji i podłoża
2010
SPIS TREŚCI
Przedmowa ………………………………………………………………. 5
1. Wstęp…………………………………………………………………… 7
2. Metoda różnic skończonych…………………………………………... 13
2.1. Wprowadzenie……………………………………………………... 13
2.2. Operatory różnicowe……………………………………………….. 14
2.3. Ocena dokładności operatorów różnicowych……………………… 20
2.4. Zastosowanie MRS w analizie statycznej belek…………………… 23
2.4.1. Przykład – belka obciążona siłą skupioną……………………. 23
2.4.2. Przykład – belka obciążona momentem skupionym…………. 28
2.5. Zastosowanie MRS w analizie statycznej płyt……………………... 32
2.5.1. Przykład – płyta obciążona równomiernie na całym obszarze 41
2.6. Wariacyjne ujęcie metody różnic skończonych (WMRS)…………. 45
2.6.1. Przykład – belka obciążona siłą skupioną…………………… 46
3. Metoda elementów skończonych……………………………………… 51
3.1. Wprowadzenie……………………………………………………... 51
3.2. Idea metody elementów skończonych……………………………... 52
3.3. Dyskretyzacja………………………………………………………. 55
3.4. Klasyfikacja elementów skończonych……………………………... 59
3.5. Równania metody elementów skończonych……………………….. 61
4 Spis treści
3.6. Funkcje kształtu……………………………………………………. 62

3.7. Agregacja…………………………………………………………... 72
3.8. Przykład – równanie różniczkowe podłoża Własowa…………..…. 81
3.9. Ogólne sformułowanie MES dla zagadnień mechaniki ciała stałego 88
4. Charakterystyka wybranych elementów skończonych……………... 95
4.1. Uwagi wstępne……………………………………………………... 95
4.2. Element belkowy…………………………………………………… 95
4.2.1. Przykład – belka z przegubem……………………………….. 106
4.3. Izoparametryczny element płyty średniej grubości wyższego rzędu 113
4.3.1. Związki podstawowe modelu płyty wyższego rzędu………... 114
4.3.2. Specyfikacja elementu płyty wyższego rzędu………………. 116
4.4. Elementy nieskończone…………………………………………….. 120
4.4.1. Zasadnicza koncepcja elementów nieskończonych…………. 121
4.4.2. Dwuwymiarowe elementy nieskończone……………………. 123
4.5. Makroelement podłoża warstwowego……………………………... 126
4.5.1. Związki podstawowe modelu podłoża warstwowego……….. 127
4.5.2. Specyfikacja makroelementu podłoża……………………..... 131
4.5.3. Uwagi dotyczące kalibracji parametrów modelu podłoża
warstwowego ……………………………………………… 138
5. Przykłady zastosowań MRS i MES w analizie statycznej
konstrukcji posadowionych na podłożu sprężystym ………………... 139
5.1. Uwagi wstępne……………………………………………………... 139
5.2. Belka spoczywająca swobodnie na podłożu Winklera – WMRS….. 142
5.3. Belka spoczywająca swobodnie na podłożu Własowa – WMRS….. 148
5.4. Belka na podłożu Winklera z więzami jednostronnymi – WMRS.... 153
5.5. Belka na podłożu dwuwarstwowym – MES……………………….. 156
5.6. Płyta wyższego rzędu na podłożu dwuwarstwowym – MES……… 165
5.6.1. Płyta obciążona równomiernie na całym obszarze………….. 171
5.6.2. Płyta obciążona centrycznie na małym obszarze……………. 174
5.7. Płyta fundamentowa na podłożu dowolnie uwarstwionym – MES .. 176
Literatura………………………….……………………………………… 183
PRZEDMOWA
Metoda elementów skończonych i metoda różnic skończonych stanowią

w obecnej dobie potężne narzędzie obliczeniowe, powszechnie wykorzystywane
w analizie wielu zagadnień techniki. Literatura omawiająca teorię obu metod
i ich aplikacji do konkretnych zagadnień jest bardzo obszerna. Spośród pozycji
książkowych warto tu wymienić pracę O.C. Zienkiewicza i R.L. Taylora
„Metoda elementów skończonych”, vol. 1, 2, (2000), pracę pod redakcją
M. Kleibera „Komputerowe metody mechaniki ciała stałego”, (1995) czy
książkę E. Rusińskiego, J. Czmochowskiego i T. Smolnickiego „Zaawansowana
metoda elementów skończonych w konstrukcjach nośnych” (2000). Ogólnie,
wśród pozycji książkowych, dają się zauważyć dwa trendy. Jeden obejmuje
prace omawiające wspomniane metody w ścisłym ujęciu matematycznym,
wymagające dobrej znajomości teorii równań różniczkowych, rachunku
wariacyjnego i analizy funkcjonalnej, drugi tworzą prace przedstawiające
szczegóły algorytmów metod na gruncie mechaniki konstrukcji. W niniejszej
książce, w której przedstawiono zarys teorii metody różnic skończonych
i metody elementów skończonych autorka starała się niejako te dwa główne
nurty zrównoważyć, dając ogólny opis metod, niekoniecznie związany
z dziedziną mechaniki konstrukcji, jakkolwiek przykłady ilustrujące, zawarte
w tej części książki, odnoszą się do zagadnień mechaniki budowli. Druga część
książki, chociaż nie wyodrębniona formalnie, ujmuje charakterystykę
wybranych elementów skończonych, zarówno tych znanych i powszechnie
stosowanych, takich, jak element belkowy, jak też szczególnych, takich jak
element płyty wyższego rzędu czy element nieskończony, dla których trudno
znaleźć nawet obszerniejsze omówienie w literaturze przedmiotu. Dotyczy to
w pełni sformułowanego i opisanego makroelementu podłoża sprężystego, który
stanowi własną autorską aplikację metody do zagadnień współpracy konstrukcji
6 Przedmowa
z podłożem dowolnie uwarstwionym. Zastosowaniu MRS i MES do tego

rodzaju zagadnień poświęcony jest ostatni rozdział książki. Zawarte w nim
przykłady, stopniowane ze względu na przyjętą idealizację masywu gruntowego
w zagadnieniach kontaktowych „konstrukcja – podłoże sprężyste”, obrazują
w konkretny sposób zastosowanie rozważanych metod w zagadnieniach
z pogranicza dwóch dziedzin: mechaniki konstrukcji i geomechaniki. W efekcie
ukazują one (ostatni przykład) możliwość analiz tych zagadnień przy
jednoczesnym, dostatecznie dokładnym modelowaniu obu substruktur: budowli
i podłoża uwarstwionego, a zarazem przy stosunkowo małym wysiłku
obliczeniowym. Jest to istotne, wciąż w powszechnym użyciu są bowiem
programy komputerowe do analizy konstrukcji na sprężystym podłożu opisanym
za pomocą najbardziej wyidealizowanego modelu-analogu Winklera. Leżące
u podstaw tego podejścia założenie, że osiadanie podłoża w danym punkcie
zależy tylko od wartości obciążenia powierzchniowego w tym punkcie ogranicza
możliwość uzyskiwania realistycznych wyników analiz do podłoża w postaci
cienkiej warstwy. Niezależnie od tego poważne trudności sprawia wiarygodne
oszacowanie współczynnika Winklera dla podłoża uwarstwionego. W tym
zakresie książka wypełnia pewną lukę między nadmiernie uproszczonym
traktowaniem podłoża, a jego zbyt wyrafinowanym modelowaniem stwa-
rzającym poważne trudności w praktycznym stosowaniu. Proponuje algorytmy
obliczeń, stanowiące kompromis pomiędzy dokładnością a czasochłonnością
obliczeń.
Książka adresowana jest do pracowników naukowych wykorzystujących
MRS i MES do tworzenia modelu numerycznego danego zagadnienia;
w szczególności zagadnienia współpracy konstrukcji z podłożem gruntowym.
Skierowana jest również do studentów wydziałów politechnik pragnących
usystematyzować swoją wiedzę z zakresu znanych i stosowanych metod
obliczeniowych. Stanowiłoby ogromną satysfakcję dla autorki, gdyby w wyżej
wspomnianym zakresie książka okazała się pomocna.
Korzystając ze sposobności autorka pragnie wyrazić wdzięczność
Recenzentom: Panu profesorowi Jerzemu Gołasiowi, a w szczególności Panu
profesorowi Maciejowi Gryczmańskiemu za życzliwość i cenne uwagi, pomocne
w uporządkowaniu trudniejszych fragmentów książki. Słowa podziękowania
kieruje również autorka pod adresem Pana Artura Łożyńskiego za znaczącą
pomoc przy komputerowym opracowaniu rysunków i istotne wskazówki
dotyczące składu tekstu.
Rozdział 1
WSTĘP
Zagadnienia fizyki matematycznej czy techniki, takie na przykład jak

poszukiwanie rozkładu pól temperatury, przemieszczenia, naprężenia i od-
kształcenia to zagadnienia brzegowe lub brzegowo-początkowe. Zagadnienia
takie opisane są najczęściej układem równań różniczkowych i warunkami na
brzegu obszaru w przypadku zagadnienia brzegowego, oraz dodatkowo,
warunkami w danej chwili t = t 0 , dla zagadnień brzegowo-początkowych
zależnych od czasu. Każde zagadnienie brzegowe można w sposób ogólny
opisać równaniem operatorowym:
A(u ) + f = 0 (1.1)
określonym w obszarze V , z warunkami brzegowymi
F (u ) = p
G (u ) = q
danymi na części odpowiednio S1 i S 2 brzegu S = S1 ∪ S 2 .

W równaniu (1.1) A oznacza pewien operator różniczkowy. Może to być np.
∂2 ∂2 ∂2
operator Laplace’a: A = 2 + 2 + 2 bądź operator biharmoniczny:
∂x ∂y ∂z
8 Wstęp
∂4 ∂4 ∂4
A= +2 + . Przez u oznaczono funkcję będącą dokładnym
∂x 4 ∂x 2 ∂y 2 ∂y 4
rozwiązaniem równania (1.1), natomiast f jest funkcją daną w obszarze V .
W ogólności zagadnienie brzegowe może mieć charakter liniowy lub nieliniowy
– odpowiednio A będzie operatorem liniowym lub nieliniowym.
Określenie funkcji u jako dokładnego rozwiązania równania (1.1) dotyczy
z reguły wąskiej klasy zagadnień, dla których operator A jest liniowy, obszar u
jest regularny, a warunki brzegowe ciągłe i elementarne. Uzyskane rozwiązania
są ścisłe w sensie matematycznym, często nazywa się je rozwiązaniami
zamkniętymi.
Z chwilą pojawienia się obszarów o złożonych kształtach czy nieciągłych
warunkach brzegowych, metody ścisłe w wyniku stosowania których uzyskuje
się funkcję u jako dokładne rozwiązanie równania operatorowego (1.1)
prowadzą do długich i żmudnych obliczeń, często też w takich przypadkach nie
ma rozwiązań w klasie funkcji elementarnych. Pozostają wówczas do dyspozycji
metody przybliżone rozwiązania równania (1.1), które kosztem zmniejszonej
dokładności obliczeń prowadzą szybciej lub w ogóle umożliwiają uzyskanie
rezultatu.
Generalnie przy przybliżonym rozwiązywaniu równania (1.1) można
wyróżnić dwa podejścia. Pierwsze z nich bazuje na oryginalnym równaniu
operatorowym, drugie na warunku minimum funkcjonału zagadnienia. W myśl
twierdzenia Weinberga analizy funkcjonalnej podejście drugie jest możliwe
wtedy, gdy operator A jest potencjalny. Istnieje wtedy funkcjonał, gradientem
którego jest oryginalny operator różniczkowy. W takim przypadku zamiast
poszukiwać rozwiązania równania (1.1) poszukuje się punktu krytycznego
funkcjonału. Zachodzi bowiem wówczas równoważność pomiędzy równaniem
operatorowym a stowarzyszonym z nim funkcjonałem. Oznacza to ogólnie rzecz
biorąc, że rozwiązaniem równania (1.1) jest funkcja, która minimalizuje
funkcjonał i odwrotnie: wyznaczenie funkcji, która minimalizuje funkcjonał
stanowi rozwiązanie równania (1.1). Podejście bazujące na funkcjonale jest
często łatwiejsze niż rozwiązywanie równania operatorowego. Jest też bardziej
dogodne z punktu widzenia zastosowań praktycznych.
Niezależnie od podejścia, metody przybliżone można podzielić na metody
analityczne, czyli takie, przy użyciu których uzyskuje się rozwiązanie równania
(1.1) w postaci funkcji, oraz numeryczne, kiedy rozwiązanie jest zbiorem
wartości funkcji u w wybranych punktach obszaru V . W obu przypadkach
rozwiązanie jest konstruowane w przestrzeni skończenie wymiarowej. Jeżeli
wymiar tej przestrzeni zostaje sukcesywnie zwiększany, wówczas błąd
aproksymacji ulega zmniejszaniu i rozwiązanie przybliżone zmierza do
rozwiązania dokładnego. W ramach wspomnianego wyżej podziału klasyfikację
poszczególnych metod przybliżonych przedstawiono schematycznie na rys. 1.1.
Wstęp 9
W metodach analitycznych traktuje się ciało jako continuum

trójwymiarowe, czyli rozpatruje się model kontynualny ciała. W metodach
numerycznych ciało traktuje się jako skończony zbiór punktów, między którymi
zakłada się funkcyjny opis zmienności poszukiwanego pola. Rozpatruje się
zatem model dyskretny. Proces dojścia do modelu dyskretnego, zwany
dyskretyzacją, może odbywać się albo na drodze matematycznej – tzw.
dyskretyzacja matematyczna, albo na drodze fizycznej przez podział
rozpatrywanego continuum na skończoną liczbę części – tzw. dyskretyzacja
fizyczna.
Rys. 1.1. Metody przybliżone rozwiązania równania A(u ) + f = 0

10 Wstęp
W rozważaniach zawartych w niniejszej książce przedstawione będą

podstawy teoretyczne oraz wybrane aspekty dwóch metod numerycznych:
metody różnic skończonych (MRS) wraz z jej wariacyjnym ujęciem (WMRS)
i metody elementów skończonych (MES). Metody te są tak skonstruowane, że
aby osiągnąć odpowiednią dokładność rozwiązań, wymagana jest duża liczba
powtarzających się obliczeń przez co niezbędnym staje się zastosowanie
komputera. Stąd też metody te nazwać można metodami komputerowymi.
Nie jest możliwe podanie jednolitego kryterium użyteczności zasto-
sowania każdej z wyżej wymienionych metod komputerowych do danego
zagadnienia. W dużej mierze jest to uwarunkowane specyfiką rozpatrywanego
problemu jak również doświadczeniem i intuicją tego, kto decyduje się na wybór
tej, a nie innej metody numerycznej.
Stwierdzić można, że najbardziej uniwersalną spośród wymienionych
metod jest metoda elementów skończonych. Jest ona wyczerpująco omówiona
w literaturze przedmiotu i podobnie jak metoda różnic skończonych całkowicie
opracowana w swoich podstawach teoretycznych. Metoda elementów brze-
gowych, która w niniejszej książce nie została omówiona, wciąż jest w fazie
rozwojowej i znajduje coraz szersze zastosowanie do rozwiązywania wielu
zagadnień naukowych i inżynierskich.
Warto również podkreślić, iż przy zastosowaniu danej metody
komputerowej wyniki odnoszą się nie do układu rzeczywistego, ale do modelu
obliczeniowego (dyskretnego), który tylko tę rzeczywistość modeluje. Mając
bowiem układ rzeczywisty, którego zachowanie chcemy zbadać, tworzymy
najpierw zbiór założeń dotyczących procesów i zjawisk, które w tym układzie
zachodzą, czyli model fizyczny. W dalszej kolejności formułujemy równania,
które wyrażają zasadnicze własności modelu fizycznego w sposób mate-
matyczny. To dopiero staje się podstawą utworzenia modelu obliczeniowego
rzeczywistego obiektu przy uwzględnieniu własności danej metody nume-
rycznej. Schematycznie przedstawia to rys. 1.2.
Właściwa interpretacja wyników obliczeń programu komputerowego jest
możliwa jedynie przy znajomości zastosowanego w tym programie modelu
obliczeniowego. Ocena adekwatności otrzymanych rezultatów może odbyć się
albo na drodze porównania z rozwiązaniem ścisłym (zazwyczaj niedostępnym),
albo na drodze weryfikacji doświadczalnej.
Nie trzeba chyba podkreślać, że programy, w których użyto różnych
modeli obliczeniowych tego samego obiektu rzeczywistego prowadzić będą do
różniących się od siebie wyników. Nawet jednak przy tych samych modelach
obliczeniowych mogą pojawiać się rozbieżności w wynikach, wynikające ze
stosowania różnych technik numerycznych, na przykład różnych procedur
rozwiązywania układu równań.
Wstęp 11
Rys. 1.2. Schemat tworzenia modeli obiektu rzeczywistego

Rozdział 2
METODA RÓŻNIC SKOŃCZONYCH
2.1. Wprowadzenie
Metoda różnic skończonych jest jedną z najstarszych metod przy-
bliżonych, stosowaną jeszcze przed powstaniem elektronicznych maszyn
cyfrowych jako dyskretna metoda rozwiązywania problemów brzegowych
opisanych równaniami różniczkowymi. Metodę tę rozszerzono następnie na
zagadnienia sformułowane w postaci wariacyjnej i stąd powstała odmiana tej
metody – wariacyjne ujęcie metody różnic skończonych (WMRS).
Metoda różnic skończonych oparta jest na dyskretyzacji matematycznej
równań różniczkowych opisujących continuum. Jej istotą jest zastąpienie
operatorów różniczkowych odpowiednimi operatorami różnicowymi (ilorazami
różnicowymi) określonymi na dyskretnym zbiorze punktów zwanym siatką
dyskretyzacyjną. Elementy tego zbioru nazywane są węzłami. W wyniku takiego
podejścia zagadnienie brzegowe sprowadzone zostaje do układu równań
algebraicznych, w którym niewiadomymi są wartości poszukiwanej funkcji
w węzłach siatki dyskretyzacyjnej.
W swojej wersji podstawowej MRS dotyczy regularnej siatki węzłów, co
stanowi znaczne ograniczenie w zastosowaniu do obszarów o złożonej
geometrii. Do mankamentów metody zalicza się również brak możliwości
lokalnego zagęszczenia siatki i trudność w łączeniu obszarów o różnych
wymiarach (na przykład 1D1 z 2D2). Niewątpliwą zaletą metody jest jej prostota
oraz duża skuteczność dla obszarów regularnych.
1
1D - jednowymiarowy (ang. 1- dimensional)
2
2D - dwuwymiarowy
14 Metoda różnic skończonych
Uogólnienia podstawowej wersji MRS opracowane w ostatnich latach

i umożliwiające uwzględnienie podziału nierównomiernego pokazują duże
możliwości metody i zwiększają jej obszar zastosowań, który staje się
porównywalny z obszarem wiodącej metody komputerowej – metody elemen-
tów skończonych.
Często stosowana jest kombinacja MRS i MES w analizach zagadnień
brzegowo-początkowych. Pierwsza używana jest wtedy do dyskretyzacji
zmiennej czasowej, druga do dyskretyzacji pól przestrzennych.
2.2. Operatory różnicowe

Operatory różnicowe stanowiące punkt wyjścia metody można otrzymać
wychodząc z definicji pochodnej funkcji w punkcie. W tym celu rozpatruje się
dowolną funkcję ciągłą y = f (x) (rys. 2.1).
Rys. 2.1. Funkcja ciągła y = f (x) i jej wartości w określonych punktach
Pochodną funkcji f (x) w punkcie x n , która oznaczona zostanie przez

y ' = f ' ( x n ) nazywa się granicę ilorazu różnicowego
dy f ( x n + ∆x) − f ( x n ) y − yn
y' = n = lim = lim n +1 (2.1)
dx ∆x →0 ∆x ∆x → 0 ∆x
2.2. Operatory różnicowe 15
Wyrażenie powyższe można zapisać w następującej formie:
dy y n+1 − y n ∆y n
n = lim = lim (2.2)
dx ∆x →0 ∆x ∆x →0 ∆x
Przybliżeniem tego wyrażenia jest równanie:
dy ∆y n
≅ (2.3)
dx ∆x
które nazywa się pierwszą różnicą przednią funkcji y w punkcie x n .
Odpowiednio pierwsza różnica wsteczna funkcji y w punkcie x n wyraża się
następująco:
dy ∆y n y − y n −1
n = = n (2.4)
dx ∆x ∆x
Pierwsza różnica centralna zawiera punkty symetrycznie położone względem

punktu x n :
dy ∆y n y n +1 − y n −1
n = = (2.5)
dx ∆x 2∆x
Wiadomo, że pochodna ma swoją interpretację geometryczną. Pochodna y ' ( x n )
jest równa tangensowi kata nachylenia stycznej do osi OX w punkcie y n
(rys 2.2). Jeżeli przyjmie się, że w punkcie x n sieczna b jest równoległa do
stycznej n , wówczas tgβ można zapisać następująco:
y n +1 − y n −1
tgβ = y ' ( x n ) = , (2.6)
2∆x
a to jest właśnie pierwsza różnica centralna.

Rys. 3.2. Styczna i sieczna do krzywej y = f (x)
Oczywiście im mniejsze ∆x , tym bardziej kąt nachylenia stycznej zbliża się do

kąta nachylenia siecznej i tym bardziej zmniejsza się błąd, którym obarczone są
ilorazy różnicowe. W ogólności pochodna w danym punkcie aproksymowana
jest najlepiej przez odpowiednią różnicę centralną, ponieważ ten schemat
różnicowy uwzględnia wartości funkcji w punktach symetrycznie położonych
względem danego punktu. Te symetrycznie położone punkty dają lepszą średnią
wartość pochodnych. Stąd też w dalszych rozważaniach uwzględniona zostanie
jedynie pierwsza różnica centralna.
Pochodne funkcji y = f ( x) wyższych rzędów zastąpione zostaną ilorazami
różnicowymi w następujący sposób:
∆y ∆y
2 2 n+ 1
− n− 1
d y ∆ yn ∆ ∆y n ∆x 2 ∆x 2
2 n → 2
= = =
dx ∆x ∆x ∆x ∆x
y n +1 − y n y n − y n −1 y n +1 − y n y n − y n −1
− − (2.7)
(1 2 ∆x + 2 ∆x) ( 2 ∆x + 2 ∆x)
1 1 1
= = ∆ x ∆x =
∆x ∆x
y n +1 − 2 y n + y n −1
=
∆x 2
∆2 y ∆2 y ∆2 y
−2 +
d4y ∆4 y n 2
∆2 ∆ y n 2 n +1
∆x 2
n
∆x 2
n −1
n → = = ∆x =
dx 4 ∆x 4 ∆x 2 ∆x 2 ∆x 2
y n + 2 − 2 y n +1 + y n − 2( y n +1 − 2 y n + y n −1 ) + y n − 2 y n −1 + y n − 2 (2.8)
= 4
=
∆x
y n + 2 − 4 y n +1 + 6 y n − 4 y n −1 + y n −2
=
∆x 4
Odpowiednie schematy różnicowe można zilustrować graficznie (rys. 2.3).

Rozważa się obecnie funkcję dwóch zmiennych f ( x, y ) . W tym przy-
padku obszar zmienności funkcji f ( x, y ) pokrywa się z siatką o bokach ∆x i
∆y (rys. 2.4) i określa wartości dyskretne w jej węzłach. Rozpatrując dowolny
przekrój y = const. otrzymuje się:
∂f ( x, y ) ∆f f i +1 − f i −1
i → i = (2.9)
∂x ∆x 2∆x
Rys. 2.3. Operatory różnicowe funkcji jednej zmiennej

Rys. 2.4. Obszar zmienności funkcji wraz z siatką dyskretyzacyjną
Z kolei dla przekroju x = const. uzyskuje się:
∂f ( x, y ) ∆f f j − fl
i → i = (2.10)
∂y ∆y 2∆y
Odpowiednio też zachodzi:
∂ 2 f ( x, y ) ∆2 f f i +1 − 2 f i + f i −1
i → i = (2.11)
∂x 2 ∆x 2 ∆x 2
oraz :
∂ 2 f ( x, y ) ∆2 f f j − 2 fi + fl
2 i → 2 i = (2.12)
∂y ∆y ∆y 2
Dla pochodnej mieszanej odpowiedni iloraz różnicowy otrzymuje się w ana-

logiczny sposób:
∆f ∆f f j +1 − f j −1 f l +1 − f l −1
2 2 j− l −
∂ f ( x, y ) ∆ f ∆ ∆f i ∆x ∆x 2∆x 2∆x
i → i = = = =
∂x∂y ∆x∆y ∆y ∆x 2∆y 2∆y
f j +1 − f j −1 − f l +1 + f l −1 (2.13)
= ,
4∆x∆y
podobnie dla pochodnej czwartego rzędu:
∆2 f ∆2 f ∆2 f
i +1 −2 i + i −1
∂ 4 f ( x, y ) ∆4 f ∆2 ∆2 f i ∆y 2 ∆y 2 ∆y 2
i → i = = =
∂x 2 ∂y 2 ∆x 2 ∆y 2 ∆x 2 ∆y 2 ∆x 2
f j +1 − 2 f i +1 + f l +1 f j − 2 fi + fl f j −1 − 2 f i −1 + f l −1
−2 + (2.14)
∆y 2 ∆y 2 ∆y 2
= =
∆x 2
f j +1 − 2 f i +1 + f l +1 − 2 f j + 4 f i − 2 f l + f j −1 − 2 f i −1 + f l −1
=
∆x 2 ∆y 2
Schematy powyższe można zilustrować graficznie (rys. 2.5).

Postępując w ten sposób można napisać schemat różnicowy dowolnej pochodnej
zwykłej lub mieszanej.
Rys. 2.5. Operatory różnicowe funkcji dwóch zmiennych

2.3. Ocena dokładności operatorów różnicowych

Rozpatruje się dowolną funkcję ciągłą y = f (x) (rys. 2.6).
Rys. 2.6. Styczna i sieczna do krzywej y = f (x) w punkcie x n
Przyjęte schematy różnicowe jedynie aproksymują rzeczywistą wartość

pochodnej:
dy
• Rzeczywista wartość pochodnej: tgβ n =
dx
dy ∆y
• Aproksymacja: = tgβ b ≅ tgβ n =
dx ∆x
Rozbieżność pomiędzy rzeczywistą wartością pochodnej a aproksymacją

maleje wraz ze zmniejszaniem się odległości ∆x . Aby ocenić tę rozbieżność,
czyli błąd wynikający z zastąpienia pochodnej ilorazem różnicowym rozwija się
daną funkcję w szereg Taylora.
W tym celu rozpatruje się funkcję jednowymiarową y = f ( x) (rys. 2.7).
Szereg Taylora funkcji f ( x + ξ ) ma postać:
∞ ξ n d n f ( x)
f (x + ξ ) = ∑ , (2.15)
n=0 n! dx n
gdzie pochodna zerowego rzędu jest równa f (x) i 0! = 1 .

2.3. Ocena dokładności operatorów różnicowych 21
Rys. 2.7. Funkcja ciągła y = f (x)
Zgodnie z rys. 2.7 dla f1 i f −1 mamy:
ξ 0 d 0 f ( x0 ) ξ 1 d 1 f ( x0 ) ξ 2 d 2 f ( x0 )
f1 = f ( x 0 + k{
∆x) = + + +
ξ 0! dx 0 1! dx1 2! dx 2
ξ 3 d 3 f ( x0 ) df ξ2 d2 f d3 f (2.16)
+ + K = f ( x0 ) + ξ x0 + x 0 + x0 + K
3! dx 3 dx 2 dx 2 dx 3
df k 2 ∆x 2 d 2 f k 3 ∆x 3 d 3 f
= f 0 + k ∆x x0 + x 0 + x0 + K
dx 2 dx 2 6 dx 3
ξ 0 d 0 f ( x0 ) ξ 1 d 1 f ( x0 ) ξ 2 d 2 f ( x0 )
f −1 = f ( x 0 − ∆
{x) = + + +
ξ 0! dx 0 1! dx1 2! dx 2
ξ 3 d 3 f ( x0 ) df ∆x 2 d 2 f
+ + K = f ( x 0 ) − ∆x x0 + x0 − (2.17)
3! dx 3 dx 2 dx 2
∆x 3 d 3 f
− x0 + K
6 dx 3
df
Schemat różnicowy dla pochodnej x 0 przedstawia się następująco:
dx
df f − f −1
x0 = 1 (2.18)
dx ∆x(1 + k )
Odejmując stronami wyrażenia na f1 i f −1 uzyskuje się:
df ∆x 2 2 d2 f
f1 − f −1 = ∆x (1 + k ) x0 + (k − 1) 2 x 0 +
dx 2 dx
3 3
(2.19)
∆x d f
+ (k 3 + 1) 3 x 0 + K
6 dx
df
Z powyższego wyrażenia można wyznaczyć x 0 . Otrzymuje się:
dx
df 1 ∆x 2 (k 2 − 1) d 2 f
x0 = ( f1 − f −1 ) + x0 +
dx ∆x (1 + k ) 2∆x (1 + k ) dx 2
∆x 3 (k 3 + 1) d 3 f f1 − f −1
+ x 0 + K = +
6∆x(1 + k ) dx 3 ∆x (1 + k )
∆x d2 f ∆x 2  k 3 + 1  d 3 f
+ (k − 1) 2 x 0 +   x0 + K
2 dx 6  1 + k  dx 3
14444444442444444444 3
R (2.20)
df
Pierwszy człon w wyrażeniu (2.20) na pochodną . x 0 to przyjęty
dx
schemat różnicowy (2.18). Pozostałe człony stanowią resztę, którą odrzuca się
zastępując pochodną jej schematem różnicowym. Tak więc przybliżone
df
wyrażenie na pochodną x 0 dane ilorazem różnicowym obarczone jest
dx
błędem:
∆x d2 f ∆x 2  k 3 + 1  d 3 f
R= (k − 1) 2 x 0 +   x0 + K (2.21)
2 dx 6  1 + k  dx 3
2.4. Zastosowanie MRS w analizie statycznej belek 23
W przypadku gdy k ≠ 1 (podział nierównomierny) R → 0 gdy ∆x → 0 , a więc

błąd maleje wraz ze zmniejszaniem siatki podziału. Dla podziału równo-
miernego, czyli gdy k = 1 , R → 0 gdy ∆x 2 → 0 .
Błąd operatora różnicowego nie jest jedynym źródłem powstawania błędu
rozwiązywania zagadnienia metodą różnic skończonych. Dochodzą tu podobnie
jak i w innych metodach numerycznych błędy wynikające z żądania spełnienia
równań różniczkowych w wybranych, a nie w każdym punkcie obszaru oraz
błędy zaokrągleń arytmetycznych.
2.4. Zastosowanie MRS w analizie statycznej belek

2.4.1. Przykład – belka obciążona siłą skupioną
Dana jest belka obustronnie utwierdzona obciążona siłą skupioną jak na
rys. 2.8. Do określenia jest rozkład momentów zginających w belce.
Rozwiązanie dokładne tego zagadnienia przedstawia rys. 2.9.
Zadanie rozwiązane zostanie w sposób przybliżony przy zastosowaniu
metody różnic skończonych. Punktem wyjścia jest równanie różniczkowe osi
odkształconej pręta:
d 4 w( x)
EJ = q ( x) (2.22)
dx 4
w którym: w(x) – funkcja ugięcia belki, q (x) – intensywność obciążenia

równomiernie rozłożonego, EJ – sztywność belki.
Rys. 2.8. Belka obustronnie utwierdzona obciążona w środku siłą skupioną

Rys. 2.9. Rozkład momentów zginających – rozwiązanie dokładne
Rozpatruje się więc zagadnienie brzegowe opisane ogólnie równaniem (1.1),

w którym obecnie:
d4
A≡
dx 4
u ≡ w(x) (2.23)
q( x)
f ≡−
EJ
Warunki brzegowe wynikające ze sposobu podparcia belki są następujące:
w( x = 0) = 0
w( x = l ) = 0
dw
ϕ ( x = 0) = x =0 = 0
dx
dw
ϕ (x = l) = x =l =0
dx (2.24)
Związek pomiędzy funkcją momentów zginających M ( x) a funkcją ugięcia

w(x) ma postać:
d 2 w( x)
M ( x) = − EJ (2.25)
dx 2
W celu uzyskania rozwiązania zadania metodą różnic skończonych dokonuje się

l
podziału belki równomierną siatką dyskretyzacyjną o boku oczka ∆x =
2
(rys. 2.10).
Rys. 2.10. Przyjęta dyskretyzacja belki
W każdym węźle siatki podziału należy napisać odpowiednie równania

różnicowe. Równania te wynikają z zastąpienia występującej w równaniu
opisującym zagadnienie (2.22) pochodnej czwartego rzędu przez odpowiedni
iloraz różnicowy. Otrzymuje się następujące równania różnicowe:
q
węzeł ‘ 1 ‘ : w1″ − 4 w1′ + 6w1 − 4 w2 + w3 = 1 ∆x 4
EJ
q2
węzeł ‘ 2 ‘ : w1′ − 4 w1 + 6 w2 − 4 w3 + w3′ = ∆x 4 (2.26)
EJ
q3
węzeł ‘ 3 ‘ : w1 − 4 w2 + 6w3 − 4 w3′ + w3 ″ = ∆x 4
EJ
W równaniach powyższych wielkości w1′ , w1″ oraz w3′ , w3 ″ są ugięciami

węzłów fikcyjnych. Aby bowiem napisać schemat różnicowy odpowiadający
czwartej pochodnej dla węzła „1” oraz „3” należy „przedłużyć” belkę o dwa
odcinki o długości ∆x , uzyskując po dwa dodatkowe węzły niezbędne
w czwartej różnicy skończonej.
Z kolei biorąc pod uwagę warunki brzegowe
w1 = 0
(2.27)
w3 = 0
wielkością poszukiwaną będzie jedynie ugięcie węzła „2”. Ugięcie to można

obliczyć z równania (2.26) 2 :
q2
w1′ − 4 w1 + 6 w2 − 4 w3 + w3′ = ∆x 4 (2.28)
EJ
′ ′
Występujące w powyższym równaniu ugięcia fikcyjne w1 oraz w3 wyrażone
zostaną poprzez ugięcia węzłów rzeczywistych przy wykorzystaniu warunków
brzegowych:
dw
ϕ 1= 0 ⇒ 1= 0
dx
(2.29)
dw
ϕ 3= 0 ⇒ 3= 0
dx
Zastępując w równaniu (2.29) pierwszą pochodną odpowiednim ilorazem

różnicowym otrzymuje się:
dw w2 − w1′
1= 0: = 0 ⇒ w1′ = w2
dx 2∆x
(2.30)
dw w3′ − w2
3= 0: = 0 ⇒ w3′ = w2
dx 2∆x
Przy uwzględnieniu związków (2.27) i (2.30) w wyrażeniu (2.28) równanie

różnicowe dla węzła „2” przyjmie postać:
q2
8w 2 = ∆x 4 , (2.31)
EJ
W równaniu (2.31) q 2 jest intensywnością obciążenia równomiernie

rozłożonego, przypadającego na węzeł „2”. Ponieważ w węźle „2” działa siła
skupiona, aby określić q 2 , należy zastąpić obciążenie rzeczywiste obciążeniem
równoważnym statycznie, zgodnie z rys. 2.11.
Rys. 2.11. Obciążenie równomiernie rozłożone przypadające na węzeł „2”
P
Tak więc q 2 = . Równanie (2.31) przybiera zatem postać:
∆x
P Pl 3
8w 2 = ∆x 3 = (2.32)
EJ 8 EJ
Stąd otrzymuje się:
Pl 3
w2 = (2.33)
64 EJ
Momenty zginające w wybranych węzłach oblicza się na podstawie zależności

(2.25), w której drugą pochodną zastępuje się odpowiednim ilorazem
różnicowym. Postępując w ten sposób uzyskuje się:
w2 − 2 w1 + w1′
M 1 = − EJ ( )
∆x 2
(2.34)
′
w3 − 2 w3 + w2
M 3 = − EJ ( )
∆x 2
Uwzględniając we wzorze (2.34) związki (2.27), (2.30) oraz uzyskane

rozwiązanie (2.33) ostatecznie otrzymuje się następujące wartości momentów
utwierdzenia M 1 i M 3 :
Pl
M1 = −
8
(2.35)
Pl
M3 = −
8
Jak widać, uzyskane w tym przypadku rozwiązanie przy wykorzystaniu metody

różnic skończonych jest identyczne jak rozwiązanie dokładne.
2.4.2. Przykład – belka obciążona momentem skupionym

Do określenia jest rozkład momentów zginających w belce obustronnie
utwierdzonej, obciążonej w środku momentem skupionym (rys. 2.12).
Rys. 2.12. Belka obustronnie utwierdzona obciążona w środku momentem skupionym
Rozwiązanie dokładne przedstawia rys. 2.13.
Rys. 2.13. Rozkład momentów zginających – rozwiązanie dokładne
Podobnie jak w przykładzie 2.4.1. punktem wyjścia jest równanie różniczkowe

osi odkształconej pręta (2.22). Warunki brzegowe określają związki (2.24).
Przyjętą w metodzie różnic skończonych dyskretyzację belki przedstawiono na

rys. 2.14.
Rys. 2.14. Przyjęta dyskretyzacja belki
Analogicznie jak w przykładzie 2.4.1. kolejnym węzłom siatki podziału

przypisane są następujące równania różnicowe:
q1
węzeł ‘ 1 ‘ : w1 ' '−4 w1 '+6 w1 − 4 w2 + w3 = ∆x 4
EJ
q
węzeł ‘ 2 ‘ : w1 '−4 w1 + 6w2 − 4w3 + w4 = 2 ∆x 4
EJ
q
węzeł ‘ 3 ‘ : w1 − 4 w2 + 6 w3 − 4 w4 + w5 = 3 ∆x 4 (2.36)
EJ
q
węzeł ‘ 4 ‘ : w2 − 4 w3 + 6 w4 − 4 w5 + w5 ' = 4 ∆x 4
EJ
q
węzeł ‘ 5 ‘ : w3 − 4 w4 + 6 w5 − 4 w5 '+ w5 ' ' = 5 ∆x 4
EJ
Wielkości w1 ' , w1 ' ' oraz w5 ' , w5 ' ' występujące w powyższych równaniach są
ugięciami węzłów fikcyjnych. Z uwagi na warunki brzegowe:
w1 = 0
(2.37)
w5 = 0
niezerowe ugięcia wystąpią jedynie w węzłach „2”, „3”, „4”. Będą one
obliczone przy wykorzystaniu równań (2.36), odniesionych do odpowiednich
węzłów. Układ równań na wyznaczenie poszukiwanych wielkości w2 , w3 , w4
ma zatem postać:
q2
31 w1 '−4 w1 + 6w2 − 4w3 + w4 = ∆x 4
EJ
q3
w1 − 4 w2 + 6 w3 − 4 w4 + w5 = ∆x 4 (2.38)
EJ
q
w2 − 4 w3 + 6 w4 − 4 w5 + w5 ' = 4 ∆x 4
EJ
Ugięcia fikcyjne w1 ' oraz w5 ' wyznacza się z warunków brzegowych:
dw w2 − w1′
ϕ 1= 0 ⇒ 1 =0⇒ = 0 ⇒ w1′ = w2
dx 2∆x
(2.39)
dw w′ − w 4
ϕ 5= 0 ⇒ 5 =0⇒ 5 = 0 ⇒ w5′ = w4
dx 2∆x
Uwzględniając zależności (2.39) w równaniach (2.38) ostatecznie otrzymuje się
q2
7 w2 − 4 w3 + w4 = ∆x 4
EJ
q
− 4 w2 + 6 w3 − 4w4 = 3 ∆x 4 (2.40)
EJ
q4
w2 − 4 w3 + 7 w4 = ∆x 4
EJ
Obciążenie przyłożone jest w węźle „2”, natomiast q 3 = q 4 = 0 . Obciążenie to

zastępuje się obciążeniem równoważnym statycznie, zgodnie z rys. 2.15.
W wyniku zastąpienia momentu skupionego parą sił, którą następnie zapisano
w postaci statycznie równoważnego obciążenia równomiernie rozłożonego
uzyskuje się:
P M
q2 = − =−
∆x 2∆x 2
P M (2.41)
q4 = =
∆x 2∆x 2
q3 = 0
Równania (2.40) przybiorą zatem postać:

M
7 w2 − 4 w3 + w4 = − ∆x 2
2 EJ
− 4 w2 + 6 w3 − 4 w4 = 0 (2.42)
M
w2 − 4 w3 + 7 w4 = ∆x 2
2 EJ
Rys. 2.15. Obciążenie statycznie równoważne obciążeniu danemu
Po rozwiązaniu układu (2.42) otrzymuje się:
Ml 2
w2 = −
192 EJ
w3 = 0 (2.43)
2
Ml
w4 =
192 EJ
Momenty utwierdzenia M 1 i M 5 oblicza się na podstawie zależności (2.25)

wyrażonej w różnicach skończonych:
w2 − 2 w1 + w1′
M 1 = − EJ ( )
∆x 2
(2.44)
w5′ − 2 w5 + w4
M 5 = − EJ ( )
∆x 2
Uwzględniając w powyższym związki (2.37) oraz (2.39), otrzymuje się:
2 w2
M
M 1 = − EJ 2
=
∆x 6
2 w4 M (2.45)
M 5 = − EJ 2
=−
∆x 6
Odpowiednie wartości dokładne momentów wynoszą:
M
M1 =
4
(2.46)
M
M5 = −
4
Rozbieżności jakie pojawiają się w tym przypadku pomiędzy rozwiązaniem

dokładnym a rozwiązaniem metodą różnic skończonych wynikają głównie ze
zbyt rzadkiej siatki dyskretyzacyjnej, co szczególnie odzwierciedla się przy
zastąpieniu momentu skupionego parą sił.
2.5. Zastosowanie MRS w analizie statycznej płyt

Obliczenia statyczne płyty izotropowej o stałej grubości i małych ugięciach
sprowadzają się do wyznaczenia funkcji ugięcia spełniającej równanie
różniczkowe:
Eh 3  ∂ 4 w( x, y ) ∂ 4 w( x, y ) ∂ 4 w( x, y ) 
 + 2 +  = q ( x, y ) (2.47)
12(1 − ν 2 )  ∂x 4 ∂x 2 ∂y 2 ∂y 4 
oraz warunki brzegowe zależne od sposobu podparcia płyty.

W równaniu (2.47) przez E oznaczono moduł sprężystości podłużnej, h jest
grubością płyty, w( x, y ) ugięciem powierzchni środkowej, ν współczynnikiem
Poissona, natomiast q ( x, y ) obciążeniem przypadającym na jednostkę
powierzchni płyty. Rozpatruje się zatem zagadnienie brzegowe (1.1), w którym:
2.5. Zastosowanie MRS w analizie statycznej płyt 33
∂4 ∂4 ∂4
A≡ +2 +
∂x 4 ∂x 2 ∂y 2 ∂y 4
u ≡ w( x, y ) (2.48)
q ( x, y )
f ≡−
D
Eh 3
gdzie D = jest sztywnością płyty na zginanie. W celu zastosowania
12(1 − ν 2 )
MRS do rozwiązania płyty rozpatruje się płytę prostokątną i dokonuje podziału
jej obszaru na elementarne podobszary o wymiarach ∆x i ∆y (rys. 2.16).
Rys. 2.16. Płyta z naniesioną siatką dyskretyzacyjną
Równanie powierzchni ugięcia płyty (2.47) można zapisać w następującej

postaci:
q ( x, y )
∇ 2 ∇ 2 w( x, y ) = (2.49)
D
gdzie: ∇ 2 ∇ 2 – operator różniczkowy, zwany bilaplasjanem (podwójny operator

Laplace’a).
Zastępując w powyższym równaniu pochodne cząstkowe odpowia-
dającymi im ilorazami różnicowymi otrzymuje się równanie różnicowe ugięcia
płyty w węźle „k”:
1 2
4
( wk −2 − 4wk −1 + 6 wk − 4 wk +1 + wk + 2 ) + ( wl −1 −
∆x ∆x ∆y 2 2
− 2 wl + wl +1 − 2wk −1 + 4 wk − 2 wk +1 + wi −1 − 2wi + wi +1 ) + (2.50)

1 q
+ 4
( wm − 4 wl + 6 wk − 4wi + wn ) = k
∆y D
Wykorzystano tutaj wcześniej wyprowadzone ilorazy różnicowe dla pochodnej

IV rzędu.
∆y
Oznaczając stosunek przez µ można przedstawić równanie (2.50)
∆x
w postaci:
µ 4 ( wk −2 − 4 wk −1 + 6 wk − 4 wk +1 + wk + 2 ) + 2 µ 2 ( wl −1 − 2 wl +
2 wk −1 + 4 wk − 2wk +1 + wi −1 − 2 wi + wi +1 ) + wm − 4wl + 6 wk i −
(2.51)
q µ 4 ∆x 4
− 4 wi + wm = k
D
W przypadku podziału równomiernego ∆x = ∆y otrzymuje się równanie

różnicowe:
20 wk − 8( wk −1 + wk +1 + wl + wi ) + 2( wl −1 + wl +1 + wi −1 + wi +1 ) +
q k ∆x 4 (2.52)
+ wk − 2 + wk + 2 + wm + wn =
D
Schematyczny zapis tego równania przedstawia rysunek 2.17.

Równania tego typu należy napisać dla każdego węzła siatki podziału.
W węzłach położonych na krawędzi płyty i w pobliżu tej krawędzi oprócz ugięć
punktów „rzeczywistych” pojawią się jeszcze w powyższych równaniach
fikcyjne ugięcia węzłów poza płytą. Ugięcia te wyznacza się na podstawie
warunków brzegowych.
Warunki brzegowe
Równanie (2.47) ugięcia płyty jest cząstkowym równaniem różniczkowym

IV rzędu, rozwiązanie w( x, y ) można więc wyznaczyć tylko dla dwóch
warunków brzegowych na każdej krawędzi. Występują tam jednak trzy siły
przekrojowe: moment zginający, moment skręcający i siła poprzeczna. Siły te
powinny być zrównoważone odpowiadającymi im siłami zewnętrznymi.
W celu otrzymania tylko dwóch warunków brzegowych dodaje się do sił

poprzecznych równoważniki momentów skręcających, otrzymując zastępcze
siły poprzeczne:
 ∂ 3 w( x, y ) ∂ 3 w( x, y ) 
Qx = −D  3
+ ( 2 − ν ) 
 ∂x ∂x∂y 2 
 ∂ 2 w( x, y ) ∂ 3 w( x, y )  (2.53)
Q y = −D 3
+ ( 2 − ν ) 
 ∂y ∂x 2 ∂y 
Rys. 2.17. Operator biharmoniczny
Odpowiednie wyrażenia na momenty zginające oraz moment skręcający

mają postać:
 ∂ 2 w( x, y ) ∂ 2 w( x, y ) 
M x = −D 2
+ ν 
 ∂x ∂y 2 
 ∂ 2 w( x, y ) ∂ 2 w( x, y )  (2.54)
M y = −D 2
+ ν 
 ∂y ∂x 2 
∂ 2 w( x, y )
M x y = −(1 − ν ) D
∂x∂y
Rozpatrzone będzie obecnie równanie różnicowe (2.52) dla węzłów

położonych w pobliżu lub na krawędzi płyty w przypadku różnych sposobów
podparcia płyty.
• Krawędź x = const. całkowicie swobodna (rys. 2.18)
Na krawędzi swobodnej muszą zanikać: moment zginający M x i zastępcza siła

poprzeczna Q x . Zgodnie zatem z zależnościami (2.53) i (2.54) otrzymuje się:
∂ 2 w( x, y ) ∂ 2 w( x, y )
Mx =0⇒ +ν =0
∂x 2 ∂y 2
(2.55)
∂ 3 w( x, y ) ∂ 3 w( x, y )
Qx = 0 ⇒ + (2 − ν ) =0
∂x 3 ∂x∂y 2
Rys. 2.18. Krawędź swobodna płyty z naniesioną siatką dyskretyzacyjną
W równaniu (2.52) dla węzła „k”:
20 wk − 8( wk −1 + wk +1 + wk + wi ) + 2( wl −1 + wl +1 + wi −1 + wi +1 ) +
q k ∆x 4
+ wk − 2 + wk + 2 + wm + wn =
D
ugięcia wk − 2 , wl −1 , wk −1 oraz wi −1 są obecnie ugięciami fikcyjnymi węzłów

położonych poza płytą. Ugięcia te wyraża się poprzez ugięcia węzłów
położonych w obszarze płyty wykorzystując warunki brzegowe (2.55),
odniesione do odpowiednich węzłów. Warunki te mają postać:
(M x ) l = 0
(M x ) k =0
(2.56)
(M x ) i = 0
(Q x ) k =0
Wyrażenia na momenty M x oraz zastępczą siłę poprzeczną Q x zastępuje się
przez odpowiednie schematy różnicowe. Przykładowo dla węzła „k” uzyskuje
się:
 ∂2w ∂2w   ∆2 wk ∆2 wk 
Mx = − D  2 + ν  k = −D  + ν =
∂y 2   ∆x 2 
k
 ∂x  ∆y 2 
= − D [ wk +1 − 2wk + wk −1 + ν ( wl − 2wk + wi )]
∂3w ∂3w   ∆3 wk
Qx k = − D  3 + (2 − ν )  k = −D  + (2.57)
 ∂x ∂x∂y 2   ∆x∆y
2
∆3 wk 
+ (2 − ν )  = − D [2wk −1 − wk − 2 − 2 wk +1 + wk + 2 +
∆x∆y 2 
+ (2 − ν )(wl +1 − 2 wk +1 + wi +1 − wl −1 + 2 wk −1 − wi −1 )]
Postępując w ten sposób, z warunków (2.56) otrzymuje się następujące

wyrażenia dla ugięć węzłów fikcyjnych:
wl −1 = 2(1 + ν ) wl − wl +1 − νwm − νwk
wk −1 = 2(1 + ν ) wk − wk +1 − νwl − νwi

wi −1 = 2(1 + ν ) wi − wi +1 − νwk − νwn
(2.58)
wk − 2 = wk + 2 − 4(3 − ν ) wk +1 + 6(2 + 2ν − ν 2 ) wk −
− 4(1 + 2ν − ν 2 )(wi + wl ) + ν (2 − ν )(wm + wn ) +

+ 2(2 − ν )( wl +1 + wi +1 )
Uwzględniając powyższe związki w równaniu różnicowym dla węzła „k”

uzyskuje się ostatecznie równanie różnicowe dla węzła położonego na krawędzi
swobodnej płyty:
2(8 − 4ν − 3ν 2 ) wk − 4(3 − ν ) wk +1 + 2 wk + 2 − 4(2 + ν )(1 − ν )(wl + wi ) +

2 q k ∆x 4 (2.59)
+ 2(2 − ν )( wl +1 + wi +1 ) + (1 − ν )( wm + wn ) =
D
Schematyczny zapis tego równania przedstawia rysunek 2.19.
Rys. 2.19. Współczynniki w równaniu różnicowym węzła położonego

na krawędzi swobodnej
• Krawędź x = const. swobodnie podparta (rys. 2.20)
Rys. 2.20. Krawędź swobodnie podparta płyty z naniesioną siatką dyskretyzacyjną
Warunki brzegowe:
∂2w ∂2w ∆2 wk −1
w = 0, M x k −1 =0⇒( +ν ) k −1 =0⇒ =0 (2.60)
∂x 2 ∂y 2 ∆x 2
Z wyrażenia (2.60) otrzymuje się:
wk − 2wk −1 − wk − 2
=0 (2.61)
∆x 2
Ponieważ zgodnie z zależnością.(2.60) 1 wk −1 = 0 , zachodzi:
wk −2 = − wk (2.62)
Postępując analogicznie jak poprzednio, uzyskuje się równanie różnicowe dla

węzła „k” leżącego w pobliżu krawędzi swobodnie podpartej. Równanie to
przyjmie postać zapisaną schematycznie na rysunku 2.21.

w pobliżu krawędzi swobodnie podpartej płyty
• Krawędź x = const. utwierdzona zupełnie (rys. 2.22)
W tym przypadku warunki brzegowe są następujące:
∂w ∆wk −1
w = 0, ϕ x = k −1 =0⇒ =0 (2.63)
∂x ∆x
Z warunku (2.63) otrzymuje się:
wk − wk − 2
= 0 ⇒ wk − 2 = wk (2.64)
2∆x
Rys. 2.22. Krawędź zamocowana płyty z naniesioną siatką dyskretyzacyjną

Schematyczny zapis równania różnicowego ugięcia płyty w węźle „k”, po

uwzględnieniu zależności (2.64), przedstawia rysunek 2.23.

w pobliżu krawędzi zamocowanej
2.5.1. Przykład – płyta obciążona równomiernie na całym

obszarze
Dana jest płyta o schemacie statycznym i obciążeniu jak na rys. 2.24. Do
wyznaczenia są ugięcia płyty.
Rys. 2.24. Płyta obciążona równomiernie na całym obszarze

W celu rozwiązania zadania metodą różnic skończonych dzieli się obszar płyty
siatką dyskretyzacyjną o boku oczka ∆x = a 4 (rys. 2.25).
Rys. 2.25. Przyjęta dyskretyzacja płyty
Z uwagi na symetrię geometrii i obciążenia obliczenia ograniczone będą do

ćwiartki płyty, zgodnie z rys. 2.26.
Warunki brzegowe
• Krawędź swobodnie podparta x = 0 , x = a
w=0
∂2w ∂2w
Mx =0⇒ +ν =0 (2.65)
∂x 2 ∂y 2
∂2w
ale wzdłuż krawędzi w = 0 więc =0
∂x 2
Rys. 2.26. Przyjęty do obliczeń obszar płyty
Z zależności (2.65) wynikają następujące równości:
w9 = w1′ = w5′ = w4 = w4& = w8 = 0

2
∂ w w1 − 2 w1′ + w1′′
2 1′ =0⇒ = 0 ⇒ w1′′ = − w1
∂x ∆x 2
∂2w w − 2w5′ + w5′′
5′ =0⇒ 5 = 0 ⇒ w5′′ = − w5
∂x 2 ∆x 2 (2.66)
∂2w w3 − 2 w4 + w3
2 4 =0⇒ = 0 ⇒ w3 = 0
∂x ∆x 2
∂2w w7 − 2 w8 + w7
8 =0⇒ = 0 ⇒ w7 = 0
∂x 2 ∆x 2
• Krawędź zamocowana y = −a
w=0
∂w (2.67)
ϕ y= =0
∂y
Zatem:
w1& = w2& = w3& = 0

∂w w&1& − w1
1&
=0⇒ = 0 ⇒ w&1& = w1
∂y 2∆x
∂w w&2& − w2 (2.68)
2&
=0⇒ = 0 ⇒ w&2& = w2
∂y 2∆x
∂w w&3& − w3
3&
=0⇒ = 0 ⇒ w&3& = w3
∂y 2∆x
Równania różnicowe
• Węzeł „1”
q1 ∆x 4
19 w1 − 8w2 + w3 − 8w1 + w&1& − 8w5 + w1 + 2w2& + 2 w6 =
D
• Węzeł „2”
20 w2 − 8w( w1 + w2& + w3 + w6 ) + 2( w1& + w3& + w5 + w7 ) + w1′ + w4 +

q 5 ∆x 4
+ w&2& + w2 =
D
• Węzeł „5”
q 5 ∆x 4
19 w5 − 8w6 + w7 − 8w1 + w1& + 2 w2 + 2w2 =
D
• Węzeł „6”
20 w6 − 8( w5 + w7 + w2 + w2 ) + 2( w1 + w3 + w1 + w3 ) + w5′ + w8 +
q 6 ∆x 4
+ w2& + w2& =
D
Po uwzględnieniu w powyższych równaniach warunków brzegowych (2.66)

i (2.68), równania różnicowe w poszczególnych węzłach przyjmą ostateczną
postać:
q1 ∆x 4
Węzeł '1' 21w1 − 8w2 − 8w5 + 2 w6 =
D
q ∆x 4
Węzeł '2' − 8w1 + 22w2 + 2 w5 − 8w6 = 2
D
(2.69)
Węzeł q ∆x 4
'5' − 16 w1 + 4 w2 + 19 w5 − 8w6 = 5
D
Węzeł q 6 ∆x 4
'6' 4 w1 − 16 w2 − 8w5 + 20 w6 =
D
gdzie: q1 = q 2 = q 5 = q 6 = q
Po rozwiązaniu układu (2.69) otrzymuje się:
q∆x 4
w1 = 0,53
D
q∆x 4
w2 = 0,8
D
(2.70)
q∆x 4
w5 = 0,67
D
q∆x 4
w6 = 0,86
D
2.6. Wariacyjne ujęcie metody różnic skończonych

(WMRS)
W odróżnieniu od klasycznej metody różnic skończonych, ujęcie waria-
cyjne bazuje nie na danym równaniu różniczkowym, ale na odpowiednim
funkcjonale. Istota metody polega teraz na zastąpieniu występujących w
funkcjonale pochodnych odpowiednimi ilorazami różnicowymi. W ten sposób
funkcjonał zostaje sprowadzony do funkcji wielu zmiennych, której argumen-
tami są dyskretne wartości poszukiwanej funkcji w wybranych punktach siatki
podziału.
W przypadku zagadnień jednowymiarowych (1-D) dany odcinek dzieli się

na krótsze odcinki, natomiast w zagadnieniach dwuwymiarowych (2-D) dany
obszar na podobszary. Całkowanie funkcjonału odnoszące się do całego obszaru
zamienia się na odpowiednie sumy. Układ równań algebraicznych, wyzna-
czający wartości funkcji w poszczególnych węzłach siatki otrzymuje się
z warunku ekstremum funkcji wielu zmiennych:
∂F (u i )
= 0 (i = 1,2, K , n) (2.71)
∂u i
gdzie F (u i ) – funkcja wielu zmiennych, otrzymana przez zastąpienie

w funkcjonale F (u ) pochodnych odpowiednimi ilorazami różnicowymi, u i –
poszukiwane wielkości dyskretne, n – liczba niewiadomych układu.
2.6.1. Przykład – belka obciążona siłą skupioną

Dana jest belka obustronnie utwierdzona jak w przykładzie 2.4.1. (rys.
2.8). Do wyznaczenia jest ugięcie belki.
Zadanie rozwiązane będzie przy wykorzystaniu wariacyjnego ujęcia metody
różnic skończonych. Punktem wyjścia jest funkcjonał energii sprężystej belki:
2
1 l  d 2 w( x) 
Eb = ∫ EJ   dx (2.72)
2 0  dx 2 
w którym EJ oznacza sztywność belki, w( x) - funkcję ugięcia belki.

Przyjęto analogiczną jak w Przykładzie 2.2.1. dyskretyzację belki (rys. 2.10).
Przyjęto ponadto, że EJ = const. na całej długości belki.
Funkcjonał (2.72) przyjmie wówczas postać:
2
1 l  d 2 w( x ) 
Eb = EJ ∫   dx (2.73)
2 0  dx 2 
Występująca w wyrażeniu (2.73) druga pochodna zastąpiona zostaje dla każdego

punktu podziału odpowiednim ilorazem różnicowym. Ilorazy te dla kolejnych
węzłów przedstawiono na rys. 2.27.
W wariacyjnym ujęciu metody różnic skończonych należy spełnić
w przyjętym układzie niewiadomych tylko geometryczne warunki brzegowe,
2.6. Wariacyjne ujęcie metody różnic skończonych (WMRS) 47
warunki statyczne spełnione są w sposób naturalny przez sam funkcjonał. Dla

rozważanego sposobu podparcia warunki brzegowe są następujące:
w1 = w5 = 0
dw w2 − w1′
ϕx 1 = 0 ⇒ 1 =0⇒ = 0 ⇒ w1′ = w2
dx 2∆x (2.74)
dw w − w2
ϕx 3 =0⇒ 3 = 0 ⇒ 3′ = 0 ⇒ w3′ = w2
dx 2∆x
Rys. 2.27. Ilorazy różnicowe dla węzłów siatki podziału
Po uwzględnieniu związków (2.74) w kolejnych węzłach (rys. 2.27), schematy

różnicowe przyjmą w rezultacie postać przedstawioną na rys. 2.28.
Rys. 2.28. Ilorazy różnicowe po uwzględnieniu warunków brzegowych

Całkowita energia potencjalna układu jest sumą energii sprężystej belki

i pracy sił zewnętrznych:
Π C = Eb + L (2.75)
gdzie L jest potencjałem sił zewnętrznych. Potencjał ten wyraża się

następująco:
l
L z = − ∫ qw( x) dx (2.76)
0
gdzie q jest intensywnością obciążenia równomiernie rozłożonego.

W analizowanym przykładzie siłę skupioną zastępuje się równoważnym
obciążeniem ciągłym jak w Przykładzie 2.4.1. Zatem q 2 = P ∆x , a praca
obciążenia zewnętrznego wyniesie:
P
L2 = − ∫ q 2 w2 dx = − ∫ w2 dx (2.77)
∆x
Zapisując wyrażenie (2.75) w sposób dyskretny otrzymuje się:
1  2 w  2  − 2 w2 
2
 2 w2 
2 
Π C = EJ ∑  2  0,5 + 
2
2
 +  2
 0,5  ∆x −
2 l  ∆x   ∆x   ∆x  
 (2.78)
P
− ∑ w2 ∆x
∆x q
W wyrażeniu powyższym przez ∑ oraz ∑ oznaczono umownie sumowanie

l q
odpowiednio po długości belki i długości obciążenia ciągłego.
Równania algebraiczne wyznaczające ugięcia kolejnych węzłów otrzymuje się
z warunku (2.71) , który obecnie przybierze postać:
∂Π C
= 0 (i = 1,2, K, n) (2.79)
∂wi
gdzie n jest liczbą węzłów o niezerowych ugięciach.

2.6. Wariacyjne ujęcie metody różnic skończonych (WMRS) 49
W rozważanym przykładzie n = 1 i zależność (2.79) sprowadzi się do jednego

równania na wyznaczenie ugięcia w2 . Równanie to ma postać:
∂Π C 1 1
= EJ ⋅ 3 ⋅ 16 w2 = P (2.80)
∂w2 2 ∆x
Stąd:
P∆x 3 Pl 3
w2 = = (2.81)
8 EJ 64 EJ
Warto zwrócić uwagę, iż oczywiście uzyskuje się ten sam jak przy
zastosowaniu klasycznej metody różnic skończonych wynik (przykład 2.4.1),
jednak przy znacznie mniejszym nakładzie pracy rachunkowej.
Rozdział 3
METODA ELEMENTÓW SKOŃCZONYCH
3.1. Wprowadzenie
Metoda elementów skończonych (MES) powstała w latach pięćdziesiątych
ubiegłego stulecia, choć jej początki można odnieść do pracy Couranta z 1942
roku. Zakres zastosowań metody koncentrował się początkowo głównie wokół
zagadnień mechaniki konstrukcji oraz problemów inżynierii lądowej i przemysłu
lotniczego. Ważnym przyczynkiem teoretycznym do rozwoju metody było
wykazanie, że jest ona w rzeczywistości pewną odmianą znanej procedury
Rayleigha-Ritza, bazującej na minimalizacji funkcjonału. Wykorzystanie
procedury minimalizacji umożliwiło zastosowanie metody w innych niż
mechanika konstrukcji dziedzinach. Zaczęto ją stosować do problemów
opisanych równaniami Laplace’a i Poissona, tj. w mechanice cieczy,
w zagadnieniach przepływu ciepła i innych. MES nabrała znamion ogólnej
metody numerycznej rozwiązywania równań czy układu równań różniczkowych
zwyczajnych lub cząstkowych.
Szczególnie dynamiczny rozwój MES datuje się od połowy lat 80.
ubiegłego wieku, głównie w związku z szybkim rozwojem elektronicznej
techniki obliczeniowej. W chwili obecnej MES jest powszechnie stosowaną
metodą numeryczną, a w zagadnieniach mechaniki konstrukcji – wiodącą, na
której bazuje zdecydowana większość profesjonalnych programów kompu-
terowych (na przykład ANSYS, COSMOS/M, ABAQUS, ADINA). O tak dużej
popularności metody zadecydowała zarówno prostota jej postaci matematycznej,
jak również szereg korzystnych cech wyrażających się m.in.: możliwością
obliczania obszarów o złożonych, a także krzywoliniowych brzegach,
52 Metoda elementów skończonych
stosowania różnych gęstości podziału na elementy czy też automatycznego

uwzględniania warunków brzegowych.
Z matematycznego punktu widzenia MES stanowi uogólnienie metody
Galerkina (bazującej na równaniu różniczkowym) i Ritza (bazującej na
funkcjonale). W obu tych metodach przyjmuje się przybliżone rozwiązanie
równania operatorowego (1.1) w postaci związku:
N
u ( x, y, z ) ≅ uˆ ( x, y, z ) = ∑ a iϕ i ( x, y , z ) (3.1)
i =1
gdzie N = 1,2, K , ∞ .
W metodzie elementów skończonych związek powyższy określony jest w
obszarze elementu, a funkcje bazowe ϕ i dobiera się tak, aby zapewniona została
ciągłość poszukiwanej funkcji wzdłuż brzegów elementów. Najczęściej funkcje
te przyjmuje się w postaci wielomianów. Zapewnia to wraz z zastosowaniem
jako elementów prostych brył (czworościany, sześciościany) lub prostych figur
geometrycznych (trójkąty, czworokąty) prostotę i efektywność obliczeń metody.
3.2. Idea metody elementów skończonych

Zasadnicza koncepcja MES, dotycząca standardowego sformułowania
metody, bazuje na założeniu, że dowolna funkcja ciągła może być
aproksymowana przez model dyskretny składający się ze zbioru kawałkami
ciągłych funkcji określonych na skończonej liczbie podobszarów. Funkcje
kawałkami ciągłe określa się wykorzystując wartości funkcji ciągłej w skoń-
czonej liczbie punktów wewnątrz obszaru jej określoności.
Proces tworzenia modelu dyskretnego można przedstawić następująco:
• obszar jest podzielony na podobszary zwane elementami skończonymi,
które połączone są ze sobą w węzłach i wspólnie odzwierciedlają kształt
rozpatrywanego obszaru
• wartości poszukiwanej funkcji w węzłach uznaje się jako niewiadome, które
mają być wyznaczone
• równanie różniczkowe opisujące problem brzegowy odniesione zostaje do
elementu, a jego rozwiązanie aproksymowane w obszarze elementu za
pomocą ciągłej funkcji przy użyciu wartości węzłowych.
W dalszej kolejności formułuje się równania dla elementów. Z chwilą
kiedy są one określone, elementy zostają „scalone” ze sobą, co od strony
fizycznej prowadzi do odtworzenia kształtu rozpatrywanego obszaru, a od strony
matematycznej do uformowania globalnego układu równań umożliwiającego
wyznaczenie wartości poszukiwanej funkcji w wybranych punktach (węzłach)
obszaru. Rozwiązanie problemu staje się w efekcie kawałkami ciągłą
3.2. Idea metody elementów skończonych 53
aproksymacją uzyskaną przez połączenie rozwiązań na wszystkich elementach.

Typową aproksymację MES dla problemu jednowymiarowego przedstawiono
na rys. 3.1.
Rys. 3.1. Typowa aproksymacja MES rozwiązania zagadnienia jednowymiarowego
Ogólnie tok postępowania w MES można ująć w następujący sposób:
1. Podział obszaru na elementy połączone ze sobą w punktach węzłowych.

2. Specyfikacja elementu. Przybliżone rozwiązanie zagadnienia dla każdego
elementu:
• wybór metody przybliżonej (Ritza, Galerkina)
• dobór funkcji interpolujących
• sformułowanie dla wybranej metody i interpolacji układu równań
algebraicznych elementu
3. Agregacja. Budowa na podstawie równań elementów jednego dla całego
obszaru układu równań algebraicznych z uwzględnieniem warunków
ciągłości między elementami.
4. Uwzględnienie w globalnym układzie równań warunków brzegowych.
5. Rozwiązanie otrzymanego układu równań. Uzyskanie dyskretnych wartości
poszukiwanej wielkości.
6. Obliczenie na podstawie wyznaczonych wartości dyskretnych innych
wielkości od nich zależnych.
Przedstawiony sposób postępowania można w sposób schematyczny zobra-
zować jak na rys.3.2.
Rys. 3.2. Tok postępowania w MES

3.3. Dyskretyzacja 55
3.3. Dyskretyzacja
Dyskretyzacja obejmuje podział obszaru na podobszary zwane elementami
skończonymi i przypisanie im zbioru węzłów. Element skończony jest
całkowicie zdefiniowany, jeżeli jednoznacznie określona jest jego geometria
(wielkość, kształt), zbiór węzłów oraz parametry węzłowe (niewiadome
w węźle). Podziału na elementy skończone dokonuje się w sposób schema-
tycznie przedstawiony na rys. 3.3.
Rys. 3.3. Sposób dokonywania podziału na elementy
Podział obszaru na elementy uwarunkowany jest geometrią obszaru oraz

przewidywanym przebiegiem poszukiwanej funkcji. Tradycyjnie w miejscach,
w których istnieje możliwość wystąpienia nagłych zmian poszukiwanej funkcji
(otwory, brzegi obszaru, regiony lokalnie obciążone, styki materiałów o różnych
własnościach) należy zagęścić siatkę podziału na elementy. Takie podejście,
najczęściej stosowane, znane jest jako tzw. wersja h metody. Bardziej

nowoczesne podejścia do tego problemu polegają na podwyższeniu rzędu
wielomianu aproksymacyjnego elementów (metoda p ) lub na zmianie
usytuowania węzłów w elemencie (metoda r ). Możliwe są także kombinacje
wyżej wymienionych metod, czyli stosowanie w rejonach nagłych zmian
poszukiwanej funkcji MES w ujęciu metody hp , hr lub hpr .
Uzyskanemu w wyniku podziału obszaru na podobszary elementowi
przypisuje się zbiór punktów węzłowych, które w MES nie mają charakteru
punktu geometrycznego. Węzeł może przemieszczać się w przestrzeni, a także
doznawać obrotu. W węźle skoncentrowana jest niejako pełna informacja
o zachowaniu się elementu i jego własnościach. Teoretycznie metoda nie
narzuca żadnych ograniczeń odnośnie umiejscowienia węzłów w elemencie.
Praktycznie jednak z uwagi na racjonalizację obliczeń stosuje się w kolejności
węzły:
• wierzchołkowe (stosowane zawsze)
• krawędziowe
• ścienne
• wewnętrzne
Przy „rozseparowaniu” elementów węzły wierzchołkowe, krawędziowe
i ścienne zostają rozdzielone na węzły przynależne do poszczególnych
elementów. Montaż elementów wymaga, aby zachodziła zgodność wielkości
węzłowych wszystkich elementów w węźle wspólnym, z uwagi na konieczność
zapewnienia ciągłości funkcji (ewentualnie jej pochodnych) między elementami.
Obrazowo przedstawia to rys. 3.4.
Montaż elementów wymaga więc (rys. 3.4) spełnienia warunków:
u33 = u41 = u12 = u21 = u5

(3.2)
v33 = v44 = v21 = v12 = v5
Kolejnym krokiem do zdefiniowania elementu jest określenie niewiadomych

w węźle, czyli przypisanie węzłowi wielkości uogólnionych poszukiwanej
funkcji. Liczba niewiadomych w węźle określana jest mianem stopni swobody
węzła (SSW) i związana jest z przyjętym modelem fizycznym i matematycznym
rozpatrywanego zagadnienia.
Istotnym aspektem dyskretyzacji jest również numeracja węzłów
i elementów. Teoretycznie może ona przebiegać w sposób dowolny. Dla danego
obszaru i podziału wszystkie możliwe numeracje węzłów prowadzą do
globalnego układu równań o tych samych wymiarach i o tej samej liczbie
członów niezerowych. Jednak numeracja węzłów i elementów wpływa na
szerokość pasma globalnej macierzy układu i na rozkład elementów zerowych
3.3. Dyskretyzacja 57
w paśmie, a więc zarówno na rezerwację pamięci, jak i na czas rozwiązywania

układu równań.
Rys. 3.4. Rozseparowanie i montaż elementów skończonych

przy uwzględnieniu warunków „zszycia” elementów
Ogólnie można wskazać, że numeracja węzłów przebiegająca wzdłuż krótszego

wymiaru danego obszaru prowadzi do zmniejszenia szerokości pasma macierzy
globalnej. Wpływ numeracji węzłów na rozkład elementów globalnej macierzy

układu i szerokość półpasma przedstawiono przykładowo dla obszaru prosto-
kątnego na rys. 3.5.
Wspomniany problem posiada jednak szerszy aspekt i najczęściej przy
rozwiązywaniu dużych zadań buduje się specjalne programy generujące siatkę
podziału na elementy (tzw. pre-processing).
Rys. 3.5. Wpływ numeracji węzłów na strukturę globalnej macierzy sztywności układu:
a) numeracja dowolna, b) numeracja wzdłuż dłuższego boku, c) numeracja wzdłuż
krótszego boku (Nikishkov G. P., http://www.u-aizu.ac.jp/~niki)
3.4. Klasyfikacja elementów skończonych 59
3.4. Klasyfikacja elementów skończonych

Z liczbą węzłów w elemencie wiąże się stopień wielomianu użytego do
aproksymacji poszukiwanej wielkości w obszarze elementu. Im większa liczba
węzłów, tym wyższy stopień wielomianu, a tym samym dokładniejsza
aproksymacja. Zarazem jednak większa globalna liczba węzłów układu
w sposób bezpośredni rzutuje na wymiar otrzymanego układu równań.
Stopień użytego wielomianu determinuje podział elementów skończonych
na:
• simpleksy – wielomian składa się z członu stałego oraz członów liniowych;
liczba współczynników wielomianu wynosi nw + 1 , gdzie nw oznacza liczbę
wymiarów przestrzeni
• kompleksy – wielomian posiada stałe i liniowe człony oraz człony
wyższych rzędów; elementy tego typu mogą mieć ten sam kształt jak
simpleksy, posiadają jednak dodatkowe węzły brzegowe lub wewnętrzne;
liczba węzłów w kompleksie LW > nw + 1
• multipleksy – wielomian zawiera człony stałe, liniowe i wyższego rzędu; w
celu zapewnienia ciągłości funkcji między elementami brzegi elementu
muszą być równoległe do osi układu odniesienia
Zastosowanie elementów wyższego rzędu (kompleksy, multipleksy)
prowadzi do lepszego odzwierciedlenia rzeczywistego rozkładu poszukiwanej
funkcji przy mniejszej w stosunku do elementów prostych (simpleksy) liczbie
elementów, na który zostaje podzielony dany obszar.
W elementach czworobocznych (sześciościennych) występują dwa
niezależne zbiory węzłów – jeden określa transformację współrzędnych (kształt
elementu), drugi – wielomian interpolacyjny. Oznaczając liczbę węzłów w tych
zbiorach odpowiednio LWK oraz LWW elementy skończone można sklasy-
fikować następująco:
• elementy subparametryczne – LWK > LWW
• elementy izoparametryczne - LWK = LWW
• elementy superparametryczne – LWK < LWW
Najczęściej stosuje się odwzorowanie izoparametryczne, które ma na celu lepsze
zbliżenie się do rozwiązań ścisłych, umożliwia ponadto uwzględnienie
zakrzywionego kształtu geometrycznego rozpatrywanego obszaru. Elementy
izoparametryczne można stosować do zagadnień opisanych funkcjonałem,
w którym występują co najwyżej pierwsze pochodne poszukiwanej wielkości.
Elementy subparametryczne mogą być stosowane do szerszej niż
izoparametryczne klasy zagadnień (w funkcjonale mogą występować dowolne
pochodne), geometria jednak elementu jest ograniczona do możliwie prostych
kształtów. Odwzorowanie superparametryczne ma charakter teoretyczny, tego
rodzaju elementy znajdują też sporadycznie zastosowanie na przykład

w modelowaniu powłok średniej grubości.
Ze względu na liczbę stopni swobody w węźle elementy można podzielić
na dwie klasy:
• Lagrange’owska – element posiada nw stopni swobody w węźle.
Stosowana, jeśli wielkościami węzłowymi są wartości poszukiwanej funkcji
(w zagadnieniach mechaniki konstrukcji na przykład element kratowy)
• Hermite’owska – element posiada więcej niż nw stopni swobody w węźle;
stosowana, jeśli wielkościami węzłowymi są nie tylko wartości
poszukiwanej funkcji, ale także jej pochodne (w mechanice konstrukcji na
przykład element belkowy).
Z punktu widzenia zastosowań w mechanice konstrukcji elementy można
również sklasyfikować następująco:
• elementy podstawowe (kratowy, prętowy, belkowy, ramowy),
• elementy ciągłe (płytowy, tarczowy, powłokowy),
• elementy specjalne (na przykład elementy umożliwiające uwzględnienie
pęknięć, elementy nieskończone, elementy umożliwiające uwzględnienie
struktury „plastra miodu”),
• makroelementy (złożone z elementów prostych).
Przykładowe elementy z wyżej wymienionych grup zostaną omówione
w rozdziale 5.
Elementy można również uszeregować ze względu na ich sposób
formułowania i specyficzne własności:
• zablokowane izoparametrycznie – wprowadzone w celu wyeliminowania
niekorzystnych zjawisk pojawiających się w elementach izoparame-
trycznych (szczególnie niższego rzędu) – na przykład „lockingu” (ekstre-
malnej sztywności); parametrami węzłowymi są jedynie przemieszczenia
• mieszane – konstruowane na bazie funkcjonału Hellingera-Reissnera.
Parametrami węzłowymi są przemieszczenia i naprężenia (wprowadzone
głównie dla modelowania materiałów nieściśliwych oraz polepszenia
zachowania się modeli przemieszczeniowych niższego rzędu)
• hybrydowe – z identyfikacją pola naprężeń wewnątrz i przemieszczenia na
brzegach, formułowane w oparciu o energię komplementarną (szczególnie
użyteczne w przypadkach, w których występują trudności w osiągnięciu
ciągłości założonego pola przemieszczeń między elementami – na przykład
cienkie płyty, powłoki
• o założonych odkształceniach – rozkład pola odkształceń wewnątrz
elementu założony według pewnej funkcji, w większości konstruowane na
bazie twierdzenia Lagrange’a o minimum energii potencjalnej układu,
szczególnie przydatne przy analizie płyt i powłok.
3.5. Równania metody elementów skończonych 61
Programy profesjonalne w swoich pakietach oferują szereg typów elementów

skończonych. Przykładowo program COSMOS/M zawiera około 200 rodzajów
elementów skończonych. Warto podkreślić, że wybór elementu do analizy
danego zagadnienia wiąże się bezpośrednio z akceptacją modelu matema-
tycznego, na bazie którego element został skonstruowany i ma bezpośrednie
implikacje w odniesieniu do uzyskanych rozwiązań.
3.5. Równania metody elementów skończonych

Punktem wyjścia do uzyskania równań MES jest transformowanie opisu
problemu brzegowego w postaci równań różniczkowych (sformułowanie mocne)
do sformułowania całkowego (globalnego). Przekształcenia tego można dokonać
przy wykorzystaniu odpowiedniej zasady wariacyjnej lub metody ważonych
residuów. W pierwszym przypadku bazuje się na słabym sformułowaniu
analitycznym lub definiuje się problem minimalizacji pewnego funkcjonału.
Metoda ważonych residuów przekształca natomiast lokalne sformułowanie
zagadnienia brzegowego w całkową postać słabą.
Odpowiednie sformułowania problemu brzegowego, a więc sformułowanie
mocne wyrażone poprzez równania różniczkowe i warunki brzegowe,
sformułowanie wariacyjne wyrażone przez funkcjonał oraz sformułowanie
słabe, określone jako uśrednione (ważone) równanie całkowe dla obszaru
determinują zastosowanie odpowiednich metod przybliżonych rozwiązywania
równań różniczkowych. Do przybliżonego rozwiązania zagadnienia
sformułowanego jako problem minimalizacji funkcjonału wykorzystuje się
metodę Rayleigha – Ritza, a w przypadku problemu sformułowanego w postaci
formy słabej stosuje się najczęściej metodę Galerkina. W obu metodach przy
podejściu analitycznym postuluje się rozwiązanie przybliżone w postaci związku
(3.1), w którym ϕ i są funkcjami bazowymi, których dobór związany jest z daną
metodą. Funkcje te rozpięte na całym analizowanym obszarze muszą spełniać
warunki brzegowe dane na jego brzegu. Przy podejściu numerycznym – MES
zarówno przy wykorzystaniu sformułowania Rayleigha – Ritza jak i Galerkina
zależność (3.1) odnosi się do poziomu elementu:
u e ( x, y, z ) ≈ uˆ e ( x, y, z ) = ∑ a I ϕ I , (3.3)
I
gdzie I oznacza liczbę węzłów w elemencie (LWE). Przybliżenie ma w tym

ujęciu postać lokalnej interpolacji i jest związane ze standardowym sposobem
formułowania metody.
Dobór odpowiednich funkcji bazowych związany jest teraz z elementem,

a więc jedynie z podobszarem rozpatrywanego obszaru. Funkcje bazowe mają
zatem charakter lokalny, tzn. są różne od zera w rozważanym obszarze
i przyjmują wartość zero w pozostałych. Stwarza to możliwość bardzo
efektywnej aproksymacji rozwiązania na obszarze o praktycznie dowolnej
geometrii.
Równania MES elementu dla odpowiednich metod można uzyskać
w następujący sposób:
• dla metody Ritza – z warunku stacjonarności funkcjonału, do którego
podstawiono rozwiązanie przybliżone (3.3)
• dla metody Galerkina – stosując następujący algorytm:
1) mnożenie równań różniczkowych przez tzw. funkcje testowe
2) całkowanie po obszarze
3) zastosowanie całkowania przez części w celu zredukowania maksy-
malnego stopnia pochodnej w równaniu
4) wprowadzenie warunków brzegowych do sformułowania
Dyskretyzując w dalszej kolejności poszukiwaną wielkość zgodnie
z zależnością (3.3) oraz przyjmując takie same funkcje testowe i bazowe,
w rezultacie na podstawie uzyskanej formy słabej otrzymuje się równania MES
w ujęciu Galerkina. Wariant ten ma charakter ogólniejszy niż podejście Ritza,
pozwala bowiem rozwiązywać także takie zagadnienia, dla których nie istnieje
sformułowanie wariacyjne
3.6. Funkcje kształtu

Chcąc znaleźć zależność funkcyjną między na przykład dwiema
wielkościami przy znajomości n wartości tych wielkości, wykorzystuje się
metody aproksymacyjne. Przy skomplikowanym rozkładzie poszukiwanej
funkcji podwyższa się stopień wielomianu aproksymacyjnego, uzyskując coraz
to lepsze „dopasowanie” krzywej do wartości danych. Wzrost jednak stopnia
wielomianu powoduje, że w pewnym momencie pojawiają się silne oscylacje
krzywej i funkcja uzyskana w wyniku aproksymacji znacząco odbiega od
rozkładu rzeczywistego. Korzystniej jest wówczas aproksymować funkcję nie
w całym zadanym obszarze, ale w podobszarach.
Niech u (x) będzie poszukiwaną funkcją jednej zmiennej, a U 1 ÷ U 5
znanymi jej wartościami w przedziale 〈 0, l 〉 . Dzieląc obszar określoności funkcji
na podobszary i dokonując w nich aproksymacji funkcji wielomianami o niskim
stopniu (na przykład liniowymi) uzyskuje się końcową aproksymację jak na
rys. 3.6. Takie jednak podejście powoduje, że na granicy podobszarów
(elementów) występują nieciągłości aproksymowanej funkcji. Nieciągłości te
3.6. Funkcje kształtu 63
można wyeliminować wykorzystując na przykład funkcje sklejane. Aby jednak

wyznaczyć niewiadome współczynniki tych funkcji, należałoby rozwiązać układ
ośmiu równań; ogólnie dla n przedziałów układ 2n równań w przypadku
aproksymacji liniowej.
Rys. 3.6. Poszukiwana funkcja u (x) aproksymowana wielomianami liniowymi uˆ ( x)

w czterech podobszarach
Zdecydowanie prostszym sposobem jest opis funkcji û w obszarze elementu

przy użyciu wartości u na jego końcach, zgodnie z poniższą zależnością:
u e ≈ uˆ e = N 1 (ξ ) u (ξ = 0) + N 2 (ξ ) u (ξ = 1) , (3.4)
gdzie ξ jest współrzędną lokalną (lokalnym układem odniesienia), związaną

z danym podobszarem (elementem), ui są wartościami węzłowymi funkcji u
w elemencie, natomiast N i są pewnymi funkcjami omówionymi w dalszych
rozważaniach.
Pomiędzy lokalnym układem odniesienia a układem globalnym
(wspólnym dla wszystkich elementów) zachodzi zależność, którą obrazuje
rys. 3.7. Zależność tę opisują wzory:
x(ξ ) = xi + li ξ
1 (3.5)
ξ ( x) = ( x − xi )
li
Z kolei rozkład funkcji û e w układzie lokalnym przedstawia rys. 3.8.

Występujące w wyrażeniu (3.4) oraz na rys. 3.8 wielkości u (ξ = 0) i u (ξ = 1)

odpowiadają dla kolejnych elementów wartościom globalnym: dla elementu
pierwszego U 1 i U 2 , dla elementu drugiego U 2 i U 3 itd. W układzie lokalnym
odpowiadają wartościom węzłowym u1 i u 2 .
Rys. 3.7. Związek między układem lokalnym a globalnym
Rys. 3.8. Aproksymacja funkcji û w obszarze elementu przy użyciu wartości

węzłowych oraz funkcji N1 i N 2
W ten sposób uzyskuje się aproksymację poszukiwanej funkcji u (x) ,

która na granicy elementów spełnia warunki ciągłości (rys. 3.9).
Funkcje N i (ξ ) zwane są funkcjami kształtu (opisują kształt rozwią-
zania w obszarze elementu) i odgrywają w MES zasadniczą rolę. Charakter
uzyskanego rozwiązania i stopień aproksymacji rozwiązania dokładnego zależą
nie tylko od wymiaru i liczby elementów, ale znacząco od przyjętych funkcji
kształtu. W układzie globalnym funkcje te zobrazowano na rys. 3.9. Przy ich
wykorzystaniu globalny rozkład funkcji uˆ ( x) można przedstawić w następującej
formie:
N
uˆ ( x) = ∑ N i ( x) U i , (3.6)
i =0
co stanowi podstawowy związek metody elementów skończonych. W zale-

żności (3.6) przez N oznaczono liczbę węzłów układu (LWU).
Rys. 3.9. Końcowa aproksymacja funkcji u (x) w obszarze 〈 0, l 〉 przy użyciu

wartości węzłowych U 1 ÷ U 5 i funkcji N i (x) , i = 0,1,..., N = LWU = 5
W zapisie lokalnym (rys. 3.8) funkcje N1 (ξ ) i N 2 (ξ ) określone są wzorami:
N1 (ξ ) = 1 − ξ
(3.7)
N 2 (ξ ) = ξ
Funkcje kształtu nie są zwykle definiowane bezpośrednio. Zamiast postaci

funkcji kształtu określa się tzw. ansatz funkcje, najczęściej w postaci
wielomianów liniowych lub kwadratowych. W rozważanym przypadku funkcja
ansatz może być tylko wielomianem liniowym o następującej postaci:
uˆ e (ξ ) = aξ + b (3.8)
Podstawiając do zależności (3.8) współrzędne punktów węzłowych elementu

otrzymuje się:
u (0) = u1 = b
(3.9)
u (1) = u 2 = a + b
Stąd natomiast wynika:
b = u1
(3.10)
a = u 2 − u1
Zatem:
uˆ e = (u 2 − u1 ) ξ + u1 = u1 (1 − ξ ) + u 2ξ = N1 (ξ ) u1 + N 2 (ξ ) u 2 (3.11)
W zapisie macierzowym zależność powyższa przybiera formę:
u 
uˆ e (ξ ) = [N1 N 2 ]  1  = Nu (3.12)
u 2 
Na podobnej zasadzie konstruuje się funkcje kształtu przy aproksymacji

wielomianami wyższych rzędów. Na rys. 3.10 przedstawiono aproksymację
uˆ e (ξ ,η ) dla funkcji ansatz będącej wielomianem kwadratowym.
Rys. 3.10. Element trójkątny i przyjęta aproksymacja funkcji uˆ e (ξ ,η )

wielomianem kwadratowym
Funkcje kształtu posiadają pewne charakterystyczne własności:

• każdy węzeł i posiada związaną z nim funkcję N i
• N i = 1 w węźle i
• N i jest niezerowa tylko w węzłach połączonych z węzłem i
• suma funkcji kształtu elementu w jego dowolnym punkcie równa jest
jedności
Funkcje kształtu powinny także spełniać warunek ciągłości i zupełności:
• ciągłość: ciągłość C n−1 na granicach między elementami
• zupełność: ciągłość C n poszukiwanej wielkości w obszarze elementu,
gdzie n jest najwyższym rzędem pochodnej występującej w równaniu
całkowym.
Zapewnienie ciągłości jest ważne zarówno z fizycznego jak i matematycznego
punktu widzenia. W pierwszym przypadku chodzi o uniemożliwienie powsta-
wania szczelin między elementami i ich nakładek. Można wykazać, że przy
spełnieniu warunków ciągłości i zupełności rozwiązanie MES zbiega się do
rozwiązania ścisłego w miarę zagęszczania siatki podziału na elementy
(w wersji h metody).
Najczęściej funkcje kształtu przyjmuje się w postaci wielomianów
Lagrange’a, Hermite’a lub uzyskanych na podstawie trójkąta Pascala. Na
rys 3.11 przedstawiono przykładowo liniowe funkcje kształtu Lagrange’a.
Rys. 3.11. Liniowe funkcje kształtu Lagrange’a:

a) w układzie globalnym, b) w układzie lokalnym
Funkcje te opisane są następującymi zależnościami:
 x − xi −1
 xi −1 ≤ x ≤ xi
 xi − xi −1
x −x
N i ( x) =  i +1 xi ≤ x ≤ xi +1 (3.13)
 xi +1 − xi
 w pozostałych przedziałach
 0

Ogólnie dla wielomianu n -tego stopnia:
n  (x − x j ) 
N i ( x) = ∏   , i = 1,2,..., ( LW + 1) (3.14)
j =1  ( x i − x j ) 
 
i≠ j
Kwadratowe funkcje kształtu Lagrange’a zobrazowano z kolei na rys. 3.12.
Rys. 3.12. Kwadratowe funkcje kształtu Lagrange’a:

a) w układzie globalnym, b) w układzie lokalnym
W przypadku zagadnienia dwuwymiarowego, dla elementu prostokątnego

funkcje kształtu Lagrange’a dla i -tego węzła oblicza się z zależności:
N i ( x, y ) = N I ( x ) N J ( y ) (3.15)
gdzie:
Wx Wy
x − xj y − yj
NI = ∏x
j =1 i − xj
, NJ = ∏y
j =1 i − yj (3.16)
j≠I j≠ J
są wielomianami dla linii poprowadzonych przez węzeł i , j (rys. 3.13)

i równoległych do osi układu odniesienia.
Rys. 3.13. Element prostokątny z liczbą węzłów Wx w kierunku osi x

i W y w kierunku osi y
Funkcje kształtu Lagrange’a dla elementu trójkątnego obrazuje rysunek

3.14.
Rys. 3.14. Liniowe funkcje kształtu Lagrange’a dla elementu trójkątnego:

a) w układzie lokalnym, b) w układzie globalnym
Omówione funkcje kształtu Lagrange’a zapewniają ciągłość poszukiwanej

wielkości (ciągłość C 0 ) na granicy elementów. Niezależnie jednak od stopnia
przyjętego wielomianu ciągłość pochodnej nie może być osiągnięta. W przy-
padkach, w których zachodzi konieczność ciągłości także pochodnej stosuje się
funkcje kształtu oparte na wielomianach Hermite’a.
Przedstawione funkcje kształtu są tzw. standardowymi funkcjami kształtu,
związanymi z interpolacją w obszarze elementu przy wykorzystaniu wartości
węzłowych poszukiwanej wielkości. Ich stosowanie ograniczone jest do wersji
tradycyjnej h MES, w której do aproksymacji wykorzystuje się najczęściej
wielomiany niskiego stopnia (liniowe lub kwadratowe), a poprawę dokładności
rozwiązania uzyskuje się przede wszystkim poprzez zagęszczanie (sukcesywne)
siatki podziału. Poprawę tę można również uzyskać poprzez zwiększenie rzędu
wielomianu aproksymacyjnego, co w wersji h wiąże się z wprowadzeniem
dodatkowych węzłów w elemencie, a co za tym idzie rekonstrukcją istniejących
już przy mniejszej liczbie węzłów funkcji kształtu.
Na innej niż standardowe funkcje kształtu koncepcji oparte są natomiast
tzw. hierarchiczne funkcje kształtu. Stosowane są one w wersji p MES. W
wersji tej poprawę dokładności rozwiązania uzyskuje się przez zwiększanie
(sukcesywne) rzędu aproksymacji poszukiwanej wielkości w obszarze elementu,
przy stałej siatce podziału na elementy. Tak więc w odróżnieniu od wersji h ,
w wersji p zakłada się, że element siatki pozostaje niezmieniony (zarówno co
do wymiaru jak i usytuowania), natomiast rząd aproksymacji wewnątrz każdego
elementu może ulegać zmianie.
W pierwszej kolejności przyjmuje się aproksymację rozwiązania rzędu p
(na przykład p = 1 , aproksymacja liniowa) przy użyciu standardowych funkcji
kształtu. Poprawę dokładności rozwiązania otrzymuje się przez dodanie do
istniejących, nowych funkcji kształtu w taki sposób, że wzrost z p do
( p + 1) rzędu aproksymacji nie zmienia ani funkcji kształtu dla rzędu p , ani
stopni swobody elementu (a więc nie zmienia liczby węzłów w elemencie).
Funkcje kształtu spełniające ten warunek nazywane są hierarchicznymi.
Podlegają one zależności:
Φ ( p =1) ⊂ Φ ( p =2 ) ⊂ ... ⊂ Φ ( p =∞ ) , (3.17)
gdzie Φ ( p ) oznacza zbiór funkcji kształtu rzędu p .

Używając hierarchicznych funkcji kształtu otrzymuje się aproksymację
nieznanego rozwiązania w formie analogicznej jak w przypadku standardowych
funkcji kształtu (por. zależność 3.6):
N
uˆ ( x) = ∑ N i ( x) U i , (3.18)
i =0
przy czym N nie jest teraz związane – jak w przypadku wersji h – z liczbą
węzłów układu (LWU). Poza tym stopnie swobody U i odpowiadające funkcjom
kształtu wyższych rzędów nie mają obecnie znaczenia fizykalnego, to znaczy nie
mogą być interpretowane jako wartości rozwiązania w węzłach. Są one
natomiast wartościami pewnych wyższych pochodnych rozwiązania między
węzłami elementu.
Przykładowo, dla zagadnienia jednowymiarowego i elementu o dwóch
węzłach (rys. 3.15) funkcje kształtu mają postać:
1
N 1 (ξ ) = (1 − ξ )
2
1
N 2 (ξ ) = (1 + ξ )
2
1 (3.19)
p
N S (ξ ) = (ξ p − 1) dla p = 2,4,6,...
p!
1
N Sp (ξ ) = (ξ p − ξ )
p! dla p = 3,5,7,...
Rys. 3.15. Element jednowymiarowy, dwuwęzłowy o parametrach węzłowych u1 i u 2

w lokalnym układzie odniesienia − 1 ≤ ξ ≤ 1
Poszukiwaną funkcję u (x) można w obszarze elementu przedstawić następująco:
d pu
u e (ξ ) ≈ uˆ e = N 1 (ξ ) u1 + N 2 (ξ ) u 2 + ∑ N Sp , (3.20)
p = 2,... dξ p
gdzie: N1 i N 2 są standardowymi (Lagrange’owskimi, por. zależność. (3.7))

funkcjami kształtu, u1 i u 2 są parametrami węzłowymi, funkcje N Sp są
funkcjami kształtu wyższych rzędów, odniesionymi do środka elementu.
Przykładowe hierarchiczne funkcje kształtu przedstawiono na rys. 3.16, na
którym także celem porównania zobrazowano standardowe funkcje kształtu.
Wersja p z hierarchicznymi funkcjami kształtu wykorzystywana jest jako
podstawa adaptacyjnej formy MES. Podejście adaptacyjne (możliwe także
w oparciu o wersję h ) ukierunkowane jest na osiągnięcie dużej (założonej)
dokładności rozwiązania w możliwie najkrótszym czasie. Wymaga to stoso-
wania zaawansowanych technik związanych z kontrolą błędu i zautoma-
tyzowaniem procesu zmiany modelu MES w tych rejonach analizowanego
obszaru, w których błąd ten przekracza wartość założoną.
Rys. 3.16. Funkcje kształtu dla elementu jednowymiarowego: a) standardowe oparte

na wielomianach Lagrange’a, b) hierarchiczne oparte na wielomianach Legendre’a
3.7. Agregacja
Globalny układ równań MES, umożliwiający wyznaczenie dyskretnych wartości
poszukiwanej funkcji (i ewentualnie jej pochodnych) w węzłach można zapisać
w postaci:
K uˆ = Q , (3.21)
gdzie K jest tzw. macierzą sztywności układu:
3.7. Agregacja 73
 k11 k12 K k1n 

 
 k 21 k22 L k2 n 
K=  , (3.22)
 sym. O M 
 
 knn 
przy czym n jest równe liczbie stopni swobody układu (SSU); wektor
uˆ e = [u1 u 2 L u n ] T jest wektorem, którego składowe są poszukiwanymi
wartościami węzłowymi układu (niewiadomymi w węzłach), a
e
Q = [Q1e Q2e L Qne T] – wektorem obciążeń węzłowych (wielkości
znanych).
Proces agregacji stanowi aspekt strony matematycznej montażu
elementów (por. rys. 3.2). Po specyfikacji elementu znane są (por. rozdział 3.5)
macierze opisujące własności elementu: macierz sztywności elementu k e i
wektor obciążeń elementu F e . Na ich podstawie, przy zapewnieniu warunków
zgodności w węzłach (w mechanice konstrukcji dla niewiadomych w węzłach
będących uogólnionymi przemieszczeniami warunek ten dotyczy zgodności
przemieszczeń węzłowych i równowagi sił w węzłach) budowana jest globalna
macierz sztywności i globalny wektor obciążeń węzłowych dla całego
dyskretyzowanego obszaru. Zachodzi zależność:
K= ∑k
e
e
(3.23)
Q= ∑Q
e
e
Występujące w zależności (3.23) sumowanie nie ma charakteru algebraicznego,

ale umowny. Globalną macierz sztywności K oraz globalny wektor obciążeń
węzłowych Q otrzymuje się przez odpowiednie „dodawanie” macierzy
i wektora obciążeń elementów, który to proces nazywany jest agregacją.
Schematycznie proces agregacji można zobrazować na przykładzie układu
złożonego z dwóch elementów (rys. 3.17).
Rys. 3.17. Układ dwóch elementów skończonych z zaznaczoną numeracją globalną

i lokalną węzłów
Zakładając stopień swobody węzła SSW=1 (poszukiwaną funkcją jest funkcja

skalarna), wektor parametrów węzłowych układu û ma następujące składowe:
uˆ = [u1 u 2 u3 u4 u5 u6 ] T , (3.24)
gdzie u1 ÷ u 6 są wartościami węzłowymi poszukiwanej funkcji u.

Element jest specyfikowany w lokalnym układzie współrzędnych, który
przedstawiono na rys. 3.18.
Rys. 3.18. Czterowęzłowy element prostokątny w układzie lokalnym
Macierz sztywności elementu ma budowę:

3.7. Agregacja 75
 k11e e
k12 e
k13 e 
k14
 
e
 e
k 22 e
k 23 e
k 24 
k =  (3.25)
e e
 k 33 k 34 
 e 
 sym. k 44 
Odpowiednio wektor F e można zapisać następująco:
Q 1e 
 
e
Q e2 
F =  (3.26)
e
Q 3 
 e
Q 4 
[ ]
Wektor parametrów węzłowych elementu: uˆ e = u 1e u e2 u 3e u e4 T .
Globalny układ równań przedstawiający udziały poszczególnych elementów
można przedstawić schematycznie w postaci:
   u1 = u11   F 1 
    
 [ ]
k1   u2 = u2   
1
  1 2
  2
  u 3 = u 4 = u1  F 
  =  (3.27)
 [ ]
k 2

 4u = u 1
3 = u 2
2  
  2   
   u5 = u3   
   u = u2   
  6 4   
przy czym uzgodniono tu wartości w węzłach wspólnych układu zgodnie z
poniższą zależnością (por. zależność 3.2):
u11 = u1
u12 = u 2
u31 = u 22 = u 4
(3.28)
u14 = u12 = u3
u32 = u5
u 42 = u 6
gdzie indeks górny oznacza numer elementu.

Tak więc układ równań systemu będzie miał postać:
(3.29)
Ogólnie przy SSW stopniach swobody węzła, liczbie węzłów w elemencie

LWE i liczbie węzłów układu LWU globalna macierz sztywności ma wymiar
N × N , gdzie N = SSW × LWU = SSU . Jest ona symetryczna, rzadka, dodatnio
określona i dla standardowych funkcji kształtu ma budowę schematycznie
przedstawioną na rys. 3.19.
Globalna macierz sztywności w programowaniu metody elementów
skończonych budowana jest bezpośrednio na podstawie tzw. tablicy incydencji,
w której określona zostaje relacja między numeracją lokalną elementów
a numeracją globalną układu. Dla rozpatrywanego układu dwóch elementów
(rys. 3.17) tablica ta ma postać:
Tablica 3.1.
Lokalne i globalne numery węzłów elementów
Element 1 Element 2
lokalnie globalnie lokalnie globalnie
1 1 1 3
2 2 2 4
3 4 3 5
4 3 4 6
3.7. Agregacja 77
Rys. 3.19. Schematyczna postać globalnej macierzy sztywności.

LE – liczba elementów, a = SSW × LWE
Proces agregacji można wówczas zapisać:
 k11
1 1
k12 1
k14 1
k13 0 0   0 0 0 0 0 0 
 1 1 1 1   
 k 21 k 22 k 24 k 23 0 0   0 0 0 0 0 0 
e k 1 1
k 42 1
k 44 1
k 43 0 0   0 0 k112 k122 k132 k142 
K= ∑ k e =  411
 k 31
1
k 32 1
k 34 1
k 33 0
+
0   0 0 2
k 21 2
k 22 2
k 23 2
k 24


0 0 0 0 0 0   0 0 2
k 31 2
k 32 2
k33 2 
k 34
   
 0 0   
2 2 2 2
0 0 0 0 0 0 k 41 k 42 k 43 k 44
K i j = Α k ea b , Q i = Α Q ea , (3.30)
∀e ∀e
gdzie: Α – symbol procedury agregacji

i – globalny numer równania
j – globalny numer składowej wektora û
a – numer równania elementu
b – numer stopnia swobody elementu
Przykładowo można zatem wprost napisać:

K 3 4 = k 41 3 + k122
(3.31)
K 4 5 = k 223
Należy zaznaczyć, że w przypadku gdy osie układów lokalnego i

globalnego nie są równoległe, przed agregacją zachodzi konieczność
transformowania macierzy sztywności elementów i wektora obciążeń
węzłowych elementów do układu globalnego.
Bez uwzględnienia warunków brzegowych otrzymana na drodze agregacji
globalna macierz sztywności jest macierzą osobliwą, a więc układ równań
systemu nie ma rozwiązania. Istnieje zatem konieczność wprowadzenia
warunków brzegowych. Dotyczy to warunków brzegowych podstawowych
(Dirichleta, w mechanice: kinematycznych warunków brzegowych). Naturalne
warunki brzegowe (Neumanna, w mechanice: statyczne warunki brzegowe) są
włączone zawsze do funkcjonału i przy minimalizacji funkcjonału (metoda
Ritza) spełnione są w sposób automatyczny (w MES w sposób przybliżony
w wersji przemieszczeniowej). W metodzie ważonych residuów (metoda
Galerkina) naturalne warunki brzegowe wprowadzone są przy przekształcaniu
równań całkowych (twierdzenie Greena).
Uwzględnienia podstawowych warunków brzegowych w globalnej
macierzy sztywności można dokonać przez podział wektora parametrów
węzłowych układu uˆ ( N ×1) na zbiór parametrów węzłowych poszukiwanych û p
oraz zbiór parametrów znanych (zadanych) û z , określonych przez warunki
brzegowe. Przy takim podziale układ równań (3.21) systemu można zapisać w
następującej postaci:
K p p K p z  uˆ p  Q p 
K  = 
K z z  uˆ z   Q z 
(3.32)
 zp
Zakładając, że zadane są jednorodne warunki brzegowe, a więc uˆ z = 0 ,

z równań powyższych uzyskuje się:
Q p = K p p uˆ p → uˆ p = K −p1p Q p (3.33)
Proces agregacji dla elementów hierarchicznych przebiega zasadniczo tak

samo jak dla elementów o standardowych funkcjach kształtu. Macierz
sztywności elementu hierarchicznego ma budowę jak na rys. 3.20 – macierz
elementu rzędu p + 1 zawiera macierz rzędu p .
3.7. Agregacja 79
Rys. 3.20. Budowa macierzy sztywności elementu hierarchicznego rzędu np
Macierz rzędu p = 1 jest macierzą sztywności elementu dla

standardowych funkcji kształtu. Odpowiadający jej wektor niewiadomych ma
składowe będące węzłowymi stopniami swobody elementu. Pozostałe składowe
wektora niewiadomych odpowiadające hierarchicznym stopniom swobody nie są
wartościami rozwiązania w węzłach, na brzegu czy w środku elementu. Stąd też
agregacja do macierzy globalnej dotyczy głównie podmacierzy macierzy
sztywności elementu k ep odpowiadającej rzędowi p = 1 . W przypadku gdy
pochodne funkcji kształtu są ortogonalne (a tak najczęściej są one dobierane w
elementach hierarchicznych), czyli:
∫ ∇N i ⋅ ∇N j dx = 0 i≠ j (3.34)
Ω
macierz sztywności k ep zawiera elementy zerowe (rys. 3.21).

Macierz ta jest prawie diagonalna, jeżeli pochodne funkcji kształtu są
wielomianami Legendre’a:
N p = ∫ Pp −1 dξ , ξ (0) = ξ (1) = 0 p≥2 (3.35)
1 dp  2
gdzie: Pp =
1
( )
ξ − 1  - wielomian Legendre’a stopnia p.
p
( p − 1)! 2 p −1
dξ 
p 
Rys. 3.21. Macierz sztywności elementu przy ortogonalnych funkcjach kształtu
Globalna macierz sztywności przy stosowaniu tej samej rodziny elementów

hierarchicznych ma postać przedstawioną schematycznie na rys. 3.22.
Rys. 3.22. Globalna macierz sztywności dla elementów hierarchicznych

3.8. Przykład – równanie różniczkowe podłoża Własowa 81
3.8. Przykład – równanie różniczkowe podłoża

Własowa
Dane jest równanie różniczkowe opisujące model podłoża Własowa:
d 2 w( x )
− 2t + k w( x) = q ( x) w Ω = (0,1) , (3.36)
dx 2
gdzie: w(x) – funkcja ugięcia

q (x) – funkcja obciążenia gruntu
2t, k – współczynniki podłoża (parametry)
Dane są następujące warunki brzegowe:
w(0) = 0
(3.37)
w(1) = 0
Do rozwiązania metodą elementów skończonych jest zagadnienie

brzegowe opisane równaniami (3.36), (3.37)
Transformowanie równań do postaci słabej – metoda ważonych

residuów
Jeżeli v(x ) jest funkcją wagową, wówczas słaba forma równania (3.36) ma
postać:
1 d 2 w( x)
∫ [−2t + k w( x) − q( x)] v( x) dx = 0 (3.38)
0 dx 2
Aby otrzymać symetryczną formę słabą (dogodną do zastosowań nume-

rycznych), do pierwszego członu równania (3.38) zostanie zastosowane
całkowanie przez części:
1 d 2 w( x) 1 dw( x ) dv ( x ) dw( x)
∫ [−2t 2
v ( x ) dx = 2t [ ∫ dx − v( x) 10 ] (3.39)
0 dx 0 dx dx dx
Uwzględniając w powyższym wyrażeniu warunki brzegowe (3.37), otrzymuje

się:
1 d 2 w( x) 1 dw( x ) dv( x )
∫ [−2t v ( x ) dx = 2t ∫ dx (3.40)
0 dx 2 0 dx dx
Zależność (3.38) przybierze więc postać:
1 dw( x) dv( x) 1
∫ (2t + k w( x) v( x)) dx = ∫ q v( x) dx (3.41)
0 dx dx 0
Wariacyjna zatem (słaba) postać równania różniczkowego (3.36) jest

następująca:
Wyznaczyć w∈V , gdzie V = {v ∈ H 01 (Ω) : v(0) = 0, v(1) = 0}, taką że:
1 dw dv 1
∫ (2t + kw v)dx = ∫ q v dx ∀v ∈V (3.42)
0 dx dx 0
H 1 (Ω) jest przestrzenią funkcji rzeczywistych określonych na Ω , które są

całkowalne w kwadracie i których pierwsze pochodne są również całkowalne w
kwadracie. W szczególności:
H 01 (Ω) = {v : Ω → R, v ∈ H 1 (Ω) i ∀x ∈ Γ v( x) = 0}
1
dv
∫ ( dx )
2
Na funkcję v nałożone jest więc minimalne ograniczenie: dx < ∞
0
(pierwsza pochodna funkcji całkowalna w kwadracie).
Formę dyskretną zależności (3.42) można zapisać w postaci:
Wyznaczyć wh ∈ Vh , gdzie Vh ⊂ V , taką że:
1 dwh dv h 1
∫ (2t + k wh v h ) dx = ∫ q v h dx ∀vh ∈ Vh (3.43)
0 dx dx 0
Rozwiązania przybliżonego poszukuje się zatem w przestrzeni skończenie

wymiarowej Vh , która jest podprzestrzenią przestrzeni rozwiązań dokładnych.
Przyjmuje się funkcję próbną wh w postaci:
LWU
wh ( x) = ∑ ai N i ( x), (3.44)
i =1
gdzie LWU jest liczbą węzłów układu.

Zgodnie z podejściem Galerkina, funkcja wagowa opisana jest ana-
logicznie:
LWU
v h ( x) = ∑ b j N j ( x), (3.45)
j =1
gdzie N i ( x) są funkcjami bazowymi (kształtu), które należy dobrać.

Podstawiając reprezentacje (3.44) i (3.45) do zależności (3.43) otrzymuje się:
1 LWU dN i LWU dN j LWU LWU

∫ [2t ( ∑ a i )( ∑ bj ) + k ( ∑ a i N i ) ( ∑ b j N j )] dx =
0 i =1 dx j =1 dx i =1 j =1
(3.46)
1 LWU
= ∫ ( ∑ b j N j ) dx
0 j =1
Po przekształceniu powyższa zależność przybierze postać:
LWU 1 dN i dN j 1
∑ a i ∫ (2t + k N i N j ) dx − ∫ qN j dx = 0,
i =1 0 dx dx 0 (3.47)
j = 1,2,..., LWU
Przyjmując oznaczenia:
1 dN i dN j 1
K ij := ∫ 2t ( + k N i N j ) dx , Qi := ∫ q N j dx ,
0 dx dx 0 (3.48)
związek (3.47) można zapisać wzorem:
LW
∑ K j i ai = Q j , (3.49)
i =1
lub przedstawić w postaci macierzowej:
Ka = Q (3.50)
Dyskretyzacja MES
Obszar Ω = (0,1) zostaje podzielony na N − 1 podobszarów (elementów)

0 = x0 < x1 < ... < x N −1 < x N = 1 zgodnie z rys. 3.23.
Rys. 3.23. Podział obszaru (0,1) na elementy skończone
Dobór przestrzeni skończenie wymiarowej Vh , w której poszukuje się

rozwiązania przybliżonego
Przyjmuje się jako przestrzeń Vh , przestrzeń kawałkami ciągłych funkcji

liniowych:
Vh = {wh : [0,1] → R : wh jest ciągła, wh [ xk , xk +1 ] jest liniowa,

(3.51)
k = 0,..., LW i wh (0) = wh (1) = 0}
Przestrzeń Vh jest przestrzenią skończenie wymiarową, a zatem posiada bazę.

Baza tej przestrzeni {N1 , N 2 ,..., N LW } dla przyjętej dyskretyzacji określona jest
następująco:
 x − x j −1

 x j − x j −1 x j −1 ≤ x ≤ x j


N j ( x) =  x j +1 − x (3.52)
x − x
 j +1 j x j ≤ x ≤ x j +1

 0 w pozostałych przedziałach
Są to więc liniowe funkcje kształtu Lagrange’a (por. zależność (3.13)) nazywane

też funkcjami „kapeluszowymi”.
Każdą funkcję należącą do Vh można przedstawić jako liniową kombi-
nację funkcji bazowych i współczynników, w szczególności rozwiązania
przybliżonego poszukuje się zgodnie z zależnością (3.45).
Dla elementu „e” (rys. 3.23) zależność ta przybiera postać:
whe = a1 N1e + a 2 N 2e , (3.53)
gdzie:
x j +1 − x x − xj
N1e = , N 2e =
x j +1 − x j x j +1 − x j
są lokalnymi funkcjami kształtu dla elementu zawierającego węzły x j i x j +1 .

Nakładając na rozwiązanie whe wymóg x = x j = wj i x = x j +1 = w j +1 ,
otrzymuje się: a1 = w j oraz a 2 = w j +1 . Stąd rozwiązanie w obszarze elementu
przybiera postać:
whe = w j N1e + w j +1 N 2e (3.54)
Równania MES dla elementu
Biorąc pod uwagę związek określający K i j (zależność (3.46)), występującą tam

całkę można przedstawić w postaci sumy całek w obszarze elementów:
1 dN i dN j x2 dN dN
j
K i j = ∫ 2t ( + k N i N j )dx = 2t [ ∫ ( i ) +K
0 dx dx 0 dx dx
, (3.55)
x LW dN i dN j
K+ ∫ ( )dx]
x LW −1 dx dx
gdzie i , j są numerami węzłów.

Podobnie:
1 x2 x LW
Q j = ∫ q N j dx = ∫ q N j dx + ... + ∫ q N j dx (3.56)
0 0 x LW −1
Uwzględniając lokalne funkcje kształtu elementu N ie , zależność (3.55) i (3.56)

można przedstawić następująco:
e
dN ie dN j
e e
Ki j =∑ ∫ dx = ∑ k iej
Ω e dx dx
(3.57)
e e
Qj = ∑ ∫ q e
N ej dx = ∑ Q ej
Ωe
W zależności powyższej przyjęto oznaczenia:

e
dN ie dN j
k iej
= ∫ dx
Ω e dx dx , (3.58)
e e e
Q j = ∫ q N j dx
Ωe
gdzie k iej i Q ej są odpowiednio macierzą sztywności i wektorem obciążeń

elementu.
Zależność (3.47) można zatem zapisać wzorem:
LW e e
∑ ai (∑ k i j ) = ∑ Q j
e e
(3.59)
i =1
Żądając, aby wartości funkcji wh w kolejnych węzłach układu (1,2,..., LWU )

wynosiły w1 , w2 ,...wLW , czyli:
wh ( x = x1 ) = w1
wh ( x = x2 ) = w2
(3.60)
......................
wh ( x = x LW ) = wLW
otrzymuje się:
a1 = w1
a 2 = w2
(3.61)
...........
a LW = wLW
Zatem zależność (3.60) przybierze postać:
LW e e
∑ wi (∑ k i j ) = ∑ Q j
e e
(3.62)
i =1
W zapisie macierzowym związek powyższy przedstawia się następująco:
Kw =Q , (3.63)
e
gdzie: K = ∑ k iej – globalna macierz sztywności, w = [ w1 , w2 ,..., wLW ]T –
e
wektor parametrów węzłowych układu, Q = ∑Q e
j – globalny wektor obciążeń
węzłowych.
Zależność (3.62) odniesiona do elementu ma postać:
LW e
∑ wie k iej = Q ej , (3.64)
i =1
gdzie: LW e – liczba węzłów elementu.

W zapisie macierzowym otrzymuje się związek:
k e w e = Qe , (3.65)
gdzie: k e – macierz sztywności elementu, w e – wektor parametrów węzłowych

elementu, Q e – wektor obciążeń węzłowych elementu.
Dla przyjętego elementu liniowego, dwuwęzłowego, zależność powyższa
przedstawia się następująco:
 k11
e e  e
k12 w1 Q1e 
 e e   e  =  e, (3.66)
k 21 k 22  w2  Q2 
gdzie wyraz k iej (i,j=1,2) (e=1,2,…,N-1) macierzy sztywności oraz składowe

Q1e i Q2e wektora obciążeń opisane są zależnością (3.58).
Równania (3.66) są równaniami MES elementu. Dla konkretnej
dyskretyzacji równania te określa się dla każdego elementu. Są one podstawą do
utworzenia równań układu (3.63). Postępując zatem w dalszej kolejności według
algorytmu przedstawionego w rozdziale 3.2, otrzymuje się rozwiązanie MES
rozpatrywanego zagadnienia.
3.9. Ogólne sformułowanie MES dla zagadnień

mechaniki ciała stałego
W mechanice ciała stałego metoda elementów skończonych może być
formułowana w ujęciu przemieszczeniowym, naprężeniowym lub mieszanym
(hybrydowym). Powszechnie stosowana jest MES w ujęciu przemieszcze-
niowym i można ją uznać za podstawę w zagadnieniach mechaniki konstrukcji.
W wersji tej, aproksymacji w obszarze elementu podlega pole przemieszczeń,
stąd niewiadomymi w węzłach (parametrami węzłowymi) są uogólnione
przemieszczenia. Pole naprężeń określane jest na podstawie obliczonego pola
przemieszczeń. Równania równowagi i warunki brzegowe w naprężeniach
spełnione są tylko w przybliżeniu (dla całego układu).
W naprężeniowej wersji MES (stosowanej na przykład do zagadnień
opisanych funkcjonałem Reissnera) parametrami węzłowymi są naprężenia. Do
aproksymacji przyjmuje się tylko takie pola naprężeń w elemencie, które
a priori spełniają równania równowagi.
Równania MES dla zagadnień mechaniki ciała stałego formułuje się na
ogólnych zasadach, przedstawionych w rozdziale 3.3. Dla wersji przemiesz-
czeniowej MES, będącej przedmiotem rozważań, równania te można uzyskać
w sposób następujący:
3.9. Ogólne sformułowanie MES dla zagadnień mechaniki 89
• w podejściu bezpośrednim, wykorzystując definicję macierzy sztywności

używaną w metodzie przemieszczeń
• w oparciu o zasadę prac wirtualnych. W przypadku stosowania zasady
prac wirtualnych ma się do czynienia ze słabym sformułowaniem równań
równowagi. Mnożąc bowiem równania konstytutywne przez funkcje
testowe, które w tym przypadku są przyjmowane jako wariacje funkcji
przemieszczeń interpretowane jako przemieszczenia wirtualne, przy
wykorzystaniu formuły Greena (całkowanie przez części) otrzymuje się
formę całkową, która w mechanice znana jest jako zasada prac wirtualnych
• w oparciu o energię potencjalną układu i sformułowanie problemu
minimalizacji funkcjonału na bazie twierdzenia Lagrange’a o minimum
energii potencjalnej układu.
Równania równowagi MES dla zagadnień mechaniki zostaną przed-
stawione przy wykorzystaniu ostatniego z wymienionych podejść. Ten też
sposób będzie wykorzystany w rozważanych w rozdziale 5 przykładach.
W celu sformułowania równań rozważa się w kartezjańskim układzie
współrzędnych xyz ciało o objętości V i brzegu S , ograniczone więzami
i poddane działaniu sił masowych i powierzchniowych (rys. 3.24).
Rys. 3.24. Ciało trójwymiarowe poddane działaniu sił masowych i powierzchniowych
W zapisie równań zostanie wykorzystana notacja macierzowa, bardziej niż

wskaźnikowa dogodna przy formułowaniu algorytmów numerycznych.
W notacji tej tensory rzędu drugiego – tensor naprężenia i odkształcenia są
odwzorowane przez wektory kolumnowe, których składowe są niezależnymi
składowymi tensorów, umieszczonymi w dowolnej (jednak konsekwentnej)
kolejności.
Odkształcając się pod wpływem powoli wzrastających obciążeń ciało

gromadzi pewien zasób energii zwanej energią potencjalną. Dla ciała
przedstawionego na rys. 3.24 wyrażenie na całkowitą energię potencjalną Π c
może być przedstawione w postaci:
1 T ~ T f dS ,
2 V∫ ∫u
Πc = σ ε dV − f B dV − ∫ u
T
S
14243 144424
V S
443 (3.67)
ENERGIA DEFORMACJI POTENCJAŁ SIŁ ZEWNETRZNYCH
SPREZYSTEJ
gdzie: σ = [σ x , σ y , σ z ,τ y z ,τ z x ,τ x y ]T – wektor stanu naprężenia,

ε = [ε x , ε y , ε z , γ y z , γ z x , γ x y ]T – wektor stanu odkształcenia,
u = [u x , u y , u z ]T – wektor stanu przemieszczenia,
~ = [u~ , u~ , u~ ]T – wektor stanu przemieszczenia powierzchni,
u x y z
f B = [ f Bx , f By , f Bz ]T – wektor sił masowych,

f S = [ f Sx , f S y , f Sz ]T – wektor sił powierzchniowych.
W warunkach równowagi statycznej ciała przy kinematycznie
dopuszczalnych polach przemieszczenia, funkcjonał energii potencjalnej osiąga
minimum. Warunkiem koniecznym istnienia ekstremum funkcjonału jest
zerowanie się wariacji funkcjonału zgodnie z zależnością:
δ Πc = 0 (3.68)
W celu wykorzystania funkcjonału (3.67) do wyprowadzenia równań MES

dyskretyzuje się obszar V , dzieląc go na N nie pokrywających się elementów.
Funkcjonał (3.68), jako skalar może być wówczas przedstawiony w postaci
sumy odpowiednich składników:
Π c = Π c(1) + Π (c2) + ... + Π c( N ) (3.69)
[ ]
Π e u e otrzymuje się odnosząc wektory występujące w zależności (3.67)
do elementu „e” i przyjmując standardową aproksymację MES pola przemie-
szczenia w obszarze elementu:
u e ( x , y , z ) = N e ( x, y , z ) δ e , (3.70)
gdzie: u ie ∈ H 1 (V e ), (i = x, y, z ) oraz:
N e ( x, y, z ) – macierz funkcji kształtu o wymiarze (3 × SSE ) ,
δ e = [δ 1e , δ 2e ,..., δ SSE
e
]T – wektor parametrów przemieszczeń węzłowych
elementu.
Związki geometryczne i fizyczne liniowej teorii sprężystości odniesione
do elementu mają postać:
εe = L ue
, (3.71)
σ e = Deε e
gdzie L i D są macierzami odpowiednio operatorów różniczkowych i modułów

sprężystości.
Uwzględniając w zależności (3.71) 1 związki (3.70) otrzymuje się:
ε e = L N eδ e = B e δ e , (3.72)
gdzie B e jest macierzą odkształceń (opisuje odkształcenia w każdym punkcie

elementu, spowodowane jednostkowym przemieszczeniem kolejnych stopni
swobody węzłów).
e
Podstawiając do funkcjonału (3.67) odniesionego do elementu Π c
związki (3.71) 2 i (3.72), otrzymuje się:
1 ~e T e e
Π ec = e T e e e T e e
∫ (σ ) ε dV − ∫ (u ) f B dV − ∫ (u ) f S dS =
2 Ve V e S e
1 e e e T e e e e e T e e
= ∫ (D B δ ) B δ dV − ∫ (N δ ) f B dV −
2 Ve Ve
(3.73)
~ 1
− ∫ (N e δ e ) T f Se dS e = ∫ (δ e ) T (B e ) T D e B e δ e dV e −
Se
2Ve
~
− ∫ (δ e ) T (N e ) T f Be dV e − ∫ (δ e ) T (N e ) T f Se dS e
Ve Se
Gdy na dyskretyzowany układ działają siły przyłożone w węzłach

P = [P1 P2 K PLWU ] T wyrażenie (3.69) można zapisać w postaci:
Π C = (∑ Π Ce ) − δ T P (3.74)
e
gdzie δ = [δ 1 δ 2 K δ SSU ] T – wektor przemieszczeń węzłowych układu.

Warunek (3.68) można obecnie zapisać w następującej postaci:
∂ ΠC ∂ ΠC ∂ ΠC
δ ΠC = δ δ1 + δ δ2 +K+ =0 (3.75)
∂δ 1 ∂δ 2 ∂δ SSU
Równość powyższa zachodząca dla dowolnych wariacji δδ 1 , δδ 2 ,..., δδ SSU

prowadzi do układu równań:
 ∂Π c 
 
∂δ
∂Π c  1 
= M =0 (3.76)
∂δ  ∂Π c 
 
 ∂δ SSE 
Ponieważ Π C jest funkcjonałem formy kwadratowej można zapisać pochodne

dla elementu „e” w następującej postaci:
∂ Π Ce
e
= k e δe + Fe (3.77)
∂δ
Minimalizujący układ równań (3.76) można zatem zapisać w postaci:
∂ ΠC
=Kδ + F =0 (3.78)
∂δ
gdzie:
K = ∑ke
e
(3.79)
F = −( ∑ F e + P )
e
Na podstawie zależności (3.73) otrzymuje się więc:
∂Π ec
e
= ∫ [(B e ) T D e B e ]δ e dV e − ∫ (N e ) T f Be dV e −
∂δ Ve Ve
~
− ∫ (N e ) T f Se dS e = k e δ e + F e
Se
gdzie:
k e = ∫ (B e ) T D e B e dV e jest macierzą sztywności elementu (3.80)
Ve
~
F e = ∫ [(N e ) T f Be dV e + ∫ (N e ) T f Se ] dS e – wektor równoważnych sił
e e (3.81)
V S
węzłowych od obciążeń masowych i powierzchniowych elementu
Równania (3.78) są równaniami MES równowagi układu .

Warto zwrócić jeszcze uwagę na aproksymację pola przemieszczenia
opisaną zależnością (3.70). Występujące w niej funkcje kształtu powinny, jak to
zostało omówione w rozdziale 3.6, spełniać warunek ciągłości i zupełności.
W odniesieniu do rozważanego problemu warunek ciągłości oznacza, że
przemieszczenia wewnątrz elementu i na jego brzegach powinny być ciągłe.
Dotyczy to przemieszczeń uogólnionych, a więc zarówno przesuwów
(translacyjne stopnie swobody) – ciągłość klasy C 0 , jak również obrotów
(rotacyjne stopnie swobody) – ciągłość klasy C 1 (jeżeli kąty obrotu wyrażają się
jako pochodne funkcji przemieszczeń).
Warunek zupełności wiąże się natomiast z wymogiem, aby funkcje
przemieszczeń elementu mogły reprezentować jego ruch sztywny oraz stan
stałych odkształceń. Elementy spełniające powyższe kryteria nazywane są
elementami dostosowanymi, w przeciwieństwie do elementów niedostoso-
wanych, dla których nie jest spełniony warunek ciągłości (na przykład niektóre
elementy płytowe), a które także są z powodzeniem stosowane.
Rozdział 4
CHARAKTERYSTYKA WYBRANYCH
ELEMENTÓW SKOŃCZONYCH
4.1. Uwagi wstępne

W rozdziale 3.4. omówiono ogólną klasyfikację elementów skończonych.
Przedstawiono również podział elementów z punktu widzenia zastosowań
w mechanice konstrukcji. Podział ten obejmował elementy podstawowe, ciągłe,
specjalne i makroelementy. W niniejszym rozdziale scharakteryzowane będą
(specyfikowane) przykładowe elementy skończone z poszczególnych, wyżej
wymienionych grup. Odpowiednio omówione będą: element belkowy,
izoparametryczny element płyty wyższego rzędu, element nieskończony oraz
makroelement podłoża warstwowego. Elementy te stanowić będą bazę
w rozważanych w rozdziale 5 przykładach.
4.2. Element belkowy

Rozważa się prostoliniową belkę o długości l, stałym przekroju
poprzecznym A i module sprężystości E, zginaną w płaszczyźnie xoy (rys. 4.1).
Zakłada się, że belka spełnia założenia teorii Eulera-Bernoulliego, a więc pomija
się wpływ sił poprzecznych na odkształcenia.
Równanie osi odkształconej belki ma znaną postać:
96 Charakterystyka wybranych elementów skończonych
d 4 w( x)
EJ = q ( x) (4.1)
dx 4
gdzie: J – moment bezwładności, w(x) – funkcja ugięcia (przemieszczeń

pionowych)
Rys. 4.1. Belka poddana działaniu obciążenia rozłożonego q (x) , prostopadłego do osi
obojętnej: a) stan deformacji, b) uogólnione siły i przemieszczenia brzegowe
Związki liniowej teorii sprężystości (3.75) przybiorą obecnie dla całego układu
postać (jednoosiowy stan naprężenia):
du ( x)
ε x = Lu =
dx (4.2)
σ x = D ε = E ( x)ε x
gdzie u (x) jest funkcją przemieszczeń poziomych związaną z funkcją ugięcia

poprzez zależność:
dw( x)
u ( x) = y (4.3)
dx
Energia deformacji sprężystej belki (por. zależność (3.68)) przyjmuje wobec

związków (4.2) i (4.3) postać:
4.2. Element belkowy 97
1 T 1 T 1l 2
Πd = ∫ σ ε dV = ∫ ε D ε dV = ∫ ∫ Eε x dA dx =
2V 2V 20A
2 2 (4.4)
1 l  d 2 w( x)  1l  d 2 w( x) 
= ∫ E  ∫ y 2
dA = ∫ EJ   dx
2 0  dx 2  A 2 0  dx 2 
l
gdzie: J = ∫ y 2 dA .
0
Przy założeniu stałego na długości belki modułu sprężystości podłużnej
E , czyli E ( x) = const. , zależność (4.4) przybierze postać:
2
EJ l  d 2 w( x) 
Πd = ∫  2
 dx

(4.5)
2 0  dx 
Potencjał sił zewnętrznych (zależność. 3.68) można zapisać następująco:
~ T f dS
L = −∫ u (4.6)
S
S
W analizowanym zagadnieniu wyraża się on następująco:
l
L = − ∫ q ( x) w( x) dx (4.7)
0
Tak więc całkowitą energię potencjalną rozpatrywanej belki opisuje zależność:
2
EJ
l
 d 2 w( x)  l
Πc = Πd + L =
2 ∫0  dx 2  dx − ∫0 q( x) w( x) dx
  (4.8)
Funkcjonał energii potencjalnej belki można również zapisać w rozszerzonej

postaci (rys.4.1):
2
EJ
l
 d 2w  l
ΠC = Πd + L =
2 ∫0  dx 2  dx − ∫0 q( x) w( x) dx + (T0 w0 −
 
(4.9)
dw0 dwl
− M0 ) + (−Tl wl + M l )
dx dx
Związki (4.8) lub (4.9) stanowią podstawę do uzyskania równań równowagi

MES elementu w oparciu o twierdzenie Lagrange’a o minimum energii
potencjalnej układu i przy zastosowaniu podejścia Ritza (por. rozdział 3.9).
Równania te można również uzyskać metodą ważonych residuów, analogicznie
jak to przedstawiono w przykładzie 3.8. Analizując zwięźle ten sposób
otrzymania równań, ogranicza się rozważania do wyprowadzenia formy słabej
równania (4.1).
Dla funkcji wagowej v(x) słaba forma równania (4.1) ma postać (por.
rozdział 3.8):
l  d 4 w( x) 
∫  EJ − q ( x )  v( x) dx = 0 (4.10)
0 dx 4 
Całkując równanie (4.10) dwukrotnie przez części otrzymuje się następujące

wyrażenie:
l d 2 w( x) d 2 v( x) l
EJ ∫ dx − ∫ q ( x) v( x) dx +
0 dx 2 dx 2 0
 d w( x)
3
d w( x) dv( x) 
2
l
(4.11)
+ EJ  3
v( x) −  0 =0
 dx dx 2 dx 
W rozpatrywanym zagadnieniu naturalne (statyczne) warunki brzegowe opisane

są zależnościami:
d 2 w( x)
EJ x=0 = M0
dx 2
momenty zginające
d 2 w( x)
EJ x =l = Ml
dx 2
(4.12)
d 3 w( x)
EJ x =0 = T0
dx 3
siły poprzeczne
d 3 w( x)
EJ x =l = Tl
dx 3
Warunki brzegowe podstawowe (geometryczne) mają natomiast postać:

w( x) x =0 =0
przemieszczenia pionowe
w( x) x =l =0
(ugięcia)
dw( x) (4.13)
x =0 =0
dx
kąty obrotu
dw( x)
x =l =0
dx
Uwzględniając w zależności (4.11) warunki statyczne (4.12), otrzymuje się:
l d 2 w( x) d 2 v( x) l
EJ ∫ dx − ∫ q ( x) v( x) dx +
0 dx 2 dx 2 0
(4.14)
dv( x) dv( x)
+ T0 v( x) l − T0 v( x) 0 − Ml l + M0 0 =0
dx dx
Na funkcję w(x) oraz v(x) nałożone są ograniczenia:
2
l  d 2w  l  d 2v 

∫ 2   < ∞ , ∫  2  < ∞ (4.15)
0  dx  0  dx 
czyli druga pochodna funkcji w(x) i v(x) całkowalne w kwadracie.

Tak więc słaba (wariacyjna) postać równania osi odkształconej belki (4.1)
jest następująca:
Wyznaczyć w ∈ V ,
 dv(0) dv(l )
gdzie V = v ∈ H 02 (Ω) : v(0) = = 0 , v(l ) = = 0} ,
 dx dx
taką, że:
d 2 w( x) d 2v( x)
l l
EJ ∫ 2
dx = ∫ q ( x) v( x) dx −Tl v( x) l −
0 dx dx 2 0
(4.16)
dv( x) dv( x)
− T0 v( x) 0 − M l l + M0 0 , ∀v ∈ V
dx dx
przy czym Ω = (0, l ) , H 2 (Ω) jest przestrzenią funkcji rzeczywistych okre-

ślonych na Ω , które są całkowalne w kwadracie aż do drugiej pochodnej.
Zależność (4.16) można zapisać w następującej formie:

Wyznaczyć w ∈ V , taką że:
B ( w, v) = l (v) , ∀ v ∈V (4.17)
gdzie:
l d 2 w d 2v
B ( w, v) = EJ ∫
0 dx 2 dx 2
(4.18)
l dv dv
l (v) = ∫ q v dx − (Tl v l − T0 v 0 − Ml l+ M0 0)
0 dx dx
Wobec symetrycznej formy biliniowej B ( w, v) można zdefiniować funkcjonał:

2
1
l
 d 2 w( x)  l
F ( w) = EJ ∫   dx − ∫ q ( x) w( x) dx +Tl w( x) l −
2 0 dx 2  0 (4.19)
dw( x) dw( x)
− T0 w( x) 0 − M l l+ M0 0
dx dx
Funkcjonał (4.19) odpowiada całkowitej energii potencjalnej rozpatrywanej

belki, opisanej zależnością (4.9).
W celu specyfikacji elementu belkowego dokonuje się dyskretyzacji
obszaru Ω = (0, l ) belki na elementy skończone zgodnie ze schematem
przedstawionym na rys. 3.3 i analizuje się wyizolowany z układu element „e”.
Element ten jest elementem prostoliniowym, dwuwęzłowym.
Z analizy funkcjonału energii potencjalnej belki (4.8), w którym występują
2 . pochodne funkcji ugięcia wynika, że aby zachowana była ciągłość funkcji
w(x) w obszarze i wzdłuż brzegu elementu, zarówno funkcje próbne wh (x)
(por. zależność (3.42)), które są aproksymacją poszukiwanej funkcji w(x) , jak
również ich pierwsze pochodne muszą być ciągłe w obszarze i wzdłuż brzegu
elementu. Potwierdza to także forma słaba (4.15), która wskazuje, że w modelu
skończenie elementowym funkcje aproksymacyjne muszą być z przestrzeni
Vh ⊂ H 2 (0, l ) . Oznacza to, że funkcje kształtu elementów muszą być klasy
C 1 (0, l ) . Aby spełnić ten wymóg, należy jako niewiadome w węźle elementu
(stopnie swobody węzła) przyjąć zarówno funkcję ugięcia w(x) jak również jej
dw( x)
pierwszą pochodną . Z drugiej strony rozpatrując model fizyczny
dx
zagadnienia, widoczne jest (rys. 4.1), że przemieszczenia dowolnego punktu
belki można opisać przez przemieszczenie pionowe (ugięcie) i kąt obrotu, który
jest pierwszą pochodną ugięcia. Tak więc w każdym węźle elementu należy
przyjąć dwa uogólnione przemieszczenia: przemieszczenie pionowe w i kąt

obrotu ϕ , uzyskując element o czterech stopniach swobody (SSE=4) (rys. 4.2).
Rys. 4.2. Element belkowy w lokalnym układzie współrzędnych − 1 ≤ ξ ≤ 1
Zależność między współrzędną globalną x i lokalną ξ (rys.4.2) wyraża się

wzorem:
x = x1 + 0,5 l e (1 + ξ ) (4.20)
Wektor przemieszczeń węzłowych elementu δ e ma następujące składowe:
[
δ e = w1 ϕ x 1 w2 ϕ x 2 ]
T
(4.21)
Przemieszczenia pionowe aproksymuje się w obszarze elementu wielomianem

3. stopnia:
w e (ξ ) = α 1 + α 2ξ + α 3ξ 2 + α 4ξ 3 (4.22)
Tak więc przemieszczenia uogólnione opisane są następująco:

 w e  1 ξ ξ2 ξ 3  α 1 
 
 e  = 0 2 4 6 2 M  = P ⋅α (4.23)
ϕ x  ξ  
 le le le  α 4 
gdzie α jest wektorem nieznanych współczynników, a P – tzw. macierzą

funkcyjną.
W wyrażeniu (4.23) wykorzystano związki:
d (•) d (•) dξ
= ⋅ (4.24)
dx dξ dx
Zależność (4.23) słuszną w obszarze elementu powinny spełniać też współ-

rzędne punktów będących węzłami elementu:
 
 1 −1 1 −1 
 w1e    α 
   2 4 6  1
0 −
ϕ x 1  
e
le le le   M 
 e =   (4.25)
 w2   1 1 1 1  M 
 e   
ϕ x 2   2 4 6  α 4 
 0 
 le le le 
czyli:
δe = C ⋅ α (4.26)
gdzie C jest tzw. macierzą współrzędnych.

Na podstawie związku (4.26) otrzymuje się:
α = C −1 δ e (4.27)
Zgodnie z podejściem Ritza przyjmuje się aproksymację poszukiwanej funkcji

w(x) w następującej postaci:
w e ( x) = N e ( x) ⋅ δ e (4.28)
Z zależności (4.23) wynika:
α 1 
 
e
w =1 ξ [ ξ 2
ξ 3
]  M  = Pw α (4.29)
α 
 4
[
gdzie: Pw = 1 ξ ξ2 ξ3 . ]
Podstawiając do związku (4.29) równanie (4.27) otrzymuje się:
w e = Pw ⋅ C −1 ⋅ δ e (4.30)
Porównanie zależności (4.28) i (4.30) prowadzi do następującego wyrażenia:
N e = Pw ⋅ C −1 (4.31)
Związek (4.31) jest podstawą do wyznaczenia macierzy funkcji kształtu

N e = [N 1 N 2 N 3 N 4 ] . Macierz odwrotna do macierzy współrzędnych C (jeśli
macierz C jest nieosobliwa) ma postać:
 1 le 1 le 
 2 −
8 2 8
 3 l 3 l 
− − e − e
C −1 = 4 8 4 8
(4.32)
 le le 
 0 − 0
 8 8 
 1 le 1 le 
−
 4 8 4 8 
Na podstawie zależności (4.31) otrzymuje się zatem:
[
Ne = 1 ξ ξ 2 ] [
ξ 3 C −1 = N1 N 2 N 3 N 4 ] (4.33)
gdzie:
1 3 1
N1 = − ξ + ξ3
2 4 4
le le l l
N2 = − ξ − e ξ 2 + e ξ 3
8 8 8 8
(4.34)
1 3 1
N3 = + ξ − ξ 3
2 4 4
le le l l
N4 = − − ξ + e ξ 2 + e ξ 3
8 8 8 8
Uzyskane funkcje kształtu N1 ÷ N 4 są funkcjami kształtu Hermite’a (rys. 4.3).

Odnosząc wyrażenie na całkowitą energię potencjalną belki (4.8) do
elementu „e” i uwzględniając przyjętą w obszarze belki aproksymację funkcji
ugięcia (4.28) otrzymuje się:
Π ec =
EJ
∫
[
 d 2 N e ( x) δ e ]  dx − ∫ q( x) N
2
e
( x) δ e dx =
2 le  dx 2  le
= ∫
[
EJ  d 2 N e ( x) e  ] 2
δ  dx − ∫ q( x) N e ( x) δ e dx = (4.35)
2 le  dx 2
 le
= ∫
[ ⋅
]
EJ 1  d 2 N e (ξ ) 4 e  l e
δ  ⋅ d ξ −
2
1
e
∫ q (ξ ) N (ξ ) δ
e le
dξ
2 −1  dξ 2 2
le  2 −1 2
W zależności powyższej wykorzystano następujące związki:
2
d 2 (•) d 2 (•)  dξ  d (•) d 2 ξ
= ⋅   + ⋅
dx 2 dξ 2  dx  dξ dx 2
(4.36)
∂x
dx = dξ
∂ξ
Na podstawie wyrażenia (4.35) i (3.77) otrzymuje się:
 dΠ ec 
 
 dw1 
 dΠ ec 
 
dΠ ec  dϕ x1 
=  e 
= k e δe + Fe (4.37)
dδ e  dΠ c 
 dw2 
 e 
 dΠ c 
 dϕ x 2 
Rys. 4.3. Funkcje kształtu elementu belkowego
Wyrazy k i j macierzy sztywności k e oraz składowe Fi wektora F e występu-

jących w zależności (4.37) przedstawiają się następująco:
8 1 d 2 N ie 
k ij = EJ 3 ∫  2
 dξ

l e −1 dξ 
(4.38)
l 1
Fi = e ∫ q (ξ ) N ie dξ
2 −1
Obliczając wyrazy k i j otrzymuje się znaną postać macierzy sztywności

elementu belkowego:
 12 6 12 6 
 l3 −
le 
2 3 2
 e le le
 4 6 2 
− 2
 le le le 
Ke =   (4.39)
 12 6
− 2
 le3 le 
 4 
 sym 
 le 
Na podstawie zależności (4.38) 3 , dla q (ξ ) = const. skierowanego w dół uzyskuje

się następującą postać wektora równoważnych sił węzłowych elementu:
T
e
 ql ql 2 ql ql e 2 
F = − e − e − e  (4.40)
 2 12 2 12 
4.2.1 Przykład – belka z przegubem

Do określenia jest rozkład momentów zginających i sił tnących w belce
obciążonej jak na rysunku 4.4.
Rys. 4.4. Belka przegubowa obciążona obciążeniem ciągłym
Do rozwiązania zadania zastosowano metodę elementów skończonych.

Przyjęty podział belki na elementy skończone przedstawiono na rys. 4.5, na
którym zaznaczono również stopnie swobody układu w globalnym układzie
odniesienia XOY.
Wektor przemieszczeń węzłowych układu w analizowanym zagadnieniu ma
następującą postać:
δ = [v1 ϕ1 v 2 ϕ 2 ϕ 2′ v3 ϕ 3 v 4 ϕ 4 ]T (4.41)
Składowe ϕ 2 i ϕ 2′ tego wektora są rotacyjnymi stopniami swobody węzła „2”.

Wynikają one z istnienia w tym węźle przegubu, co implikuje nierówność
ϕ 12 ≠ ϕ 12 .
Rys. 4.5. Podział belki na elementy skończone
Analiza na poziomie elementu
Macierze sztywności poszczególnych elementów przybiorą zgodnie z zależno-

ścią (4.39) formę:
 12 6 12 6 
 3 2
− 3 
 a a a a2 
 4 6 2 
−
K = K = K = EJ 
1 2 3 a a2 a 
(4.42)
12 6 
 − 2
 a3 a 
 sym 4 
 a 
Wektor równoważnych sił węzłowych ma dla elementu „1”, na którym jest

przyłożone obciążenie równomiernie rozłożone następujące składowe (por.
zależność 4.40)):
T
1  qa qa 2 qa qa 2 
F = − − −  (4.43)
 2 12 2 12 
Składowe wektorów równoważnych sił węzłowych dla pozostałych elementów,

a więc wektorów F 2 oraz F 3 są – w wyniku braku obciążenia na odpowiednich
elementach – równe zero.
Zgodnie z zależnością (4.37) dla poszczególnych elementów otrzymuje się:
 qa 
− 
 v 11   2 
   qa  2
−
ϕ 1  
1
dΠ 1C 
ELEMENT 1 = k 1   +  12  (4.44)
dδ 1 1 qa
v 2   − 
 1  2 
ϕ 2   qa 2 
 
 12 
 v12  0
   
dΠ C2  2 0
2 ϕ 1   
ELEMENT 2 =K  +  (4.45)
dδ 2 2
 v 2  0
 2  
ϕ 2  0
 v13  0
   
ELEMENT 3 dΠ C3  3 0
3 ϕ 1   
=K  +  (4.46)
dδ 3 3
 v 2  0
 3  
ϕ 2  0
W powyższych zależnościach wskaźniki górne występujące po lewej stronie

sumy odnoszą się do numeru elementu.
Analiza na poziomie układu
Równania równowagi układu, zgodnie z zależnością (4.78) mają postać:
K (9 x 9) ⋅ δ (9 x1) = F(9 x1) (4.47)

Macierz sztywności układu oraz wektor sił węzłowych układu otrzymane na

3 3
drodze agregacji K = ∑ k e , F = ∑ F e przedstawiają zależności (4.48) i (4.49).
e =1 e =1
(4.48)
T
 qa qa 2 qa qa 2 
F = − − − 0 0 0 0 0 (4.49)
 2 12 2 12 
W zależności (4.48) uwzględniono warunki zgodności przemieszczeń:
v11 = v1 v 22 = v13 = v3
ϕ 11 = ϕ 1 ϕ 22 = ϕ 13 = ϕ 3
(4.50)
v12 = v12 = v 2 v 23 = v 4
ϕ 12 = ϕ 2 ϕ 32 = ϕ 4
ϕ 12 = ϕ '2
W analizowanym zadaniu warunki brzegowe wynikające ze sposobu podparcia

są następujące:
v1 = v3 = v 4 = 0 (4.51)
Uwzględnienie powyższych warunków w zależności (4.49) i (4.50) prowadzi

w rezultacie do ostatecznej postaci równań równowagi układu:
(4.52)
Po przyjęciu danych: a = 1 m i q = 1 KN / m oraz po rozwiązaniu układu (4.52)

otrzymuje się następujące wartości przemieszczeń węzłowych:
 ϕ1   − 0,37498
   
 v 2   − 0,33333
ϕ  − 0,29168 1
    (4.53)
δ =  '2  =  ⋅
ϕ
 2  0, 41667  EJ
ϕ   0,16667 
 3  
ϕ 4   − 0,08333
Obliczenie sił wewnętrznych
Siły wewnętrzne oblicza się na poziomie elementu zgodnie z zależnością:
F ′ e = −F e + k e δ e (4.54)
Element „e1”:
 qa   12 6 12 6 
   3 2
− 3 
2 
 a a a a2  v
 T11   2 4 6 2  1  0,5
 1  qa   − ϕ   
 M 1   12   a a2 a  1  1 =  0 
 1  =  +
6  v   
EJ (4.55)
 12
 T2   qa   − 2   2  EJ 0,5
M 1   2   a3 a  ϕ 2   0 
 2 
qa 2   sym 4  
−   a 
 12 
Element „e2”:
 12 6 12 6 
 3 2
− 3 
 a a a a2 
 T12  0 4 6 2 v2  − 0,4996
 2    −  '  
M 1  0  a a2 a  ϕ 2  1 =  0 
=
 2   +
6   v  EJ  0,4996 
EJ  (4.56)
12
 T2  0  − 2  3   
M 2  0  a3 a  ϕ 3  − 0,4998
 2 4 
 sym
 a 
Element „e3”:
 12 6 12 6 
 3 2
− 3 
 a a a a2 
 T13  0 4 6 2   v3   0,5 
 3    − ϕ   
M 1  0  a a2 a  ⋅  3  ⋅ 1 =  0,5 
 3  =   + EJ  12 6      (4.57)
T
 2   0
 − 2   v 4  EJ − 0,5
M 3  0  a3 a  ϕ 4   0 
 2 4 
 sym
 a 
Obliczone uogólnione siły węzłowe przedstawiono dla poszczególnych

elementów na rys. 4.6, natomiast końcowy rozkład momentów zginających i sił
tnących w belce zobrazowano na rys. 4.7. Celem porównania otrzymanych
metodą elementów skończonych rezultatów obliczeń na rysunku tym
zamieszczono w nawiasach odpowiednie wartości rozwiązania dokładnego
analizowanej belki.
Rys. 4.6. Uogólnione siły węzłowe dla elementów
Rys. 4.7. Rozkład momentów zginających i sił tnących w belce

4.3. Izoparametryczny element płyty średniej grubości 113
4.3. Izoparametryczny element płyty średniej

grubości wyższego rzędu
Cechą elementów izoparametrycznych jest to, że te same funkcje
interpolacyjne użyte są do opisu zarówno geometrii elementu, jak i rozkładu
pola wewnątrz elementu. Liniowe elementy izoparametryczne (w których do
aproksymacji zastosowano funkcje liniowe) mają brzegi prostoliniowe. Brzegi
elementów z kwadratową zmiennością funkcji kształtu i elementów wyższych
rzedów są zakrzywione, co umożliwia wykorzystanie ich do modelowania
obszarów o brzegach krzywoliniowych. Ważną cechą elementów izopara-
metrycznych jest także to, że funkcje kształtu można wyznaczyć w sposób
bezpośredni, bez konieczności – jak w elementach subparametrycznych –
obliczania macierzy odwrotnej do macierzy współrzędnych, czyli macierzy C −1
(por. zależność 4.31)).
Wykorzystanie elementów izoparametrycznych w obliczeniach płyt
możliwe jest zasadniczo dla płyt średniej grubości, czyli płyt, w których
uwzględnia się wpływ sił poprzecznych i odkształceń postaciowych z nimi
związanych na deformację płyty. Dla płyt cienkich, czyli spełniających
założenia teorii Kirchhoffa zastosowanie elementów izoparametrycznych jest
praktycznie niemożliwe, co wynika z trudności związanych z transformacją
współrzędnych, dla której w przypadku płyt cienkich wymagana jest ciągłość
klasy C 1 .
Do obliczeń płyt średniej grubości metodą elementów skończonych
stworzono szereg elementów, głównie na bazie teorii Mindlina, rzadziej teorii
Reissnera. W chwili obecnej elementy skończone dyskretyzujące płytę Mindlina
tworzą bibliotekę szeregu profesjonalnych programów MES. Należy podkreślić,
że w zakresie rozkładu pola przemieszczenia wzdłuż grubości płyty zarówno
teoria płyt cienkich, Mindlina i Reissnera są teoriami o tym samym rzędzie
aproksymacji. W teoriach tych przyjmuje się bowiem, że wspomniane pole
przemieszczenia zmienia się w sposób liniowy. Uściślenie stanu płytowego
pojawia się w tzw. teoriach wyższego rzędu [13]. Przyjmuje się w nich, że
przemieszczenia zmieniają się wzdłuż grubości płyty nieliniowo, co jest zgodne
z ich rzeczywistym rozkładem, o ile działają siły tnące.
Jedną z teorii wyższego rzędu, interesującą z punktu widzenia zastosowań
praktycznych, jest teoria płyt średniej grubości Kujawskiego [8]. Poniżej
przedstawione będą związki podstawowe tego modelu płyty, a następnie
wyspecyfikowany będzie element skończony dyskretyzujący płytę średniej
grubości Kujawskiego.
4.3.1. Związki podstawowe modelu płyty wyższego rzędu

Rozpatruje się jednorodną, izotropową płytę o stałej grubości, poddaną działaniu
obciążenia powierzchniowego p( x, y ) (rys. 4.8).
Rys. 4.8. Deformacja przekroju poprzecznego płyty
W teorii Kujawskiego pole przemieszczenia określone jest następująco:
uα = z ϕ α0 + f p ϕ α1
(4.58)
0
u z = w = w( x, y ) α = x, y
gdzie uα , u z są jego składowymi odpowiednio w kierunku osi x, y i z ,

natomiast ϕ α0 , w 0 oraz ϕ 1α są niewiadomymi funkcjami zmiennych x, y ;
f p = z [1 − (5 z 2 ) / 3h 2 ] .
Pole przemieszczenia określone jest więc jako suma przemieszczeń zerowego
( zϕ α0 , w 0 ) i pierwszego ( f pϕ 1α ) rzędu. Można zauważyć, że przyjmując
w wyrażeniu (4.59) ϕ 1α = 0 otrzymuje się pole przemieszczenia, analogiczne do
założeń teorii Reissnera (zerowego rzędu).
Na podstawie zależności (4.58) oraz związków Cauchy’ego otrzymuje się
następujące wyrażenie na składowe pola odkształcenia:
1 1
ε αβ = ε 0 1
αβ + ε αβ = z (ϕ α0, β + ϕ 0β , α ) + f p (ϕ 1α , β + ϕ 1β , α )
2 2
1 1 (4.59)
ε α z = ε α0 z + ε 1α z = (ϕ α0 + w, α ) + f p1 ϕ 1α
2 2
εz z =0
15 z 2
gdzie f p1 = (1 − ) , β = x, y .
3 h2
Pole naprężenia przyjmuje natomiast postać związków:
0 1 E 1 0 0 ν 0
σ αβ = σ αβ+σ αβ= [ z (ϕ α , β + ϕ β ,α ) + δ α β zϕ γ ,γ +
1 +ν 2 1 −ν
1 ν ν (1 +ν )
+ f p (ϕ 1α , β + ϕ 1β , α ) + δ α β f p ϕ 1γ , γ + δα β σ z z ] (4.60)
2 1 −ν E (1 − ν )
E
σ α z = σ α0 z + σ 1α z = [(ϕ α0 + w, α ) + f p1 ϕ 1α
2(1 + ν )
Składową σ z z tensora naprężenia można otrzymać z układu różniczkowych

równań równowagi continuum trójwymiarowego. W teorii Kujawskiego
przyjmuje się analogiczną jak w teorii Reissnera postać wyrażenia określającego
naprężenia σ z z :
3p 2 z 1 z 3
σzz =− [ − + ( ) ] (4.61)
4 3 h 3 h
W zależnościach (4.59) i (4.60) przecinkiem oznaczono pochodne cząstkowe,

przez δ α β oznaczono deltę Kroneckera, natomiast γ = x, y .
Energia deformacji sprężystej płyty (por. zależność (3.67)) przyjmuje
wobec związków (4.59), (4.60) i (4.61) następującą postać wyrażoną
w przemieszczeniach:
D0 1 −ν 0
ΠK = ∫∫ { (ϕ α,β ϕ α0 , β + ϕ 0
β ,α ϕ 0
β ,α ) +ν δα β ϕ 0
γ, γ ϕ α0 , β +
2 Ωp 2
(4.62)
4 1 −ν
+ [ (ϕ 1α , β ϕ 1α , β + ϕ 1β , α ϕ 1β , α ) + ν δ α β ϕ 1γ , γ ϕ 1α , β ] +
21 2
3(1 − ν ) 2(1 − ν )
+ 2
(ϕ α0 + w,0α ) 2 − 2
(ϕ 0
α ϕ 1α + ϕ 1α w, α ) +
2h h
12(1 − ν ) 6 p ν (1 +ν ) 1 1
+ 2
ϕ 1α ϕ 1α + δ α β (ϕ 0
α, β + ϕ α , β )} dΩ p
3h 5 Eh 42
2 Eh 3
gdzie D 0 = , Ω p – obszar płyty.
3(1 − ν 2 )
Potencjał sił zewnętrznych w analizowanym zagadnieniu wyraża się
następująco (por. zależność (3.68)):
L = − ∫∫ p ( x, y ) w( x, y )dΩ p
Ωp
(4.63)
Całkowitą energię potencjalną rozpatrywanej płyty opisuje więc zależność:
ΠC = Π K + L (4.64)
gdzie Π K i L określone są odpowiednio przez związki (4.62) i (4.63).
4.3.2. Specyfikacja elementu płyty wyższego rzędu

Wyrażenie (4.64) stanowi podstawę do uzyskania równań równowagi
elementu skończonego. Równania te zostaną wyprowadzone na bazie
twierdzenia Lagrange’a oraz przy zastosowaniu podejścia Ritza. W tym celu
obszar płyty Ω p podzielony zostaje na elementy skończone, a analizie poddany
zostaje wyizolowany z układu element „e”.
Jak wynika z funkcjonału (4.63) energii deformacji sprężystej płyty
wyższego rzędu, w którym występują 1. pochodne funkcji ugięcia i kątów
obrotu, w modelu skończenie elementowym funkcje aproksymacyjne muszą być
z przestrzeni Vh ⊂ H 1 (Ω p ) . Oznacza to, że funkcje kształtu elementu powinny
być klasy C 0 (Ω p ) , a to umożliwia zastosowanie dla płyty wyższego rzędu
elementu izoparametrycznego.
Wobec powyższego jako element skończony przyjęto ośmiowęzłowy
element z serendipowskiej rodziny elementów zakrzywionych (rys. 4.9).
Jako stopnie swobody węzła przyjęto ugięcie powierzchni środkowej w
oraz cztery kąty obrotu: dwa zerowego rzędu ϕ 0x , ϕ 0y i dwa pierwszego rzędu
ϕ 1x , ϕ 1y . Tak więc wektor przemieszczeń w r -tym węźle elementu ma postać:
[
δ er = ϕ 0
xr ,ϕ 0
yr , wr , ϕ 1x r , ϕ 1y r ] T
r = 1,2,...,8 (4.65)
Wektor przemieszczeń wewnątrz elementu określony jest natomiast następująco:
ue = ϕ [ 0
x ϕ 0
y w ϕ 1x ϕ 1y ] T
(4.66)
Rys. 4.9. Przyjęty element skończony dla płyty wyższego rzędu,

(( ξ ,η ) – współrzędne lokalne, x , y , z – współrzędne globalne)
Geometrię elementu izoparametrycznego opisują zależności:
x = Nx
(4.67)
y = Ny
gdzie: N = [N1 N 2 L N 8 ] , natomiast składowymi wektorów x i y są odpo-

wiednie współrzędne globalne punktów węzłowych elementu.
Funkcje kształtu dla elementu o 8 węzłach opisane są zależnością:
N r = 0,25(1 + ξ 0 )(1 + η 0 )(ξ 0 + η 0 − 1) dla r = 1,2,3,4
N r = 0,5(1 + ξ 0 )(1 − η 2 ) dla r = 5,7

2
(4.68)
N r = 0,5(1 − ξ )(1 + η 0 ) dla r = 6,8
gdzie ξ 0 = ξ ξ r ,η 0 = η η r są nowymi zmiennymi, za pomocą których

wszystkie funkcje kształtu można opisać w jednolity sposób.
Zgodnie z podejściem Ritza, związek miedzy przemieszczeniami
wewnątrz elementu a przemieszczeniami węzłowymi ma postać:
ϕ 0x 
 
ϕ 0y 
 
u e =  w  = [N 1 , N 2 , K , N 8 ] δ e (4.69)
 1
ϕ x 
 1
ϕ y 
przy czym:
N r = Nr I5 (4.70)
gdzie I 5 jest macierzą jednostkową o wymiarze (5 x5) .

Odnosząc wyrażenie na całkowitą energię potencjalną płyty (4.64) do elementu,
przy uwzględnieniu związków (4.69) oraz zależności
 N i , x  N 
−1  i , ξ 
 =J  
 N i, y   N i ,η 
(4.71)
1 1
∫∫ (•) ( x, y ) dx dy = ∫ ∫ (•)[ x(ξ ) y (η ) ] det J dξ dη
Ωe −1−1
w której J jest macierzą Jakobiego:
N ,ξ 
J=  [x y ] (4.72)
 N ,η 
uzyskuje się:
 ∂Π C 
 0 
 ∂ϕ x 
 ∂Π C 
 0 
 ∂ϕ y 
∂Π Ce  ∂Π C 
= = k eK δ e + F e
∂w 
(4.73)
∂δ e  
 ∂Π C 
 1
 ∂ϕ x 
 ∂Π C 
 1
 ∂ϕ y 
gdzie macierz sztywności elementu płyty wyższego rzędu k eK można

przedstawić w następującej formie:
k 11
K k 12
K L k 18
K

 22 28 
kK L kK 
k eK = (4.74)
 O M 
  ( 40 x 40)
 S Y M k 88
K 
Przykładowe elementy podmacierzy k rKm (5 x 5) określają związki:
1 −ν 3(1 − ν )
k rKm11 = D 0 ∫∫ ( N r , x N m, x + N r , y N m, y + 2
N r N m ) dΩ e
Ωe
4 4h
1 −ν (4.75)
k rKm12 = D 0 ∫∫ (ν N r , x N m, y + N r , y N m , x ) dΩ e
Ωe
4
r , m = 1,2,...,8
Wektor F e (zależność 4.73) przedstawia się następująco:
F e = [F1 F2 K F8 ] T (4.76)
gdzie:
Fr = [ 0 0 Fr 0 0 ]T (4.77)
oraz:
Fr = ∫∫ p ( x, y ) N r dΩ p
Ωp
(4.78)
Należy nadmienić, że cechą charakterystyczną sformułowanych dla płyt

średniej grubości elementów skończonych, w których do aproksymacji pola
przemieszczenia zastosowano funkcje spełniające ciągłość klasy C 0 , jest
nadmierna sztywność w przypadku niewielkiej grubości płyty. Efekt ten znany
w literaturze pod nazwą „lockingu” związany jest z koniecznością spełnienia
przy małej grubości płyty warunków ε x z = ε y z = 0 , odpowiadających założe-
niom płyt cienkich Kirchhoffa. W przypadku natomiast wymienionych
elementów wraz ze zmniejszaniem się ich grubości, występujące w macierzy
sztywności wyrazy odpowiadające deformacji stycznej płyty stają się coraz
bardziej dominujące, co powoduje, że otrzymane wyniki numeryczne są błędne.
W celu ograniczenia tego niekorzystnego zjawiska wskazuje się na konieczność
stosowania technik ogólnego lub selektywnego zredukowanego całkowania.
Warto podkreślić, że przedstawiony element płyty wyższego rzędu nie wymaga
stosowania wyżej wspomnianych technik, wykazując przy pełnym całkowaniu
(3 x3 punktów Gaussa) i małej grubości płyty dobrą zgodność z odpowiednimi
rezultatami uzyskanymi dla płyty Kirchhoffa.
4.4. Elementy nieskończone

W wielu problemach naukowych i inżynierskich pojawia się potrzeba
analizowania pola w obszarach nieograniczonych. Zastosowanie MES w jej
klasycznej postaci do tego typu zagadnień sprowadza się do wyznaczenia
dyskretnych wartości pola w odpowiednio dużym, jednak skończonym obszarze,
przyjętym przy założeniu, że poszukiwane wartości pola na jego brzegach
praktycznie zanikają. Prowadzi to do rozbudowanej siatki elementów i dużej
liczby niewiadomych układu równań MES. Z drugiej strony przyjęcie a priori
skończonych rozmiarów analizowanego obszaru w przypadku jego nieogra-
niczoności może znacząco wpływać na otrzymane wyniki obliczeń. W celu
rozszerzenia możliwości zastosowania MES do analizy pola w obszarach
nieograniczonych – bez konieczności dyskretyzowania dużego obszaru
i wystąpienia związanych z tym niedogodności – sformułowano elementy
nieskończone. Elementy te, chociaż utworzone stosunkowo niedawno, zostały
już umieszczone i z powodzeniem wykorzystywane w niektórych profe-
sjonalnych programach MES (na przykład ABAQUS 2002).
4.4. Elementy nieskończone 121
4.4.1 Zasadnicza koncepcja elementów nieskończonych

W wyznaczaniu pola w obszarach nieograniczonych szczególnie
przydatne – z uwagi na łatwość implementacji numerycznej i uzyskiwane
wyniki – są odwzorowane elementy nieskończone („mapped infinite elements”)
[9]. Idea tych elementów polega na tym, że geometria elementu jest określona
przez parametryczne odwzorowanie elementu rodzimego przy użyciu funkcji
kształtu wykazujących osobliwość dla ξ lub η = 1 w przypadku elementu
pojedynczo nieskończonego lub dla ξ i η = 1 w przypadku elementu podwójnie
nieskończonego. Przemieszczenia natomiast w obszarze elementu są aproksy-
mowane przez standardowe funkcje kształtu dla elementu izoparametrycznego.
Ideę odwzorowanych elementów nieskończonych można zobrazować,
prowadząc rozważania na przykładzie jednowymiarowego elementu nieskoń-
czonego (rys. 4.10). Element ten można odwzorować na element opisany w
lokalnym układzie współrzędnych − 1 ≤ ξ ≤ 1 . Zależność między współrzę-
dnymi globalnymi i lokalnymi przedstawia się następująco:
x 
x(ξ ) = [M 1 M2]  1 (4.79)
x2 
Rys. 4.10. Jednowymiarowy element nieskończony
Występujące w wyrażeniu (4.79) funkcje odwzorowujące M 1 i M 2 opisane są

następująco:
− 2ξ
M1 =
1−ξ
(4.80)
1+ ξ
M2 =
1−ξ
natomiast x1 , x 2 są odpowiednimi współrzędnymi globalnymi punktów węzło-

wych 1, 2 elementu.
W elementach nieskończonych geometria oraz rozkład pola prze-
mieszczenia są odniesione do tego samego punktu lub zbioru punktów zwanych
biegunami elementu. Na rys. 4.10 oznaczono je przez O. Bieguny powinny
znajdować się na zewnątrz elementu, a więc x0 < x1 . Zachodzą ponadto
zależności:
x1 = x 0 + a 0
(4.81)
x2 = x0 + a0 + a
Zakłada się przykładowo, że poszukiwaną funkcją jest funkcja φ ,

a wektor przemieszczeń węzłowych elementu ma postać:
δ e = [δ 1 δ 2 δ 3 ] T = [φ 1 φ 2 φ 3 ]T (4.82)
Zmienność funkcji φ w obszarze elementu można przedstawić za pomocą

standardowego wyrażenia :
φ = [N 1 N2 N3 ] δe (4.83)
gdzie
N 1 = −0,5(1 − ξ )
N 2 = (1 − ξ )(1 + ξ ) (4.84)
N 3 = 0,5(1 + ξ )
Na podstawie zależności (4.79) i (4.81) otrzymuje się:
− r + a0 + a
ξ= (4.85)
− r + a0 − a
gdzie r oznacza odległość dowolnego punktu wewnątrz elementu od bieguna O.

Przy wykorzystaniu związku (4.85) zależność (4.83) przyjmie postać:
a + a0 a
−1 1− 0
a 4 a
φ= ⋅ r ⋅φ 1 + ⋅ r ⋅φ
2 +
a0 − a + r a 0 − a a 0 − a − r a0 − a
−1 −1
r r
(4.86)
a0 a + a0
−1 −1
+ r ⋅ r ⋅φ 3
a0 − a a −a
−1 0 −1
r r
Przyjmując w powyższym wyrażeniu a 0 = a , otrzymuje się:
a a2
φ = φ 3 + (−φ 1 + 4φ 2 − 3φ 3 ) + (2φ 1 − 4φ 2 + 2φ 3 ) 2 (4.87)
r r
Dla zależności (4.86) i (4.87) zachodzi:
lim φ = φ 3 (4.88)
r →∞
Tym samym więc warunek, że funkcja φ osiąga w nieskończoności wartość

zero można automatycznie spełnić przez rozpatrzenie w wyrażeniu (4.83) tylko
„skończonych” węzłów 1 i 2. Na podstawie związków (4.86) i (4.87) widoczne
jest również znaczenie umiejscowienia bieguna O. Należy podkreślić, że dobór
położenia bieguna może być zrealizowany jedynie na drodze empirycznej, przy
uwzględnieniu tak geometrycznych jak i fizycznych własności danego zadania.
4.4.2. Dwuwymiarowe elementy pojedynczo i podwójnie

nieskończone
Przedstawiona w podrozdziale 4.4.1 zasadnicza koncepcja odwzorowanych
elementów nieskończonych znajduje odzwierciedlenie w przypadku elementów
dwuwymiarowych, pojedynczo lub podwójnie nieskończonych. Przykłady
takich elementów przedstawiono na rys. 4.11. Dla elementów tych przyjęto
założenie, że a 0 = a . Element pojedynczo nieskończony przedstawiony na
rys. 4.11 a) jest pięciowęzłowym, serendipowskim elementem, pochodzącym od
ośmiowęzłowego serendipowskiego elementu skończonego.
Współrzędne globalne punktów 5 i 4 tego elementu określone są zależnością:
y 4 = 2 y3 − y 0
(4.89)
y 5 = 2 y1 − y 0
Rys. 4.11. Dwuwymiarowe elementy nieskończone:

a) element pojedynczo nieskończony, b) element podwójnie nieskończony
Geometrię elementu opisuje wyrażenie:
x=Mx
(4.90)
y=My
gdzie:
M = [M 1 M2 L M5]
x = [x1 x2 L x5 ] T (4.91)
y = [ y1 y 2 L y5 ] T
Przyjęte funkcje odwzorowujące mają postać:

(1 − ξ )(−1 − ξ − η )
M1 =
1 −η
2(1 − ξ 2 )
M2 =
1 −η
(1 + ξ )(−1 + ξ − η )
M3 = (4.92)
1 −η
(1 + ξ )(1 + η )
M4 =
2(1 − η )
(1 − ξ )(1 + η )
M5 =
2(1 − η )
Można zauważyć, że dla funkcji (4.92) zachodzi zależność:
η = 1 ⇒ y (ξ ,η ) → ∞ (4.93)
Element podwójnie nieskończony przedstawiony na rys. 4.11 b) jest z

kolei czterowęzłowym, izoparametrycznym elementem Lagrange’a, pocho-
dzącym od elementu skończonego Lagrange’a o dziewięciu węzłach. Analo-
gicznie jak w przypadku elementu pojedynczo nieskończonego, współrzędne
globalne punktów 4 i 2 określone są następująco:
x 2 = 2 x1 − x 0
(4.94)
y 4 = 2 y1 − y 0
Geometrię elementu przedstawia zależność (4.90), w której:
M = [M 1 M2 K M4]T
x = [x1 x2 x3 x4 ] T (4.95)
y = [ y1 y2 y3 y4 ] T
Składowe macierzy M mają obecnie postać:

4ξη
M1 =
(1 − ξ )(1 − η )
(1 + ξ )(−2η )
M2 =
(1 − ξ )(1 − η )
(4.96)
(1 + ξ )(1 + η )
M3 =
(1 − ξ )(1 − η )
− 2ξ (1 + η )
M4 =
(1 − ξ )(1 − η )
Dla przyjętych funkcji odwzorowujących (4.96) zachodzi:
ξ = 1 ∧ η = 1 ⇒ x(ξ ,η ) → ∞ ∧ y (ξ ,η ) → ∞ (4.97)
4.5. Makroelement podłoża warstwowego

Współpraca konstrukcji z dowolnie uwarstwionym podłożem gruntowym
jest problemem złożonym, nawet jeżeli przyjmie się założenie, że grunty
poszczególnych warstw są materiałami izotropowymi i liniowo-sprężystymi.
Dość powszechnie stosowane do opisu ośrodka gruntowego w zagadnieniach
interakcji konstrukcja-podłoże modele analityczne (Winklera, Pasternaka,
Reissnera) okazują się tu być daleko niewystarczające. Modelują one ośrodek
gruntowy często w oderwaniu od jego rzeczywistych własności mechanicznych,
nie umożliwiają również uwzględnienia tej zasadniczej cechy geometrycznej
podłoża, jaką jest niejednorodność strefowa w kierunku pionowym. Przy
zróżnicowanych modułach ściśliwości poszczególnych warstw może to skut-
kować błędnymi wynikami naprężeń i odkształceń. Realistyczne oszacowania
tak osiadań jak i deformacji konstrukcji spoczywającej na geometrycznie
skomplikowanym układzie warstw wymagają pełnej, trójwymiarowej analizy
MES, która jest z reguły zadaniem czasochłonnym i kosztownym.
Alternatywną grupę w zakresie modelowania ośrodka gruntowego
w zagadnieniach interakcji, pozwalającą na osiągnięcie pewnego kompromisu
między dokładnością a czasochłonnością obliczeń stanowią modele numeryczne
podłoża z efektywnym modelem Kolara-Nemeca (K-N) [6] na czele. Symulacja
zachowania podłoża przez ten model pozwala na znaczne uproszczenie procesu
obliczeniowego, przy wyrażeniu wszystkich niezbędnych własności podłoża,
które w analizie układu konstrukcja-podłoże warstwowe nie mogą być
pominięte.
4.5. Makroelement podłoża warstwowego 127
4.5.1. Związki podstawowe modelu podłoża warstwowego

Punktem wyjścia do rozważań jest układ n stref obliczeniowych o zmiennych
grubościach h (i ) ( x, y ) , (i = 1,2,K , n) , który spoczywa na podłożu nieod-
kształcalnym (rys. 4.12).
Rys. 4.12. Podłoże dowolnie uwarstwione
Grunty poszczególnych warstw traktowane są jako materiały liniowo-sprężyste

i niejednorodne, wyspecyfikowane w ogólności za pomocą modułów Younga
i Kirchhoffa, odpowiednio E (i ) = E (i ) ( z ) oraz G (i ) = G (i ) ( z ) . Wszystkie
granice warstw geologicznych muszą być pokryte powierzchniami z = H (i )
stref obliczeniowych, lecz nie odwrotnie.
Istotą modelu podłoża warstwowego, który oparty jest na tzw. efek-
tywnym modelu Kolara-Nemeca są ograniczenia nałożone na przemieszczenia:
uα(i ) ( x, y , z ) ≅ 0 (α = x, y )
u z(i ) ( x, y , z ) = w (i ) ( x, y , z ) = w(i ) ( x, y ) ⋅ f (i ) ( x, y, z ) + (4.98)
+ w(i +1) ( x, y ) ⋅ f (i +1) ( x, y, z )
W zależności (4.98) wskaźniki ujęte w nawiasy nie podlegają konwencji

sumacyjnej. W zależności tej poszczególne symbole oznaczają:
uα(i ) , u z(i ) – przemieszczenia punktów warstwy "i" w kierunkach osi współ-

rzędnych x, y, z ; w(i ) , w(i +1) – ugięcia odpowiednich płaszczyzn granicznych
z = H (i −1) , przy czym H ( 0 ) = 0 ; f (i ) , f (i +1) – dane funkcje interpolacyjne
definiujące pionowy rozkład przemieszczenia w (i ) wewnątrz warstwy "i" .
Funkcje f (i ) oraz f (i +1) powinny spełniać graniczne warunki geometryczne:
0 z = H (i −1) 1 z = H (i )
 
f (i ) = 1 z = H (i ) f (i +1) = 0 z = H (i −1) (4.99)
 
 
W modelu K-N funkcje te opisane są następująco (rys. 4.13):
H (i ) − z z − H (i −1)
f (i ) = , f (i +1) = i = 1,2,..., n − 1
H (i ) − H (i −1) H (i ) − H (i −1)
(4.100)
H (i ) − z
f (i ) = i=n
H (i ) − H (i −1)
Moduł sprężystości E (i ) ( z ) w obrębie warstwy uzależniony jest od modułu

przypisanego powierzchni środkowej i -tej warstwy E (i ) w sposób następujący:
E (i ) ( z ) = E (i ) [1 + kξ (i ) ] , −1≤ k ≤1 (4.101)
z
gdzie: ξ (i ) = (i )
[ z − H (i −1) ] − 1 jest współrzędną lokalną i -tej warstwy
h
(rys. 4.14).
Rys. 4.13. Wykresy funkcji opisujących rozkład ugięć wzdłuż grubości podłoża
Rys. 4.14. Globalny i lokalny układ współrzędnych i -tej warstwy
Przyjmując, że warstwy mogą mieć zmienną grubość, zakłada się, że

powierzchnie graniczne są płaszczyznami dowolnie nachylonymi w stosunku do
płaszczyzny xoy (rys.4.12). Można je więc opisać równaniami:
H (i ) ( x, y ) = a (i ) + b (i ) x + c (i ) y + d (i ) xy (4.102)
gdzie a, b, c, d są parametrami, które wyznacza się na podstawie współrzędnych

punktów leżących na odpowiedniej powierzchni granicznej.
Na podstawie zależności (4.98) oraz związków Cauchy’ego otrzymuje się
następujące wyrażenia opisujące składowe pola odkształcenia:
(i )
ε αβ =0
(i ) 1
ε α z= ( w(i ) f (i ) + w(i +1) f (i +1) ) ,α (4.103)
2
ε (ziz) = ( w(i ) f (i ) + w(i +1) f (i +1) ) , z
Związki fizyczne opisujące przyjęty ośrodek można w przypadkach anizotropii

zawartych w przedziale <izotropia, ortotropia> zapisać w postaci:
(i ) (i ) (i )
σ ij=G ij ε i j + δi j c (jik) ε (i )
kk i, j = x, y , z (4.104)
przy czym δ i j jest deltą Kroneckera.

(i )
Występujące w równaniach (4.106) parametry materiałowe G ij oraz c (ij k) i -tej
warstwy tworzą macierze:
G (xi x) G (xi )y G (xi z)   0 K c (xi z) 

   
G (i ) = M , c (i ) =  M  (4.105)
 sym G (ziz)   sym 0 
 
W przypadku ośrodka izotropowego parametry te wyrażają się przez techniczne

stałe wzorami:
E (1 − ν )
G(i i ) =
(1 + ν )(1 − 2ν )
E
G(i j ) = i≠ j (4.106)
2(1 + ν )
Eν
G(i j ) = i≠ j
(1 + ν )(2 − ν )
Na podstawie związków (4.103) i (4.104) otrzymuje się pole naprężenia

o następujących składowych:
(i )
σ αβ = δ αβ c β(i )z ( w(i ) f (i ) + w(i +1) f (i +1) ) , z
(i ) 1 (i )
σ α z= G α z ( w(i ) f (i ) + w(i +1) f (i +1) ) ,α (4.107)
2
σ (ziz) = G α(i )3 ( w(i ) f (i ) + w(i +1) f (i +1) ) , z
Wyrażenie na energię sprężystą (por. zależność (3.67)) przyjmuje w rozpa-

trywanym zagadnieniu następującą postać:
1 n (i ) (i ) ( i ) (i ) (i )
ΠW = ∑ ∫∫∫ (σ z z ε z z + 2σ α z ε α z ) dV , (4.108)
2 (i ) =1 V (i )
która wynika z przyjętego założenia (4.98) 1 .

Podstawiając do zależności (4.108) związki (4.103) i (4.107) otrzymuje się:
1 n (i )
ΠW = ∑ ∫∫ {G z z [( w( i ) f (i ) + w(i +1) f (i +1) ) , z ] 2 +
2 (i ) =1 h(i )
(4.109)
1
+ G α(i )z [(w(i ) f (i ) + w(i +1) f (i +1) ) ,α ] 2 } dz (i ) dΩ
2
gdzie przez ∫ oznaczono umownie całkowanie w granicach grubości odpowie-

h(i )
dniej warstwy.
Całkowita energia potencjalna przy uwzględnieniu pracy sił zewnętrznych
ma postać:
Π C = ΠW + L (4.110)
gdzie Π W określa zależność (4.109), natomiast potencjał sił zewnętrznych

opisany jest następująco:
n
(i ) (i ) (i )
L= ∑ ∫ q ( xα ) ⋅ w ( x, y, z ) dV (4.111)
( i ) =1 V ( i )
4.5.2. Specyfikacja makroelementu podłoża

Związek (4.110) stanowi podstawę do uzyskania równań MES. W tym
celu ciągły model podłoża zostaje podzielony wzajemnie prostopadłymi płasz-
czyznami pionowymi na elementy skończone (rys. 4.15) [10].
Makroelement podłoża jest zbiorem płaskich elementów skończonych
pokrywających powierzchnie graniczne warstw (rys. 4.16).
Funkcjonał (4.111) można przedstawić w postaci sumy energii poten-
cjalnych poszczególnych warstw:
1 (1) ( 2) (n)
ΠW = (Π W +ΠW +K+ Π W ) (4.112)
2
gdzie:
(i ) (i )
ΠW = ∫∫ ∫ {G zz [( w(i ) f (i ) + w(i +1) f (i +1) ) , z ] 2 +
Ω h(i )
(4.113)
1 (i )
+ G α z [(w(i ) f (i ) + w(i +1) f (i +1) ) ,α ] 2 } dz (i ) dΩ
2
Rys. 4.15. Dyskretyzacja MES podłoża warstwowego
Rys. 4.16. Makroelement podłoża warstwowego
Analogicznie można przedstawić potencjał sił zewnętrznych (4.111):
L = L(1) + L( 2 ) + K + L( n ) (4.114)
gdzie:
L(i ) = ∫∫ ∫ q (i ) ( x, y ) [ w(i ) ( x, y ) f (i ) ( x, y, z ) +
Ω h(i )
(4.115)
+ w(i +1) ( x, y ) f (i +1) ( x, y, z )] dz (i ) dΩ
(i )
Widoczne jest, że w funkcjonale Π W i -tej warstwy występują 1. pochodne
funkcji ugięcia, a więc w modelu skończenie elementowym funkcje aproksy-
macyjne powinny być z przestrzeni Vh ⊂ H 1 (Ω).
Dla i -tej warstwy (rys. 4.17) przyjmuje się liniowy rozkład ugięć wzdłuż
grubości warstwy, opisany funkcjami kształtu:
N (i ) = 0,5(1 − ξ (i ) )
(4.116)
N (i +1) = 0,5(1 + ξ (i ) )
Rys. 4.17. Liniowy rozkład przemieszczeń w obrębie i -tej warstwy
Tak więc funkcje interpolacyjne f (i ) i f (i +1) określone są następująco:
f (i ) = N (i )
(4.117)
f (i +1) = N (i +1)
Ugięcie i -tej warstwy można zatem zapisać zgodnie z zależnością (por.

zależność (4.98)):
w (i ) ( x, y , z ) = w(i ) ( x, y ) N (i ) + w(i +1) ( x, y ) N (i +1) =

 w(i )  (4.118)
[
= N (i ) ]
N (i +1)  
w(i +1) 
W obrębie powierzchni granicznych i -tej warstwy przyjmuje się czterowęzłowe

elementy prostokątne o SSW=1 (rys. 4.18).
Rys. 4.18. Elementy skończone usytuowane na powierzchniach granicznych

i -tej warstwy
Wektor przemieszczeń węzłowych dla i -tej warstwy ma zatem postać:
[
δ (i ) = δ (i ) δ (i +1) ]T (4.119)
[
δ ( i ) = w(i )1 w(i ) 2 K w( i ) r ] T
gdzie: (4.120)
[
δ ( i +1) = w( i +1)1 w(i +1) 2 K w(i +1) r ] T
, r=4
Wektor przemieszczeń węzłowych makroelementu (rys. 4.16) przedstawia się

następująco:
[
δ M = δ (1) δ ( 2) K δ ( n−1) ] = [δ
T
(1) δ ( 2) K δ (n) ]T (4.121)
Zgodnie z podejściem Ritza związek między przemieszczeniami wewnątrz

elementu a przemieszczeniami węzłowymi opisuje zależność:
 w(i )1 
w 
 (i ) 2 
w(i ) = [N 1 N2 N3 N4 ]   = N δ (i )
 w(i ) 3 
 w( i ) 4 
 
 w(i +1)1 
w 
 (i +1) 2 
w(i +1) = [N 1 N2 N3 N 4 ]  = N δ (i +1) (4.122)
 w(i +1) 3 
w(i +1) 4 
 
gdzie:
N 1 = 0,25(1 − ξ )(1 − η )
N 2 = 0,25(1 + ξ )(1 − η )
(4.123)
N 3 = 0,25(1 + ξ )(1 + η )
N 4 = 0,25(1 − ξ )(1 + η )
Podstawiając związki (4.117) i (4.122) do wyrażeń (4.113) oraz (4.115)

odniesionych do elementu „e” otrzymuje się:
Π (wi ) e = ∫∫ ∫ {G (i )
z z [(N δ ( i ) N (i ) + N δ (i +1) N (i +1) ) , z ] 2 +
Ω e h(i )
(4.124)
1 (i )
+ G α z [(N δ (i ) N (i ) + N δ (i +1) N (i +1) ) ,α ] 2 } dz (i ) dΩ e
2
oraz:
e
L(i ) = ∫∫ (i ) (i )
∫ q [N δ (i ) N (i ) + N δ (i +1) N (i +1) ] dz dΩ
e
(4.125)
Ω e h(i )
Biorąc pod uwagę zależność:
z
ξ (i ) = (i )
[ z − H (i −1) ] − 1 (4.126)
h
otrzymuje się:
∂z h (i )
dz = ⋅ dξ (i ) = dξ ( i ) (4.127)
∂ξ (i ) 2
Zatem:
(i ) e (i )
1 h (i )
ΠW = ∫∫ e
∫ (•) dz dΩ = ∫∫ ∫ (•) dΩ
e
dξ (i ) (4.128)
Ω e h(i ) Ω e −1
2
oraz:
e 1 h (i )
L(i ) = ∫∫ (i ) e
∫ (•) dz dΩ = ∫∫ ∫ (•) dΩ
e
dξ ( i ) (4.129)
Ω e h(i ) Ω e −1
2
Dokonując ścisłego całkowania analitycznego zależności (4.128) i (4.129)

względem zmiennej ξ (i ) , co jest zadaniem nieskomplikowanym wobec postaci
(i ) e
funkcji kształtu (4.116), otrzymuje się funkcjonał Π w i potencjał sił
e
zewnętrznych L (i ) zależne jedynie od współrzędnych x, y . Tak więc:
(i ) e (i ) e
ΠW ( x, y , z ) → Π W ( x, y )
e e (4.130)
L ( i ) ( x , y , z ) → L ( i ) ( x, y )
Funkcjonał energii potencjalnej makroelementu przedstawia się zgodnie

z zależnością (4.114) następująco:
(1) e ( 2) e ( n) e
Π WM = Π W +ΠW +K+ Π W (4.131)
Podobnie można przedstawić potencjał sił zewnętrznych makroelementu (por.

zależność (4.114)):
e e ( n)e
LM = L (1) + L ( 2 ) + K + L (4.132)
Na podstawie zależności (3.77) uzyskuje się:
 ∂Π CM 
 
 δ (1) 
 M 
∂Π CM ∂ (Π WM + LM )  ∂Π C  M M M
= = δ  = kW δ + F (4.133)
∂δ e ∂δ e  M
( 2)

 M 
 ∂Π C 
 δ (n) 
 
gdzie macierz sztywności makroelementu k WM można przedstawić w nastę-

pującej formie:
 k 11
W k 12
W L k 1Wn 
 22 2n 
kW L kW
k WM = 
(4.134)
 O M 
 nn

 sym kW 
ik
Podmacierze k W (i, k = 1,2, K , n) macierzy k WM opisuje wyrażenie (w którym
(i )
nie uwzględniono całkowania analitycznego względem zmiennej ξ ):
ik
n
(i ) 1 (i )
kW r m = ∫∫ ∑ ∫ {G zz N r N m N (i ), z N ( k ), z + G α z ( N m, α N ( k ) +
Ω e ( i ) =1 h ( i )
2
(4.135)
+ N m N ( k ), α )( N r , α N (i ) + N m N (i ),α )} dz (i ) dΩ e
r , m = 1,2, K ,4
Wektor F M przedstawia się natomiast następująco:
F M = F(1) [ F( 2 ) K F( n ) ]T (4.136)
gdzie:
[
F(i ) = F(i )1 F(i ) 2 K F(i ) 4 ] T
(4.137)
oraz:
F(i ) r = ∫∫ ∫ q (i ) N (i ) N r dz (i ) dΩ e r = 1,2, K ,4 (4.138)
Ωe h (i )
Warto nadmienić, że w odróżnieniu od przyjętego w modelu K-N rozkładu

modułu sprężystości wzdłuż grubości warstwy (por. zależność (4.103)) można
przyjąć opis bardziej zbliżony do rzeczywistości, zgodnie z poniższym
wyrażeniem:
(i )
E (i ) ( z ) = E 0 + E 1(i ) ( z ) (4.139)
przy czym współczynniki E 0 i E1 wyznacza się dla odpowiedniej warstwy na

podstawie liniowej regresji wartości modułów uzyskanych z pomiarów na
różnych głębokościach.
4.5.3. Uwagi dotyczące kalibracji parametrów modelu

podłoża warstwowego
Praktyczny sens stosowania makroelementu podłoża warstwowego polega
na redukcji analizy MES zagadnienia trójwymiarowego do dwóch wymiarów
przy wykorzystaniu znanych funkcji rozkładu przemieszczeń pionowych wzdłuż
grubości warstw. W racjonalnym makroelemencie warstwowym uwzględnia się
też znaną (na przykład [4]) zmienność modułów Younga E (i ) i Kirchhoffa G (i )
na głębokości. Innym aspektem tej operacji jest znaczne przyśpieszenie procesu
obliczeniowego. Zamiast trójwymiarowej siatki elementów skończonych
podłoża o trzech stopniach swobody w każdym węźle stosuje się zbiór płaskich
elementów na granicach warstwy o jednym stopniu swobody.
Wymienione cechy modelu wskazują na jego dużą przydatność w
liniowych analizach zagadnień współpracy konstrukcji z podłożem gruntowym.
Względnie miarodajna ocena użyteczności modelu wymaga jednak odniesienia
się do problemu kalibracji jego parametrów. W modelu podłoża warstwowego
parametrami tymi są moduł odkształcalności i moduł Kirchhoffa warstwy: E (i )
i G (i ) . Możliwość wiarygodnej identyfikacji tych parametrów stwarzają
współczesne badania penetracyjne gruntu. Wymienić tu można sondowanie
statyczne CPTU (zalecane przez system europejskich norm geotechnicznych
ujętych w EUROKODZIE 7), badania dylatometrem (testy DMT) względnie
badania presjometryczne [4]. Przy wykorzystaniu sondy CPTU wartości
parametrów geotechnicznych mierzone są w sposób ciągły na całym przelocie
otworu badawczego. Uzyskane z sondowań dane pozwalają na precyzyjne
określenie zarówno modułu odkształcalności na całej badanej głębokości, jak
również granic warstw geotechnicznych o zbliżonych parametrach wytrzyma-
łościowych. Eliminuje to uciążliwe poszukiwanie zależności teoretycznych
zmian modułu odkształcalności, które z reguły nie odzwierciedlają jego
rzeczywistego rozkładu, jako że pomija się w nich wpływ prekonsolidacji
(OCR). Szczegółowy opis sposobów pozyskiwania realistycznych oszacowań
modułu odkształcalności podłoża znaleźć można w pracy [4].
Rozdział 5
PRZYKŁADY ZASTOSOWAŃ MRS I MES
W ANALIZIE STATYCZNEJ
KONTRUKCJI POSADOWIONYCH
NA PODŁOŻU SPRĘŻYSTYM
5.1. Uwagi wstępne

W niniejszym rozdziale przedstawione będą przykłady zastosowań metody
różnic skończonych i metody elementów skończonych w rozwiązywaniu
zarówno prostych zagadnień współpracy konstrukcji z podłożem sprężystym,
jakimi są belki na podłożu Winklera, jak i zagadnień złożonych, do których
niewątpliwie zaliczyć można płyty wyższego rzędu na podłożu warstwowym.
Przyjęte do opisu zachowania się ośrodka gruntowego pod obciążeniem modele
podłoża sprężystego, a mianowicie model Winklera, warstwa Własowa i model
K-N są reprezentantami trzech zasadniczych grup modeli podłoża stosowanych
w zagadnieniach interakcji konstrukcja – podłoże, tj. modeli prostych, podsta-
wowych i złożonych. Ilustruje to schemat klasyfikacyjny modeli podłoża
przedstawiony na rys. 5.1 [3].
Pierwsza grupa modeli podłoża – analogi globalne, której przedsta-
wicielem jest model Winklera obejmuje schematy podłoża, które nie repre-
zentują masywu gruntowego w ścisłym sensie, a jedynie są zastępczymi
układami mechanicznymi, konstruowanymi w ten sposób, aby rozkład ugięć
powierzchni dość dobrze naśladował krzywe doświadczalne. Model Winklera,
który jest zestawiony z wzajemnie niepowiązanych sprężyn o jednakowej
140 Przykłady zastosowań MRS i MES w analizie statycznej
charakterystyce sprężystej, cechuje lokalność i proporcjonalność związku

„obciążenie – osiadanie”, co oznacza daleko idącą idealizację rzeczywistego
masywu. W rezultacie model ten nie daje pewnych oszacowań, a jego dobra
zgodność z rezultatami badań doświadczalnych i bardziej precyzyjnych analiz
numerycznych jest często przypadkowa.
Rys. 4.1. Schemat klasyfikacyjny modeli podłoża sprężystego stosowanych

w zagadnieniach interakcji konstrukcja - podłoże
Warstwa Własowa należy do grupy modeli podstawowych, w skład której

wchodzą także modele półprzestrzeni sprężystej i warstwy klasycznej.
W odróżnieniu od analogów globalnych, podstawowe równania modeli tej grupy
nie są zależnościami założonymi a priori, lecz wynikają z rozwiązania
5.1. Uwagi wstępne 141
zagadnień klasycznej teorii sprężystości dla ciała izotropowego i jednorodnego.

W wyniku tego parametry E , ν tych modeli są średnimi parametrami
materiałowymi, które można wyznaczyć w laboratorium lub w terenie. Modele
podstawowe wykazują niezłą jakościową zbieżność do wyników pomiarów,
a model Własowa jest szczególnie interesujący z uwagi na łatwość
implementacji numerycznej i możliwość metodycznego doboru grubości
rozpatrywanej warstwy sprężystej.
Trzecia grupa modeli – modele złożone, stanowi rozwinięcie modeli
podstawowych w kierunku możliwie najlepszych aproksymacji rzeczywistego
masywu gruntowego. W ramach liniowej teorii sprężystości należy wyróżnić
tutaj model K-N, który stanowi energetyczny odpowiednik warstwy Własowa.
Umożliwia on rozpatrzenie dowolnego uwarstwienia podłoża, jak również
zmienności modułu odkształcalności wraz z głębokością. Ścisłe potraktowanie
tego problemu wymagałoby zastosowania modelu strefy uwarstwionej,
uwzględniającej różne parametry warstw oraz warunki zszycia na ich stykach,
co prowadziłoby w rezultacie do trójwymiarowej analizy MES. Uwzględnienie
uwarstwienia w modelu K-N jest nieporównywalnie prostsze, a przy tym
wprowadzenie podziału na warstwy, nawet w przypadku podłoża jednorodnego
zwiększa dokładność uzyskanych wyników w porównaniu z warstwą Własowa.
Warto nadmienić, że model Winklera, pomimo swoich mankamentów jest
w dalszym ciągu dość powszechnie wykorzystywany w analizie zagadnień
współpracy konstrukcji z podłożem gruntowym, głównie z uwagi na łatwość
implementacji numerycznej. Jednak tylko wtedy, gdy odkształcalna warstwa
gruntu ma małą grubość, model ten wydaje się być odpowiedni, dając wyniki
zbieżne z rezultatami dla warstwy klasycznej. W chwili obecnej model Winklera
wykorzystywany jest m.in. w takich programach MES jak SAP 2000, Robot,
PlaTo 4.0 czy ABC płyta – wersja podstawowa.
Analizowane w niniejszym rozdziale modele podłoża, włącznie
z przedstawionym w rozdziale 4 modelem podłoża warstwowego, cechuje
sprężystość i liniowość związków „naprężenie-odkształcenie”. Powszechnie
wiadomo, że związki te dla ośrodka gruntowego charakteryzują się silną
nieliniowością; w świetle nowych badań trójosiowych i badań dynamicznych –
także w zakresie małych i bardzo małych odkształceń [4]. Współcześnie
dostrzega się potrzebę wdrażania modeli uwzględniających niesprężystość
i silną nieliniowość ośrodka gruntowego do analiz obliczeniowych geotechniki.
Nie oznacza to jednak całkowitej rezygnacji z sensownych modeli liniowo-
sprężystych, takich jak omówiony wcześniej model warstwowy typu K-N. Jest
też inny ważny powód, by stosować i doskonalić takie modele. Analizy MES
i WMRS zagadnień liniowej sprężystości stanowią wszak powtarzalne, bazowe
ogniwo wszystkich procedur przyrostowo-iteracyjnych, stosowanych w ana-
lizach zagadnień fizycznie nieliniowych, charakterystycznych dla ośrodków
hiposprężystych, sprężysto-plastycznych, itp.
5.2. Belka spoczywająca swobodnie na podłożu

Winklera - WMRS
Rozpatruje się belkę o stałym na długości module sprężystości E i
przekroju poprzecznym A spoczywającą swobodnie na podłożu sprężystym
(rys. 5.2). Zakłada się, że w płaszczyźnie kontaktu nie występują siły tarcia,
a więzy są dwustronne. Poszukuje się rozkładu ugięć i momentów zginających
w belce.
Rys. 5.2. Belka spoczywająca swobodnie na podłożu sprężystym
Równanie osi odkształconej belki (4.1) ma dla rozważanego zagadnienia postać:
d 4 w( x)
EJ = q ( x) − q s ( x) (5.1)
dx 4
gdzie q s (x) jest odporem gruntu (reakcją podłoża).

Warunki brzegowe dla belki spoczywającej swobodnie na podłożu wyrażają się
w zanikaniu sił wewnętrznych (tj. momentu zginającego i siły poprzecznej) na
końcach swobodnych belki. Zachodzi zatem:
d 2w
M x =0 = − EJ x =0 =0
dx 2
d 2w
M x =3l = − EJ x =3l =0
dx 2
(5.2)
d 3w
Q x =0 = − EJ x =0 =0
dx 3
d 3w
Q x =3l = − EJ x =3l =0
dx 3
Rozpatruje się więc zagadnienie brzegowe (1.1), w którym:

5.2. Belka spoczywająca swobodnie na podłożu Winklera 143
d4
L≡
dx 4
(5.3)
f ≡ −( q − q s )
u ≡ w( x)
przy danych warunkach brzegowych opisanych równaniami (5.2).

Model Winklera należy do analogów globalnych, istotą których jest
symulacja zależności „obciążenie-ugięcie” powierzchni podłoża za pomocą
zastępczego układu mechanicznego (sprężyn z membraną, warstwy
ścinanej itp.). Analogi globalne abstrahują od rzeczywistej budowy masywu
gruntowego, oraz od związków „naprężenie-odkształcenie” w jego punktach.
W modelu tym związek między siłą (reakcją podłoża) a ugięciem ma charakter
liniowo-proporcjonalny:
q s = k 0 w( x) (5.4)
gdzie k 0 jest współczynnikiem proporcjonalności (parametrem modelu).

Przy uwzględnieniu związku (5.4) równanie (5.1) przybiera postać:
d 4 w( x)
EJ = q ( x ) − k 0 w( x ) (5.5)
dx 4
Zależność (5.5) jest równaniem osi odkształconej belki spoczywającej na

podłożu Winklera. Funkcjonał równoważny równaniu (5.5), określający
całkowitą energię potencjalną belki na podłożu Winklera i będący sumą energii
potencjalnej belki (2.72), pracy podłoża sprężystego i pracy sił zewnętrznych ma
postać:
2
1 l  d 2 w( x )  1 l
Π C = Π b + Π s + L z = EJ ∫   dx + k 0 ∫ [w( x)] 2dx −
2 0  dx 2  2 0
  (5.6)
l
− ∫ q w( x) dx
0
Funkcjonał (5.6) stanowi podstawę do zastosowania do rozwiązania sformu-

łowanego zagadnienia brzegowego wariacyjnej metody różnic skończonych.
W tym celu dokonuje się podziału belki równomierną siatką dyskretyzacyjną o
boku oczka ∆x = l 2 (rys. 5.3).
Rys. 5.3. Belka z naniesioną siatką dyskretyzacyjną
Dyskretna postać wyrażenia (5.6), po zastąpieniu dla każdego węzła drugiej

pochodnej odpowiednim ilorazem różnicowym przybierze postać:
1  7  wk +1 − 2 wk + wk −1   1
2
7 2
Π C = EJ  ∑  2
 ∆x k  + k 0 ∑ wk ∆x k −
2  k =1  ∆x   2 k =1
 (5.7)
7
− ∑ q k wk ∆x k
k =1
gdzie:
q1 = q 2 = q 3 = q
q 4 = q5 = q 7 = 0
(5.8)
P
q6 =
∆x
W zależności (5.7) ∆x k dla węzłów skrajnych „1” i „7” jest równe 0,5∆x . Dla
pozostałych węzłów ∆x k = ∆x (k = 2,3,...,6) .
Ugięcia fikcyjne w0 i w8 wyznacza się z warunków brzegowych (5.2) 1, 2 :
d 2w w2 − 2 w1 + w0
M1 = 0 ⇒ 2 1 =0⇒ ⇒ w0 = 2 w1 − w2
dx ∆x 2
(5.9)
d 2w w8 − 2 w7 + w6
M7 = 0⇒ 7 =0⇒ = 0 ⇒ 2 w7 − w 6
dx 2 ∆x 2
Równania algebraiczne na wyznaczenie ugięć kolejnych węzłów otrzymuje się

z warunku (2.79), który obecnie przybiera postać:
∂Π C
=0 (i = 1,2, K ,7) (5.10)
∂wi
Przykładowo dla i = 1 (węzeł „1”) na podstawie zależności (5.10) otrzymuje się:
∂Π C 1  w3 − 2 w2 + w1 1  1
= EJ  2 ⋅ 2
⋅ 2  ⋅ ∆x + k 0 2 w1 ⋅ 0,5∆x −
∂w1 2  ∆x ∆x  2 (5.11)
− q ⋅ 0,5∆x = 0
Stąd:
EJ 1 q
4
( w3 − 2 w2 + w1 ) + k 0 w1 = (5.12)
∆x 2 2
Postępując analogicznie dla kolejnych węzłów siatki podziału uzyskuje się układ
równań, który w zapisie macierzowym przedstawia się następująco:
 D1 −2 1 0 0 0 0   w1  q 2
  w   q 
 D2 −4 1 0 0 0  2  
 D3 −4 1 0 0   w3   q 
    
 D3 −4 1 0   w4  = B ⋅  0  (5.13)
 D3 −4 1   w5   0 
    
 D2 − 2  w6  P ∆x 
 sym    
 D1  w7   0 
gdzie:
∆x 4 B
B= , D1 = 1 + k0 , D2 = 5 + B k0 , D3 = 6 + B k 0
EJ 2
Dla danych:
∆x = 1 [ m]
EJ = 1,432 ⋅ 10 5 [N m2 ]
k 0 = 10 7 [ N / m] (5.14)
3
q = 2 ⋅ 10 [ N / m]
P = 10 4 [N ]
otrzymuje się następujące rozwiązanie układu (5.13):
 w1  1,009 
w   
 2  1,378 
 w3  1,777 
    −4
w4  = 2,280 ⋅ 10 [ m] (5.15)
 w   2,975
 5  
w6  3,792 
   
w7  4,455
Na podstawie uzyskanych wartości ugięć (5.15) oblicza się według zależności
(2.25) zapisanej różnicowo dyskretne wartości momentów zginających:
wk +1 − 2 wk + wk −1
M k = − EJ k = 2,3, K ,6 (5.16)
∆x 2
Otrzymuje się:
M 2   − 429 
 M   − 1489 
 3   
M 4  = − 2749 [ Nm} (5.17)
 M   − 1749 
 5  
M 6   2205 
Zobrazowaniem wykonanych obliczeń są wykresy linii ugięcia i momentów

zginających analizowanej belki przedstawione odpowiednio na rys. 5.4 i 5.5.
Rys. 5.4. Rozkład ugięć w belce na podłożu Winklera
Rys. 5.5. Rozkład momentów zginających w belce na podłożu Winklera

5.3. Belka spoczywająca swobodnie na podłożu

Własowa – WMRS
Analizuje się analogiczną jak w podrozdziale 5.2 belkę na podłożu
sprężystym, przy czym jako model podłoża rozpatruje się warstwę Własowa
opisaną równaniem różniczkowym (3.36). Biorąc pod uwagę przyjęty model
podłoża równanie osi odkształconej belki (5.1) przybierze obecnie postać:
d 4 w( x)  d 2 w( x) 
EJ = q ( x) −  − 2t + k w( x) (5.18)
dx 4  dx 2

Parametry k i t modelu Własowa, charakteryzujące pracę podłoża odpowiednio

na ściskanie i ścinanie wyrażają się poprzez parametry materiałowe w sposób
następujący:
h
k = ∫ G z z ( f , z ) 2 dz
0
(5.19)
1h
t = ∫ G f 2 dz
20
E (1 − ν )
gdzie h jest grubością warstwy, G z z = (por. zależność (4.106)),
(1 + ν )(1 + 2ν )
natomiast f jest funkcją definiującą pionowy rozkład przemieszczenia
w obrębie warstwy opisaną zależnością (4.100) dla i = 1 .
Można zauważyć, że podłoże Własowa jest szczególnym przypadkiem
omówionego w podrozdziale 4.5.1 podłoża warstwowego, tj. opisują go związki
dla modelu warstwowego, w którym i = 1 .
W odróżnieniu od modelu Winklera, podłoże Własowa lepiej odzwierciedla
rzeczywistą nieckę osiadań pod fundamentem, co zobrazowano dla przypadku
obciążenia równomiernie rozłożonego na rys. 5.6.
Rys. 5.6. Porównanie rozkładu osiadań dla modelu Winklera i Własowa

z krzywą doświadczalną
5.3. Belka spoczywająca swobodnie na podłożu Własowa 149
Widoczne jest, że w przypadku podłoża Winklera ugięcia rozkładają się tylko

w obrysie fundamentu, podczas gdy dla modelu Własowa ugięcia zanikają
w znacznej odległości od fundamentu. Wynika stąd konieczność uwzględnienia
w obliczeniach konstrukcji spoczywającej na podłożu Własowa oddziaływania
podłoża poza obszarem konstrukcji.
Funkcjonał całkowitej energii potencjalnej układu belka-podłoże Własowa
odpowiadający wyrażeniu (5.6) przyjmie obecnie postać:
2
1 l  d 2 w( x) 
Π C = Π b + Π Wł + L= ∫
2 0  dx 2
1
{
 dx + ∫ k [w( x)] 2 +
  2 lp
(5.20)
2
 dw( x)   l
+t   dx − ∫ q w( x) dx
 dx   0
gdzie przez ∫ oznaczono umownie całkowanie po długości belki i długości

lp
podłoża poza belką.

Funkcjonał (5.20) stanowi punkt wyjścia do zastosowania WMRS do
rozwiązania belki posadowionej na podłożu Własowa. Dokonany w tym celu
podział siatką dyskretyzacyjną belki i podłoża poza belką przedstawiono na
rys. 5.7. Przyjęto przy tym, że ugięcia podłoża zanikają w odległości dwukrotnie
większej od przyjętego boku oczka siatki podziału.
Rys. 5.7. Obszar belki i podłoża poza belką z naniesioną siatką dyskretyzacyjną
Zastępując w zależności (5.20) pochodne odpowiednimi ilorazami różni-

cowymi otrzymuje się dyskretną postać funkcjonału energii potencjalnej układu
belka - podłoże Własowa:
1  7  w − 2w + w  2  1 9
Π C = EJ  ∑  k +1 k
2
k +1
 ∆x k  + k 0 ∑ wk2 ∆x kp +
2  k =1  ∆x   2 k = −1
(5.21)
1  9  w − wk −1   7
2
+ t  ∑  k +1  ∆x kp  − ∑ q k wk ∆x k
2  k =−1  2∆x   k =1

gdzie ∆x1 = ∆x 7 = ∆x −p1 = ∆x9p = 0,5∆x , natomiast dla pozostałych węzłów
∆x k = ∆x kp = ∆x .
Ugięciami fikcyjnymi dla belki, a więc ugięciami węzłów nie należących do
obszaru belki są ugięcia w0 i w8 . Wyznacza się je z warunków brzegowych
analogicznych do warunków (5.9). Otrzymuje się wówczas:
w0 = 2w1 − w2
(5.22)
w8 = 2 w7 − w6
Intensywność obciążenia q1 ÷ q 7 dla kolejnych węzłów określa zależność (5.8),

ponadto zachodzi również:
w− 2 = w−1 = w9 = w10 = 0 (5.23)
co jest wynikiem przyjętego założenia odnośnie zanikania ugięć poza belką.

Uwzględniając w zależności (5.21) związki (5.22) i (5.23) otrzymuje się:
1  6  wk +1 − 2wk + wk −1   1
2
8
ΠC = EJ  ∑  2
 ∆x k  + k 0 ∑ wk2 ∆x kp +
2 k =2  ∆x   2 k =0

2 2 2
1  7  w − wk −1   w   w 
+ t  ∑  k +1  ∆x kp +  0  0,5∆x +  1  ∆x +
2  k =1 2∆x   2∆x   2∆x  (5.24)
 − w7 
2
 − w8 
2 
+  ∆x +   0,5∆x  − (q1 w1 0,5 + q 2 w2 +
 2∆x   2∆x  
+ q 3 w3 + q 6 w6 ) ∆x
Zależność (5.24) oraz zależność (5.10), w której przyjmuje się i = 0,1, K ,10
prowadzi do układu równań, umożliwiającego wyznaczenie ugięć kolejnych
węzłów siatki podziału. Układ ten przedstawia związek (5.25), w którym
przyjęto następujące oznaczenia:
5.3. Belka spoczywająca swobodnie na podłożu Własowa 151
A1 = k 0 + 0,375s, A2 = −0,25s, A3 = 1 + k 0 B + 0,5R, A4 = 1 − 0,25R

A5 = −0,25 R, A6 = 5 + k 0 B + 0,5R, A7 = 6 + k 0 B + 0,5R, A8 = −0,25R
przy czym:
∆x 4 t ∆x 2 t
B= , R= , s= 2
EJ EJ ∆x
 A1 0 A2 0 0 0 0 0 0   w0   0 
 0 A3 − 2 A4 0 0 0 0  
0   w1    0,5q 
  
 A5 − 2 A6 − 4 A4 0 0 0 0  w2   q 
    
 0 A4 − 4 A7 − 4 A4 0 0 0   w3   q 
 0    
0 A4 − 4 A7 − 4 A4 0 0  w4  = B  0  (5.25)
    0 
 0 0 0 A4 − 4 A7 − 4 A4   w5   
 0 0 0 0 A4 − 4 A6 − 2 A8  w6   P ∆x 
    
 0 0 0 0 0 A4 − 2 A3 0   w7   0 
 
 0 0 0 0 0 0 A2 0 A1   w8   0 
Przyjmując dane (5.14) oraz parametr podłoża t = 0,5k 0 otrzymuje się po

rozwiązaniu układu (5.25) następujące wartości ugięć kolejnych węzłów:
w0   1,422 
 w  10,524 
 1  
w2  13,504 
   
 w3  16,408 
    −5
w4  = 19,557  ⋅ 10 [m] (5.26)
 w   23,539 
 5  
w6  27,588
   
w7  29,286
 w8   2,904 
Obliczone na podstawie rozwiązania (5.26) i zależności (5.16) momenty

zginające przedstawiają się natomiast następująco:
M 2   − 108 
 M   − 350 
 3   
 4 
M = − 1129  [ Nm] (5.27)
 M   − 96 
 5  
M 6   3366 
Uzyskane wyniki ugięć i momentów zginających przedstawiono na rysunkach

5.8 i 5.9
Rys. 5.8 Rozkład ugięć belki i podłoża poza belką
Rys. 5.9. Rozkład momentów zginających w belce

5.4. Belka na podłożu Winklera z więzami jednostronnymi 153
Na rysunkach (5.8) i (5.9) zamieszczono także – celem porównania – odpowie-

dnie wartości otrzymane w podrozdziale 5.2 dla belki na podłożu Winklera.
Widoczny jest istotny wpływ przyjętego modelu podłoża na rozkład ugięć
i momentów zginających w belce.
5.4. Belka na podłożu Winklera z więzami

jednostronnymi – WMRS
Rozważa się belkę jak w podrozdziale 5.2 przyjmując obecnie, że więzy są
jednostronne. Oznacza to, że podłoże sprężyste zdolne jest przenosić tylko
obciążenia ściskające i że belka może odrywać się od podłoża. Zagadnienie to,
które w swojej istocie jest zadaniem nieliniowym można numerycznie rozwiązać
postępując drogą kolejnych iteracji. Iteracje te praktycznie dotyczą zwalniania
więzów w punktach, w których wystąpiły ujemne wartości ugięć podłoża, czyli
w miejscach, gdzie konstrukcja wykazuje tendencje do odrywania się od
podłoża. Tym samym więc obszar kontaktu konstrukcji z podłożem zostaje
sprowadzony do obszaru rzeczywistej współpracy podłoża z konstrukcją, jako że
ośrodek gruntowy nie „pracuje” na rozciąganie.
Przyjmuje się, że belka jest obciążona w środku rozpiętości siłą skupioną
(rys. 5.10). Do określenia jest rozkład ugięć w belce przy założeniu więzów
jednostronnych.
Rys. 5.10. Belka na podłożu Winklera
Do rozwiązania zadania zastosowano wariacyjną metodę różnic skończonych.

Przyjęto analogiczną jak w podrozdziale 5.2 dyskretyzację belki (rys. 5.3). Stąd
też równania na wyznaczenie ugięć kolejnych węzłów określa zależność (5.13),
w której obecnie zmianie ulega wektor prawych stron w wyniku innego niż
w podrozdziale 5.2 sposobu obciążenia belki. Wektor ten określony jest
następująco:
0  0 
0  0 
   

 0  
 0 
Fp = B⋅ = B⋅  (5.28)
q 4   P ∆x 
0  0 
   
 0   0 
∆x 4
gdzie B = .
EJ
Przyjmując dane:
EJ = 17 ⋅ 10 4 [ KNm 2 ]
P = 30 [ KN ]
l = 6 [ m] (5.29)
∆x = 1 [m]
k 0 = 25 ⋅ 10 4 [ KN / m 2 ]
i rozwiązując układ równań (5.13) przy uwzględnieniu zależności (5.28)

otrzymuje się następujące wartości ugięć belki na podłożu Winklera:
 w1  − 0,085
w   
 2   0,076 
 w3   0,299 
    −4
w4  =  0,536  ⋅ 10 [m] (5.30)
 w   0,299 
 5  
w6   0,076 
   
w7  − 0,085
Należy zwrócić uwagę, że uzyskane wartości ugięć (5.30) dotyczą więzów

dwustronnych. Widoczne jest, że w węźle „1” i „7” belka wykazuje tendencje do
odrywania od podłoża (ujemne wartości ugięć). W kolejnym kroku należy więc
zwolnić więzy w tych punktach, czyli rozpatrywać w węzłach „1” i „7”
oddzielnie dyskretne równania ugięć belki i podłoża. Prowadzi to do
modyfikacji układu równań (5.13), który przybiera obecnie formę:
5.4. Belka na podłożu Winklera z więzami jednostronnymi 155
 1 −2 1 0 0 0 0 0 0   w1   0 
− 2 D 2 − 4 1 0 0 0 0 0    0 
   w2   
 1 − 4 D3 − 4 1 0 0 0 0   w3   0 
    
 0 1 − 4 D3 − 4 1 0 0 0   w4   P ∆x 
 0 0 1 − 4 D3 − 4 1 0 0   w  = B ⋅  0  (5.31)
  5   0 
 0 0 0 1 − 4 D2 − 2 0 0   w6   
 0 0 0 0 1 −2 1 0 0  w   0 
  7   
 0 0 0 0 0 0 0 k0 0  w1 
p
 0 
  p  0 
 0 0 0 0 0 0 0 0 k0   w7 
przy czym
∆x 4
B= , D 2 = 5 + Bk 0 , D3 = 6 + Bk 0 ,
EJ
natomiast w1 ÷ w7 oznaczają ugięcia belki, z kolei w1p i w7p ugięcia

odpowiednich węzłów podłoża.
Po rozwiązaniu układu (5.31) otrzymuje się następujące wartości :
 w1  − 0,262
   
 w2   0,02 
 w   0,302 
 3  
 w4   0,555 
    −4
 w5  =  0,302  ⋅ 10 [m] (5.32)
 w   0,02 
 6  
 w7  − 0,262
 p  
w1   0 
w p   0 
 2
Uzyskane wartości ugięć (5.32) belki na jednostronnym podłożu Winklera

zobrazowano na rys. 5.11. Na rysunku tym zamieszczono także wyniki
rozwiązania analitycznego [12] (w obliczu symetrii geometrii i obciążenia – dla
połowy rozpiętości belki).
Rys. 5.11. Rozkład ugięć belki na podłożu Winklera z więzami jednostronnymi
Przy przyjętej do obliczeń stosunkowo rzadkiej siatce dyskretyzacyjnej

widoczna jest dobra zgodność rozwiązania wariacyjną metodą różnic
skończonych belki na jednostronnym podłożu Winklera z rozwiązaniem
analitycznym.
5.5. Belka na podłożu dwuwarstwowym - MES

Rozpatrywana w poprzednich przykładach belka posadowiona jest teraz na
podłożu, które jest układem dwóch poziomych warstw o różnych grubościach
i różnych parametrach materiałowych. Belka poddana jest na całej długości
działaniu obciążenia równomiernie rozłożonego o intensywności q (rys. 5.12).
Poszukuje się rozkładu ugięć w belce.
5.5. Belka na podłożu dwuwarstwowym 157
Rys. 5.12. Belka na podłożu dwuwarstwowym
Całkowitą energię potencjalną układu belka – podłoże dwuwarstwowe

opisuje zależność:
Π C = Π b + ΠW + L , (5.33)
w której funkcjonał energii potencjalnej belki Π b oraz potencjał sił

zewnętrznych L określają odpowiednio związki (4.5) i (4.7).
Energia potencjalna podłoża warstwowego Π W (por. zależność (4.109))
przyjmuje dla rozważanego układu dwóch warstw postać:
2
(i )   +
ΠW =
1 2
∑ ∫ ∫
2 (i ) =1 l h (i )
G{ (w f
y y  (i ) (i )

+ w f
(i +1) (i +1) , y 

)
(5.34)
2
1
+ G x(iy)  w f +w
( f (i +1)   dy (i ) dx
)
 ( i ) ( i ) ( i + 1) , x 
2 
gdzie parametry G y(iy) i G x(iy) opisane są zgodnie ze wzorem (4.106) następująco:
G y(iy) =
(
E (i ) 1 − ν (i ) ) , G x(iy) =
E (i )
(5.35)
(1 + ν (i )
)(1 − 2ν )(i )
(
2 1 + ν (i ) )
Poszukuje się funkcji ugięcia w( i ) , (i = 1,2) , które minimalizują funkcjonał

(5.33) i spełniają warunki:
w(1) ( x) ≡ w( x)
lim w(1) ( x) = 0 (5.36)
x →∞
lim w( 2 ) ( x) = 0
x →∞
przy czym w(x) jest ugięciem belki, w(1) i w( 2) są ugięciami odpowiednich

powierzchni granicznych podłoża.
Do rozwiązania sformułowanego zagadnienia brzegowego wykorzystuje
się metodę elementów skończonych. Belka i podłoże dwuwarstwowe zostają
podzielone liniami pionowymi, w wyniku czego zdyskretyzowany układ składa
się z elementów belkowych i makroelementów podłoża, będących zbiorem
dwóch liniowych elementów dwuwęzłowych (rys. 5.13).
Rys. 5.13. Makroelement podłoża warstwowego
Wektor przemieszczeń węzłowych makroelementu opisany zależnością (4.123)

ma obecnie dwie składowe:
[
δ M = δ (1) δ ( 2 ) ]T (5.37)
[
przy czym δ (1) = w(1)1 w(1) 2 ]T
[
, δ ( 2 ) = w( 2 ) 1 w ( 2 ) 2 ] T
.
Zgodnie z zależnością (4.122) związek między przemieszczeniami
makroelementu i przemieszczeniami węzłowymi ma postać:
 w(1)1 
w(1) = [N 1 N 2 ]   = N δ (1)
w(1) 2 
(5.38)
 w( 2 )1 
w ( 2 ) = [N 1 N 2 ]   = N δ ( 2)
w( 2 ) 2 
gdzie funkcje kształtu N 1 i N 2 dla elementu liniowego, dwuwęzłowego opisane

są w układzie lokalnym (rys. 5.13) w następujący sposób:
N 1 = 0,5(1 − ξ )
(5.39)
N 2 = 0,5(1 + ξ )
Funkcjonał (5.33) odniesiony do układu element belkowy – makroelement

podłoża ma postać:
Π Ce = Π be + Π WM + Le (5.40)
Podstawiając do równania (5.40) przyjęte aproksymacje: dla elementu

belkowego (4.28) i dla makroelementu (5.38) oraz uwzględniając warunek
tożsamości ugięć (5.36) 1 otrzymuje się na podstawie związku (3.77) następujące
wyrażenie:
d Π Ce
e
= k be δ e + kWM δ e + Fbe (5.41)
dδ
[
gdzie δ e = w1r ϕxr ]
w( 2 ) r T , r = 1,2 oraz:
k 1W1 k 12
W

k WM = 
k W
21
kW22

przy czym:
 kW1111 0 kW1112  0 kW1211 kW1212 

   
kW11 = 0 0 0  kW12 = 0 0 0 
 11   
kW 21 0 kW1122  12
0 kW 21 kW1222 
(5.42)
 0 0 0  0 0 0 
   
kW21 =  kW2111 0 kW2112  kW22 = 0 kW22 kW2212 
 21   
kW 21 0 kW2222  22
0 kW 21 kW2222 
Wektor Fb dla obciążenia równomiernie rozłożonego odpowiada wektorowi Fbe

określonemu przez zależność (4.40), powiększonemu o składowe zerowe
odpowiednio do wektora δ e . Analogiczne postępowanie należy zastosować w
celu uzyskania macierzy k be ze związku (4.39).
ik
Podmacierze k W (i = 1,2) otrzymuje się na podstawie zależności (4.135).
Przedstawiają się one następująco:
h (1)
ik
kW rm =∫ (G y(1y) N r N m ∫ N (1) , y N ( k ) , y dy (1) +
l 0
h (1)
1
+ G x(1y) N m , x N r , x ∫ N (1) N ( 2) dy (1) +
2 0
(5.43)
H
+ G y( 2y) N r N m ∫ N ( 2 ) , y N ( k ) , y dy ( 2)
+
h (1)
H
+
1 ( 2)
2
G x y N m , x N r , x ∫ N ( 2) N ( k ) dy ( 2) dx ) (r , m = 1,2)
h(1)
Do opisu rozkładu przemieszczeń pionowych wzdłuż grubości warstw (por.

zależność (5.34)) przyjmuje się funkcję liniową w pierwszej i nieliniową
w drugiej warstwie (rys. 5.14):
h (1) − y y
f (1) = (1)
, f ( 2) = (1)
y ∈ 〈 0, h (1) 〉
h h
(5.44)
shγ ( H − y )
f (1) = 0 , f ( 2) = y ∈ 〈h , H 〉
(1)
shγ h ( 2)
przy czym γ jest współczynnikiem zanikania przemieszczeń ( γ = 1,25 H ),

e x − e−x
sh = .
2
Zgodnie z wyrażeniem (4.117) zachodzi:
N (1) = f (1)
(5.45)
N ( 2) = f ( 2)
Rys. 5.14. Przyjęte funkcje rozkładu przemieszczeń pionowych

wzdłuż grubości warstwy
Uwzględniając związek (5.45) w zależności (5.43) oraz obliczając

analitycznie całki w granicach odpowiednich warstw otrzymuje się następującą
ik
postać poszczególnych podmacierzy k W rm :
 G y(1y) 1 h (1)
k 1W1 r m = ∫  (1) N r N m + G x(1y) N r , x N m , x dx )
l h 2 3

 G y(1y) 1 (1) h (1)
k 1W2 r m = k W
21
= ∫  N N + Gx y N r , x N m , x dx )
l h
rm (1) r m
2 6

 G y(1y) (5.46)
22  1 (1) h (1)
kW rm = ∫  (1) r m N N + Gx y N r ,x N m ,x +
l h 2 3

1 ( 2)
+ G y( 2y) A(h ( 2 ) ) N r N m +
2
)
G x y B( h ( 2 ) ) N r , x N m , x dx
sh2γ h ( 2) + 2γ h ( 2) sh 2γ h ( 2) − 2γ h ( 2)
gdzie: A(h ( 2) ) = , B( h ( 2 ) ) = .
4γ sh 2 γ h ( 2) 4γ sh 2 γ h ( 2)
Zależność (5.41) po uwzględnieniu następującego wektora parametrów

węzłowych, w którym zgrupowano parametry dotyczące i -tego węzła:
δ ie c = wi [ ϕxi w( 2 ) i ] T
(5.47)
można zapisać w postaci:
dΠ Ce δ1e c 
= (
k be + k WM )  e c  + Fb
e
(5.48)
dδ ec δ 2 
gdzie:
 12 6 12 6 
 2 2
0 − 3 2
0 
 le le le le 
 4
0 −
6 2
0 
 le le 2 le
 
 0 0 0 0 
 
k be =   (5.49)
 12 6
− 0 
 le 3 le 2 
 4 
 0 
 le 
 
 sym 0 
 
 k1111
0 12
k11 11
k12 0 12 
k12
 
 0 0 0 0 0 
 22
k11 21
k12 0 22 
k12
k WM = 11 12
 (5.50)
 k 22 0 k 22 
 0 0 
 
 sym 22
k 22 
oraz:
[
Fbe = ql e 2 ql e 2 12 0 ql e 2 − ql e 2 12 0 ] T
(5.51)
Związek (5.48) i (3.78) stanowi podstawę do otrzymania równań

równowagi MES układu N elementów belki spoczywającej na podłożu
dwuwarstwowym. W celu przeprowadzenia obliczeń dokonano podziału belki i
podłoża w obszarze i poza obszarem belki na elementy skończone zgodnie z rys.
5.15. Założono przy tym, że ugięcia podłoża zanikają w odległości 1 m od
końców belki. Zrozumiałym jest, że przyjęty podział w praktyce obliczeniowej
jest daleko niewystarczający. Celem jednak prześledzenia kolejnych etapów
rozwiązania poprzestano na przyjętej dyskretyzacji.
Rys. 5.15. Przyjęty do obliczeń obszar belki i podłoża wraz z podziałem

na makroelementy „e1” i „e2”.
Przyjęto następujące dane:
EJ = 10 2 [ MNm 2 ]
l = 2 [ m]
E (1) = 10 [ MPa] 
(1) 
ν = 0,35  glina
(1) 
G = 3,7 [ MPa] 
(1)
h = 2 [ m]
( 2)
E = 36 [ MPa]

ν ( 2)
= 0,3  piasek

G ( 2 ) = 13,85 [ MPa]
h ( 2 ) = 6 [ m]
γ =4 (5.52)
q = 0,012 [ MN / m]
Wobec przyjętych danych (5.52) macierze (5.49) i (5.50) przedstawiają się

następująco:
 1200 600 0 − 1200 600 0 
 400 0 − 600 200 0 
 
 0 0 0 0 
k be =   (5.53)
 1200 − 600 0 
 400 0 
 
 sym 0 
 4,367 0 − 0,67 − 1,517 0 − 2,183

 0 0 0 0 0 

 32,471 − 2,183 0 9,939 
k WM =  (5.54)
 4,367 0 − 0,67 
 0 0 
 
 32,471 
Wektor (5.51) ma natomiast składowe:
Fbe = [ 0,006 0,001 0 0,006 − 0,001 0 ]T (5.55)
Otrzymany na drodze agregacji globalny układ równań MES, po uwzględnieniu

warunków brzegowych (por. rys. 5.15)
w1 = w( 2 )1 = 0, ϕ x3 = 0 (5.56)
ma następującą postać:
1208,734 600 − 1,34 − 1201,517 − 2,183  w2  0,006

  
 400 0 − 600 0   ϕ x 2   0,001
   
 64,942 − 2,183 9,939  w( 2 ) 2  =  0  (5.57)
 
 1204,367 − 0,67   w3  0,006
   
 sym 32,471   w( 2) 3   0 
Rozwiązując układ (5.58) otrzymuje się:

 w2  1,2305 
ϕ   
 x 2  0,0228
    −3
 w( 2 ) 2  = 0,0531 ⋅ 10 [ m] (5.58)
 w  1,2441 
 3   
 w( 2 ) 3  0,0921
Uzyskane dyskretne wartości ugięć belki i podłoża zobrazowano na rys. 5.16.
Rys. 5.16. Ugięcia belki i podłoża
5.6. Płyta wyższego rzędu na podłożu

dwuwarstwowym - MES
Rozważa się jednorodną, izotropową płytę wyższego rzędu spoczywającą
swobodnie na podłożu dwuwarstwowym, poddaną działaniu obciążenia
prostopadłego do powierzchni środkowej q ( x, y ) (rys. 5.17). Analogicznie jak
w poprzednich przykładach przyjmuje się założenie, że w płaszczyźnie kontaktu
nie występują siły tarcia, a płyta nie może odrywać się od podłoża. Poszukuje się
funkcji ugięcia powierzchni środkowej płyty.
Rys. 5.17. Płyta spoczywająca na podłożu dwuwarstwowym
Dla rozpatrywanego układu płyta-podłoże dwuwarstwowe funkcjonał całkowitej

energii potencjalnej jest sumą energii deformacji sprężystej płyty, energii
sprężystej podłoża i pracy sił zewnętrznych:
Π C = Π K + ΠW + L (5.59)
gdzie Π K , Π W i L określone są odpowiednio przez zależności (4.62), (4.109)

i (4.63).
Poszukiwane funkcje ugięcia w(i ) ( x, y ) , które minimalizują funkcjonał (5.59)
powinny spełniać warunki analogiczne do opisanych przez związek (5.36):
w(1) ( x, y ) ≡ w( x, y )
lim w(1) ( x, y ) = 0
x , y →∞ (5.60)
lim w( 2 ) ( x, y ) = 0
x , y →∞
Celem zastosowania MES do rozwiązania sformułowanego zagadnienia

brzegowego górna powierzchnia podłoża, tj. powierzchnia z = 0 (rys. 5.17)
zostaje podzielona na elementy, analogiczne dla płyty i podłoża. Elementem
płytowym jest 8-węzłowy element serendipowski, wyspecyfikowany w pod-
rozdziale 4.3, natomiast makroelement podłoża składa się obecnie z dwóch
8-węzłowych elementów izoparametrycznych, co obrazuje rys. 5.18.
5.6. Płyta wyższego rzędu na podłożu dwuwarstwowym 167
Rys. 5.18. Makroelement podłoża dwuwarstwowego
Wektor przemieszczeń węzłowych makroelementu podłoża opisany zależ-

nością (5.37) ma w rozpatrywanym zagadnieniu składowe:
[
δ (1) = w(1)1 w(1) 2 KK w(1) 8 ] T
(5.61)
δ ( 2) = [w ( 2)1 w( 2) 2 KK w( 2) 8 ] T
Związek (4.122) przybiera natomiast formę:
w(1) = [N 1 N 2 KK N 8 ] δ (1) = N δ (1)

(5.62)
w( 2) = [N 1 N 2 KK N 8 ]δ ( 2 ) = N δ ( 2 )
przy czym funkcje kształtu N 1 ÷ N 8 określone są wzorami (4.68).

Macierz sztywności makroelementu ma zasadniczo podobną budowę jak
ik
przedstawia to zależność (5.42), przy czym obecnie podmacierze k W (i = 1,2)
mają następującą postać (por. związek (4.135)):
h (1)
ik
kW rm = ∫ G ( (1)
zz N r N m ∫ N (1) , z N ( k ) , z dz (1) +
Ω 0
h (1)
1 (1)
+ G α z N m , α N r , α ∫ N (1) N ( 2 ) dz (1) +
2 0
H
( 2)
+G zz N r N m ∫ N ( 2) , z N ( k ) , z dz ( 2) +
h (1)
(5.63)
H
1
+ G α( 2z) N m , α N r , α ∫ N ( 2) N ( k ) dz ( 2) dΩ
2
) r , m = 1,2, K 8
h (1)
α = x, y
gdzie przez ∫ oznaczono umownie całkowanie w obszarze płyty i podłoża,

Ω
natomiast parametry G z z i Gα z opisuje zależność (4.106).
Przyjmując funkcje f (1) i f ( 2) definiujące rozkład przemieszczeń wzdłuż
grubości podłoża według wyrażenia (5.44), przy uwzględnieniu zależności
(5.45) i wykonaniu całkowania analitycznego otrzymuje się podmacierze (5.63)
w następującej formie:
 G (z1z) 1 h (1) 
k 1W1 r m = ∫  (1) N r N m + G α(1)z Nr , α Nm, α  dΩ
Ω h 2 3 
 
 G (z1z) 1 h (1) 
k 1W2 r m 21
= kW = ∫  − (1) N r N m + G α(1)z N r, α Nm, α  dΩ

Ω
rm
2 6 
 h
(5.64)
 G (z1z) 1 h (1)
22
kW = ∫  (1) N r N m + G α(1)z Nr , α Nm, α + G ( 2)
A(h ( 2 ) ) N r N m
Ω h
rm zz
2 3

1 ( 2) 
+ G α z B ( h ( 2 ) ) N r , α N m , α  dΩ
2 
gdzie A(h ( 2 ) ) i B (h ( 2 ) ) opisane są jak w wyrażeniu (5.46).
W podrozdziale 5.5, w którym analizowano belkę na podłożu dwu-

warstwowym, oddziaływanie podłoża poza belką uwzględniono dzieląc obszar
podłoża wykraczający poza obrys belki na elementy skończone. Przyjęcie
rozmiarów tego obszaru ma z reguły charakter intuicyjny, co powoduje, że
warunek (5.36) 2, 3 może być spełniony jedynie z dużym przybliżeniem. Bardziej
dokładne spełnienie tego warunku umożliwiają omówione w podrozdziale 4.4
elementy nieskończone. Elementy te w rozpatrywanym przykładzie wyko-
rzystuje się do uwzględnienia oddziaływania podłoża dwuwarstwowego poza
płytą. Pięciowęzłowe, izoparametryczne elementy pojedynczo nieskończone i

czterowęzłowe, izoparametryczne elementy podwójnie nieskończone (por. 4.4.2)
stosuje się odpowiednio dla krawędzi swobodnych i naroży płyty. Elementy te
pokrywają powierzchnie graniczne warstw w makroelemencie podłoża (rys.
5.19).
Rys. 5.19. Umiejscowienie elementów nieskończonych w zdyskretyzowanym układzie

płyta - podłoże
Dla elementu pojedynczo nieskończonego geometria określona jest przez

zależność (4.90), natomiast pole przemieszczeń w obszarze elementu opisuje
wyrażenie:
 w(1)1 
w 
 (1) 2 
w(1) = [N 1 N 2 K N 5 ]   = N δ (1)
 M 
 w(1) 5 
 
(5.65)
 w( 2 )1 
w 
 ( 2) 2 
w( 2) = [N 1 N 2 K N 5 ]   = N δ( 2)
 M 
 w( 2 ) 5 
 
gdzie N 1 ÷ N 5 są funkcjami kształtu określonymi przez związek (4.68).
Odpowiednio dla elementu podwójnie nieskończonego geometrię i aproksy-

mację ugięcia przedstawia zależność (4.95) oraz wyrażenie:
 w(1)1 
w 
 (1) 2  ~ ~
w(1) = [N 1 N 2 L N 4 ]   = N δ (1)
 M 
w(1) 4 
 
(5.66)
 w( 2 )1 
w 
 ( 2) 2  ~ ~
w( 2) = [N 1 N 2 L N 4 ]   = N δ( 2)
 M 
w( 2 ) 4 
 
przy czym funkcje kształtu N 1 ÷ N 4 opisane są następująco:
N 1 = 0,25 (ξ 2 − ξ ) (η 2 − η ) N 3 = (1 − ξ 2 )(1 − η 2 )
(5.67)
N 2 = 0,5 (1 − ξ 2 ) (η 2 − η ) N 4 = 0,5(ξ 2 − ξ )(1 − η 2 )
Minimalizacja całkowitej energii potencjalnej (5.59) elementu płyty
wyższego rzędu na podłożu dwuwarstwowym prowadzi (por. zależność (3.77))
do związku, który dla obszaru kontaktu płyty z podłożem można zapisać
w postaci:
∂Π Ce k eK 0 δ e  kˆ 1W1 kˆ 1W2   δˆ (1)  F e 

=    +    + 
∂δ e  0 0  0  kˆ W
21
kW 22
 δ ( 2 )   0 
( 48 x 48) ( 48 x1) ( 48 x 48 ) ( 48 x1) ( 48 x1)
(5.68)
e
gdzie k W , δ e i F e określone są odpowiednio przez związki (4.74), (4.65)
22
i (4.76), natomiast podmacierz k W wyraża się zgodnie z zależnością (5.64).
W równaniach (5.68) składowe wektora δ ( 2) są zgodne z opisanymi przez
związek (5.61) 2 , a składowe wektora δ̂ (1) przedstawiają się następująco:
δˆ (1) = [0 0 w1 0 0 0 0 w2 0 0 0 0 K w8 0 0] T (5.69)
przy czym uwzględniono tu warunek (5.60) 1 tożsamości ugięć płyty i górnej

powierzchni ograniczającej podłoże.
Podmacierz kˆ 1W1 , występująca w wyrażeniu (5.68) ma wymiar (40 x 40)

i wynika z powiększenia podmacierzy k 1W1 opisanej przez zależność (5.64)
o elementy zerowe, odpowiednio do budowy wektora δ̂ (1) (por. związek (4.69)).
Analogicznie, podmacierze kˆ 1W2 i kˆ W21 są powiększonymi o wyrazy zerowe
podmacierzami k 1W2 i k W21 o wymiarach odpowiednio (40 x 8) i (8 x 40) .
Uwzględnienie oddziaływania podłoża poza płytą prowadzi natomiast
w odniesieniu do elementu do równań, które ogólnie można zapisać
w następującej postaci:
∂Π WS k 1W1 s k 1W2 s   δ (s1)   F(e1) s 

=  21 
22  s   e 
+
δ (S1)  k W s  
s   ( 2)   
∂ S  
kW δ F( 2) s  (5.70)
δ ( 2) 
ik
gdzie podmacierze k W s (i, k = 1,2) określone są przez związek (5.64),
~
w którym obecnie macierz funkcji kształtu N jest równa N lub N
(por. wyrażenie (5.65) i (5.66)) w zależności od tego, czy element jest
pojedynczo czy podwójnie nieskończony. Wektory δ (s1) , δ (s2) występujące
~
w równaniach (5.70) odpowiadają natomiast wektorom δ (i ) lub δ (i ) (i = 1,2)
opisanym w zależnościach (5.65) i (5.66). Przez F(ei ) s oznaczono z kolei
wektory obciążeń węzłowych elementów usytuowanych na odpowiednich
powierzchniach granicznych warstw podłoża. Są one różne od wektora
zerowego jedynie w przypadku, gdy w obszarze podłoża poza płytą jest
przyłożone obciążenie.
Związki (5.68) i (5.70) stanowią podstawę obliczeń MES płyty na
dwuwarstwowym podłożu sprężystym.
5.6.1. Płyta obciążona równomiernie na całym obszarze

Kwadratowa płyta o boku a = 10 m i grubości h = 0,2594 m poddana jest na
całym obszarze działaniu obciążenia równomiernie rozłożonego o intensywności
q = 0,01 MPa . Dane materiałowe płyty i ośrodka gruntowego są następujące:
E p = 21000 [ MPa]
ν p = 0,15
E (1) = E ( 2 ) = 9,1 [ MPa]
ν (1) = ν ( 2) = 0,3
h (1) = 4 [ m]
h ( 2 ) = 6 [ m] (5.71)
γ = 0,125
Do obliczenia są ugięcia płyty.

Z uwagi na symetrię obliczenia ograniczono do ćwiartki płyty, którą
podzielono na 5 x 5 elementów skończonych (rys. 5.20).
Rys. 5.20. Podział ćwiartki płyty i podłoża na elementy
Pomimo jednorodności podłoża w kierunku pionowym dokonano – w celu

zwiększenia dokładności rozwiązań – podziału podłoża na dwie warstwy, co
zilustrowano na rys. 5.21.
Rys. 5.21. Schemat uwzględnienia warstw podłoża w analizowanym zadaniu
Otrzymane wyniki ugięć przedstawiono dla trzech charakterystycznych

punktów A, B, C płyty (rys. 5.20) w tablicy 5.1. W tablicy tej porównano także
uzyskane wartości z innymi rozwiązaniami numerycznymi. Widoczna jest dość
dobra zgodność otrzymanych rezultatów z wynikami zamieszczonych analiz
numerycznych, przy jednak stosunkowo mniejszym (w zakresie analiz MES)
nakładzie pracy obliczeniowej (por. liczbę elementów).
Tablica 5.1.
Ugięcia płyty spoczywającej na podłożu sprężystym.
Liczba elementów podziału ćwiartki płyty
MES-elementy MES-elementy MES/MEB [1] Półprzestrzeń

Punkt nieskończone podłoża poza sprężysta [2]
płyty [10] płytą [7]
2 warstwy 3 warstwy 1 warstwa 1 warstwa
Ugięcie w x 10 −3 [m]
A 6,18 6,62 6,47 7,3

B 3,97 4,63 4,62 4,5
C 2,25 2,89 2,95 2,8
Liczba elementów ¼ płyty
36 144 32 36
5.6.2. Płyta obciążona centrycznie na małym obszarze

Kwadratowa płyta o boku a = 5 m i grubości h = 0,5 m obciążona jest w środku
na niewielkim poletku obciążeniem równomiernie rozłożonym o intensywności
q = 1 MPa (rys. 5.22).
Rys. 5.22. Płyta obciążona centrycznie na niewielkim poletku
Dane płyty i podłoża są następujące:
b = 0,25 [m]
E p = 20000 [ MPa]
ν p = 0,167
E (1) = E ( 2 ) = 40 [ MPa]
(5.72)
ν (1) = ν ( 2) = 0,3
h (1) = 3 [m]
h ( 2 ) = 7 [ m]
γ = 0,125
Do obliczenia są ugięcia płyty.

Analogicznie jak poprzednio podłoże zostaje podzielone na dwie warstwy
zgodnie ze schematem przedstawionym na rys. 5.21. Wobec symetrii geometrii
płyty i obciążenia obliczenia zostają ograniczone do ćwiartki płyty, dla której
przyjmuje się podział na elementy skończone przedstawiony na rys. 5.23.
Rys. 5.23. Podział ćwiartki płyty na elementy skończone
Celem porównania rozwiązań obliczono również analogiczną płytę spoczy-

wającą na podłożu Winklera. Współczynnik k 0 tego modelu podłoża (por.
podrozdział 5.2) obliczono na podstawie wzoru [11]:
2 Es
k0 =
 H (5.73)
a log1 + 2 
 a
gdzie: E s = E (1) = E ( 2 ) , H – całkowita grubość warstwy sprężystej, tj.

(1) ( 2)
H =h +h .
Uzyskane dla podłoża dwuwarstwowego i modelu Winklera wyniki ugięć
płyty zobrazowano dla przekroju C − C (por. rys.5. 23) na rys. 5.24.
Rys. 5.24. Rozkład ugięć płyty i podłoża poza płytą
Widoczne jest, że maksymalne ugięcia płyty na podłożu Winklera są o 23%

mniejsze w stosunku do odpowiednich ugięć płyty na podłożu dwuwar-
stwowym.
5.7. Płyta fundamentowa na podłożu dowolnie

uwarstwionym - MES
Płyta o wymiarach 48,4 x 21,0 x 0,6 m będąca fundamentem cztero-
kondygnacyjnego budynku o mieszanej konstrukcji składającej się ze ścian
zewnętrznych o grubości 0,5 m i przestrzennych żelbetowych ram szkieletu
wewnętrznego (rys. 5.25) posadowiona jest na podłożu o dużej niejednorodności
poprzecznej i znacznie zróżnicowanych modułach sprężystości poszczególnych
warstw [5]. W podłożu tym ponadto na małej głębokości występuje soczewka
miękkiej gliny pylastej o skomplikowanym kształcie, której grubość zmienia się
znacząco od 5,9 m do 0,0 m . Obrazuje to profil geologiczny podłoża, przed-
stawiony schematycznie w postaci ciągu jednakowo odległych przekrojów
poprzecznych gruntu na rys. 5.26. Obciążenie płyty stanowią siły skupione
5.7. Płyta fundamentowa na podłożu dowolnie uwarstwionym 177
przekazywane na płytę od słupów oraz obciążenie równomiernie rozłożone o

intensywności q = 332 KN / m przekazywane od ścian usztywniających (rys.
5.25). Należy określić rozkład ugięć płyty.
Rys. 5.25. Płyta fundamentowa
Do rozwiązania przedstawionego zagadnienia zginania płyty na podłożu

uwarstwionym wykorzystuje się MES w ujęciu algorytmu sformułowanego
w podrozdziale 5.6. Podłoże przy uwzględnieniu miejscowo usytuowanej
soczewki gliny pylastej zostaje podzielone na 4 warstwy o zmiennej grubości,
odpowiednio do danego profilu geologicznego. Opis zmienności grubości
warstw wzdłuż osi x i y przyjmuje się według zależności (4.102). Przyjmuje
się ponadto stały współczynnik Poissona dla poszczególnych warstw, ogólnie
zawarty w granicach ν = 0,3 ÷ 0,35 . Moduł sprężystości 1. warstwy traktuje się
jako stały na grubości warstwy E (1) = 50 MPa , natomiast w przypadku warstwy
2. i 3. zakłada się, że moduł ten zmienia się liniowo wzdłuż grubości zgodnie z
zależnością (4.139).
Rys. 5.26. Profil geologiczny podłoża
Występujące w zależności (4.139) współczynniki E 0(i ) i E1(i ) ( i = 2,3 ) wyznacza

się dla odpowiedniej warstwy na podstawie liniowej regresji wartości modułów
uzyskanych z pomiarów na różnych głębokościach. Moduły te zostały
oszacowane na podstawie metody ścieżek naprężeń [5]. W przypadku gliny
pylastej wykorzystano doświadczalny wariant Lambego metody, dla piasku
średniego posłużono się analityczną wersją metody ścieżek naprężeń,
wykorzystując znany wzór Duncana-Changa [5]. Otrzymane wartości modułów
odkształcalności na różnych głębokościach, jak również uzyskane na ich
podstawie przy wykorzystaniu metody najmniejszych kwadratów wyrażenia
opisujące zmianę modułu wzdłuż grubości odpowiedniej warstwy, tj. przedsta-
wiono na rys. 5.27.
E ( 2 ) ( z ) = 5,47 − 0,43z
(5.74)
E (3) ( z ) = 43,4 + 6,27 z
Płyta zostaje podzielona na 70 elementów skończonych, odpowiednio do

siatki kolumn szkieletu budynku. Ponadto siatka podziału na elementy zawiera
elementy pod ścianami i w obszarze wsporników (rys. 5.28).
W celu uwzględnienia wpływu sztywności ścian na płytę przyjęto, że

sztywności elementów pod ścianami są znacznie większe w porównaniu do
sztywności pozostałych elementów. Sztywność płyty D = 370 MNm , (moduł
sprężystości płyty E p = 2 ⋅ 10 4 MPa i współczynnik Poissona ν p = 0,167 ),
podczas gdy w obszarze pod ścianami sztywność tę przyjmuje się jako równą
10 4 MNm . Wobec bardzo małego stosunku grubości do rozpiętości płyty
h a = 0,03 jako element płytowy przyjmuje się element płyty Kirchhoffa o 16
stopniach swobody (znany w literaturze jako element BFS16).
Rys. 5.27. Moduły odkształcalności poszczególnych warstw
Makroelement podłoża składa się obecnie z 4 płaskich elementów pokry-

wających powierzchnie graniczne warstw. Dla elementu 1. powierzchni
granicznej przyjmuje się aproksymację ugięć jak dla elementu płytowego. Dla
pozostałych elementów aproksymacja ta jest zgodna z zależnością (4.122).
W wyniku tego makroelement posiada 28 stopni swobody: 16 na powierzchni
kontaktu i 4 na każdej pozostałej powierzchni. Ponadto, w odróżnieniu od
poprzednich rozważań (por. podrozdział 5.5 i 5.6), w których funkcje f (i )
definiujące rozkład ugięć wzdłuż grubości warstwy przyjęto a priori, funkcje te
przyjmuje się obecnie jako funkcje kształtu (por. zależność 4.116):
f (i ) = N (i ) = 0,5(1 − ξ (i ) )
(5.75)
f (i +1) = N (i +1) = 0,5(1 + ξ (i ) ) i = 1,2,3,4
Otrzymany w rezultacie model dyskretny analizowanego układu płyta

fundamentowa-podłoże dowolnie uwarstwione składa się z 70 elementów
skończonych, 38 elementów nieskończonych i 130 węzłów. Globalna macierz

sztywności ma wymiar N = 784 przy szerokości półpasma b = 78 i liczbie
niezerowych elementów w półpaśmie wynoszącej 53372. W wyniku rozwią-
zania układu równań otrzymuje się wartości węzłowe ugięć płyty jak również
wartości węzłowe ugięć kolejnych powierzchni granicznych warstw podłoża.
Uzyskane wyniki ugięć płyty zobrazowano na rys. 5.29, natomiast na rys. 5.30
przedstawiono przykładowo rozkład ugięć powierzchni granicznej czwartej
warstwy podłoża.
Rys. 5.29. Rozkład ugięć w płycie fundamentowej w x 10 −2 [ m]
−2
Rys. 5.30. Rozkład ugięć powierzchni granicznej czwartej warstwy w( 4) x 10 [ m]
w obszarze i poza obszarem płyty
Maksymalne ugięcie płyty występuje w węźle 65 (rys. 5.28) i wynosi

wmax = 6,78 ⋅ 10 −2 m . Widoczny jest również (rys. 5.30) wpływ soczewki
miękkiej gliny pylastej na lokalne zwiększenie osiadań warstw granicznych
podłoża, co odzwierciedla się nawet na dużej głębokości (czwarta warstwa
podłoża).
LITERATURA
Literatura podstawowa
[1] CICHOŃ Cz.: Metody obliczeniowe. Wybrane zagadnienia. Wyd.
Politechniki Świętokrzyskiej, Kielce, 2005, s. 273
[2] CICHOŃ Cz., CECOT W., KROK J., PLUCIŃSKI P.: Metody
komputerowe w liniowej mechanice konstrukcji. Wyd. Politechniki
Krakowskiej, Kraków, 2002, s. 421
[3] GELFAND I. M., FOMIN S.W.: Rachunek wariacyjny. PWN, Warszawa,
1979, s. 268
[4] GŁAZUNOW J.: Metody wariacyjne. Wyd. Elbląskiej Uczelni
Humanistyczno-Ekonomicznej, Elbląg, 2005, s. 404
[5] GRABARSKI A., WRÓBEL I.: Wprowadzenie do metody elementów
skończonych. Wyd. Politechniki Warszawskiej, Warszawa, 2008, s. 126
[6] GRZYMKOWSKI R., KAPUSTA A., NOWAK I., SŁOTA D.: Metody
numeryczne. Zagadnienia początkowo-brzegowe. Gliwice, 2009, s. 194
[7] HUNTER P., PULLAN A.: FEM/BEM Notes. Univ. of Auckland, Dep. of
Eng. Sci., New Zeland, 2001, s. 145
[8] HUTTON D.: Fundamentals of finite element analysis. Mc-Graw Hill,
New York, 2004, s. 488
[9] KLEIBER M. ( red.): Komputerowe metody mechaniki ciał stałych. PWN,
Warszawa, 1995, s. 728
[10] MICHLIN S.G., SMOLICKI C.L.: Metody przybliżone rozwiązywania
równań różniczkowych i całkowych. Warszawa, 1972, s. 383
[11] RAKOWSKI G.: Metoda elementów skończonych. Wybrane problemy.
Wyd. Politechniki Warszawskiej, Warszawa, 1996, s. 154
184 Literatura
[12] RAKOWSKI G., KACPRZYK Z.: Metoda elementów skończonych

w mechanice konstrukcji. Wyd. Politechniki Warszawskiej, Warszawa
2005, s. 452
[13] RUSIŃSKI E., CZMOCHOWSKI J., SMOLNICKI T.: Zaawansowana
metoda elementów skończonych w konstrukcjach nośnych, Wyd.
Politechniki Wrocławskiej, Wrocław 2000, s. 443
[14] SEGERLIND L.J.: Applied finite element analysis. J. Wiley and Sons,
1976, s. 422
[15] STRANG J., FIX G.: An analysis of the finite element method. Prentce-
Hall, Englewoods Cliffs N. J., 1973, s. 349
[16] ZIENKIEWICZ O.C., TAYLOR R.L.: The finite element metod.
Butterworth-Heinemann, Oxford, vol. 1,2000, s. 689
Literatura szczegółowa
[1] ALMEIDA V.S., PAIVA J.B.: A mixed BEM – FEM formulation for
layered soil – superstructure interaction. Eng. Analysis with Boundary
Elements, 28, 2004, s. 1111-1121
[2] FRASER R.A., WARDLE L.J.: Numerical analysis of rectangular rafts on
layered foundations. Geotechnique, 26, 4, 1976, s. 613-630
[3] GRYCZMAŃSKI M.: Analytical and numerical subsoil models for soil –
foundation interaction problems. Studia Geotechnica et Mechanica, 16,
3-4, 1994, s. 29-72
[4] GRYCZMAŃSKI M.: Modele podłoża gruntowego stosowane w
projektowaniu. XX Ogólnopolska Konf. „Warsztat pracy projektanta
konstrukcji”, Wisła-Ustroń, 2005, s. 159-208
[5] GRYCZMAŃSKI M., SADECKA L.: Analysis of raft foundation –
arbitrarily layered subsoil interaction problem. Proc. of X European Conf.
on Soil Mechanics and Foundation Engineering, Florence, 1991, s. 209-
212
[6] KOLAR V., NEMEC I.: Modelling of soil – structure interaction, Elsvier,
North-Holland, 1989, s. 333
[7] KOLAR V., NEMEC I.: Komplexni automatizace vypoctu interakce
technickych objektu. W Vyuzite vypoctovej techniki pri statickom reseni
stavebnych konstrukcji, CSVTS, Brno, XI, 1983, s. 103-128
[8] KUJAWSKI J.: Analiza grubych płyt i tarcz metodą elementów
skończonych. Wyd. Politechniki Białostockiej, Białystok, 1979, s. 225
[9] MARQUES M.M., OWEN D.R.J.: Infinite elements in quasi – static
materially nonlinear problems. Int. J. Num. Meth. Eng., 18, 4, 1984, s.
739-751
Literatura 185
[10] SADECKA L.:A finite/infinite element analysis of thick plates on a

layered foundation. Comp. and Struct., 76, 5, 2000, s. 603-610
[11] SELVADURAI A. P. S.: Elastic analysis of soil – foundation interaction.
Elsvier, Amsterdam, 1979, s. 543
[12] TORBACKI W.: Numerical analysis of beams on unilateral elastic
foundation. Arch. of Materials Science and Eng., 29, 2, 2008, s. 109-112
[13] WANG C. M., REDDY J. N., LEE K. H.: Shear deformable beams and
plates. Elsvier, 2000, s. 296
187
METODA RÓŻNIC SKOŃCZONYCH I METODA ELEMENTÓW SKOŃCZONYCH

W ZAGADNIENIACH MECHANIKI KONSTRUKCJI I PODŁOŻA
Streszczenie
W książce przedstawiono zarówno zarys dwóch uznanych i szeroko stosowanych

metod rozwiązywania zagadnień brzegowych: metody różnic skończonych (MRS) i
metody elementów skończonych (MES), jak również ich zastosowania w zagadnieniach
mechaniki konstrukcji współpracującej z podłożem gruntowym. W opisie MRS i MES
nawiązano do nowoczesnych ujęć matematycznych formułowania teorii tych metod,
przedstawiając ogólne koncepcje metod, nie przytłaczające jednak swoim formalizmem
matematycznym. Praca zawiera pięć rozdziałów i bibliografię zawierającą ogółem 28
pozycji, w tym 16 pozycji książkowych. Trzy pierwsze rozdziały zawierają klasyfikację
metod przybliżonych rozwiązywania zagadnień brzegowych, opis MRS ujmujący
sposób tworzenia schematów różnicowych, ocenę dokładności operatorów różnicowych
oraz aplikację metody i problemy z tym związane w analizie statycznej belek i płyt.
W części książki obejmującej te rozdziały przedstawiono także ideę MES, ogólne zasady
dyskretyzacji obszarów rozwiązania, klasyfikację elementów skończonych pod wzglę-
dem różnorodnych kryteriów oraz sposoby formułowania równań MES zarówno przy
podejściu Ritza jak i Galerkina. Znaczną uwagę poświęcono funkcjom kształtu jako
aspektowi metody mającemu istotny wpływ na otrzymane rozwiązanie. Przedstawiono
przykłady zastosowań MES w analizie statycznej belek i przy rozwiązywaniu równań
różniczkowych II rzędu.
Część druga książki obejmująca rozdziały 4 i 5 zawiera charakterystykę wybra-
nych elementów skończonych, takich jak element belkowy, element płyty wyższego
rzędu, element nieskończony i makroelement podłoża warstwowego oraz liczne
przykłady zastosowań MRS i MES w analizie belek i płyt spoczywających na podłożu
sprężystym przy uwzględnieniu różnych modeli podłoża.
188
FINITE DIFFERENCE METHOD AND FINITE ELEMENT METHOD

IN THE STRUCTURAL AND SOIL MECHANICS PROBLEMS
Summary
In the book the acknowledged and widely used methods for solving boundary-
value problem: Finite Difference Method (FDM) and Finite Element Method (FEM)
have been presented together with their application to the soil-structure interaction
problems. In the description of these methods one refer to the up-to-date mathematical
formulations of the basic ideas of FDM and FEM.
The book contains five chapters and bibliography which include 28 items, in this 16
books. Chapters one to three comprises the classification of approximate methods of
solving the boundary-value problem, the way of forming the difference schemes, the
estimation of difference operator errors and application of the FDM method to the static
analysis of beams and plates. In this part of the book the FEM fundamentals has been
also presented. The main rules of discretization, the finite element classification from
different point of view and the ways of forming the FEM equations based on Ritz and
Galerkin approach have been described. The significant attention has been put on the
shape functions as the aspect of the method which essentially influenced the obtained
results. Two FEM examples have been presented: static analysis of beam and solving
second order differential equation.
The second part of the book enclosing chapters four and five contains
specification of some finite elements, such as beam element, higher-order plate element,
infinite element and macroelement of layered subsoil. It contains also examples of
application of FDM and FEM to the problem of beam and plate resting on the elastic
foundation. This foundation has been modeled by different schemes, from most simple
Winkler foundation to more complicated layered foundation.

SiM - z.258 - Pop PDF

Uploaded by

Copyright:

Available Formats

SiM - z.258 - Pop PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

SiM - z.258 - Pop PDF

Uploaded by

Copyright:

Available Formats

Studia i Monografie z.

Metoda różnic skończonych i metoda

3.6. Funkcje kształtu……………………………………………………. 62

Metoda elementów skończonych i metoda różnic skończonych stanowią

z podłożem dowolnie uwarstwionym. Zastosowaniu MRS i MES do tego

Zagadnienia fizyki matematycznej czy techniki, takie na przykład jak

określonym w obszarze V , z warunkami brzegowymi

danymi na części odpowiednio S1 i S 2 brzegu S = S1 ∪ S 2 .

W metodach analitycznych traktuje się ciało jako continuum

Rys. 1.1. Metody przybliżone rozwiązania równania A(u ) + f = 0

W rozważaniach zawartych w niniejszej książce przedstawione będą

Rys. 1.2. Schemat tworzenia modeli obiektu rzeczywistego

Uogólnienia podstawowej wersji MRS opracowane w ostatnich latach

2.2. Operatory różnicowe

Rys. 2.1. Funkcja ciągła y = f (x) i jej wartości w określonych punktach

Pochodną funkcji f (x) w punkcie x n , która oznaczona zostanie przez

Wyrażenie powyższe można zapisać w następującej formie:

Przybliżeniem tego wyrażenia jest równanie:

Pierwsza różnica centralna zawiera punkty symetrycznie położone względem

a to jest właśnie pierwsza różnica centralna.

Rys. 3.2. Styczna i sieczna do krzywej y = f (x)

Oczywiście im mniejsze ∆x , tym bardziej kąt nachylenia stycznej zbliża się do

Odpowiednie schematy różnicowe można zilustrować graficznie (rys. 2.3).

Rys. 2.3. Operatory różnicowe funkcji jednej zmiennej

Rys. 2.4. Obszar zmienności funkcji wraz z siatką dyskretyzacyjną

Z kolei dla przekroju x = const. uzyskuje się:

Dla pochodnej mieszanej odpowiedni iloraz różnicowy otrzymuje się w ana-

podobnie dla pochodnej czwartego rzędu:

Schematy powyższe można zilustrować graficznie (rys. 2.5).

Rys. 2.5. Operatory różnicowe funkcji dwóch zmiennych

2.3. Ocena dokładności operatorów różnicowych

Rys. 2.6. Styczna i sieczna do krzywej y = f (x) w punkcie x n

Przyjęte schematy różnicowe jedynie aproksymują rzeczywistą wartość

Rozbieżność pomiędzy rzeczywistą wartością pochodnej a aproksymacją

gdzie pochodna zerowego rzędu jest równa f (x) i 0! = 1 .

Rys. 2.7. Funkcja ciągła y = f (x)

Zgodnie z rys. 2.7 dla f1 i f −1 mamy:

Odejmując stronami wyrażenia na f1 i f −1 uzyskuje się:

W przypadku gdy k ≠ 1 (podział nierównomierny) R → 0 gdy ∆x → 0 , a więc

2.4. Zastosowanie MRS w analizie statycznej belek

w którym: w(x) – funkcja ugięcia belki, q (x) – intensywność obciążenia

Rys. 2.8. Belka obustronnie utwierdzona obciążona w środku siłą skupioną

Rys. 2.9. Rozkład momentów zginających – rozwiązanie dokładne

Rozpatruje się więc zagadnienie brzegowe opisane ogólnie równaniem (1.1),

Warunki brzegowe wynikające ze sposobu podparcia belki są następujące:

Związek pomiędzy funkcją momentów zginających M ( x) a funkcją ugięcia

W celu uzyskania rozwiązania zadania metodą różnic skończonych dokonuje się

Rys. 2.10. Przyjęta dyskretyzacja belki

W każdym węźle siatki podziału należy napisać odpowiednie równania

W równaniach powyższych wielkości w1′ , w1″ oraz w3′ , w3 ″ są ugięciami

wielkością poszukiwaną będzie jedynie ugięcie węzła „2”. Ugięcie to można

Zastępując w równaniu (2.29) pierwszą pochodną odpowiednim ilorazem

Przy uwzględnieniu związków (2.27) i (2.30) w wyrażeniu (2.28) równanie

W równaniu (2.31) q 2 jest intensywnością obciążenia równomiernie

Rys. 2.11. Obciążenie równomiernie rozłożone przypadające na węzeł „2”

Stąd otrzymuje się:

Momenty zginające w wybranych węzłach oblicza się na podstawie zależności

Uwzględniając we wzorze (2.34) związki (2.27), (2.30) oraz uzyskane

Jak widać, uzyskane w tym przypadku rozwiązanie przy wykorzystaniu metody