1. Sposoby sformułowania zagadnień brzegowych (lokalne, globalne) i koncepcje

budowy rozwiązań przybliżonych.

Zagadnienie (problem) brzegowe: dany jest obszar

Ω

, w którym poszukiwane jest

rozwiązanie, układ równań różniczkowych cząstkowych oraz warunki początkowo –
brzegowe nałożone na zbiór punktów należących do brzegu

∂Ω obszaru.

Ω

∂Ω

( )

P x

W rozważanym obszarze poszukiwana jest funkcja

( )

u x w każdym punkcie ( )

P x . Można

stosować następujące sformułowania zagadnień brzegowych:

• Sformułowanie lokalne (mocne, silne): szukane jest rozwiązanie układu równań

różniczkowych w każdym z punktów obszaru osobno:

u

f

dla P

u g dla P

=

∈

=

∈

L

B

Ω

∂Ω

∈∂Ω

gdzie

są operatorami różniczkowymi odpowiednio w obszarze i na jego

brzegu. Równanie

nosi nazwę warunków brzegowych. Jeżeli

są one nałożone na funkcję (tzn.

i

L B

u g dla P

=

B

1

≡

B

), noszą nazwę podstawowych warunków

brzegowych Dirichleta, natomiast dowolna kombinacja warunków brzegowych
złożona z pochodnych nosi nazwę nosi nazwę naturalnych warunków brzegowych
Neumanna.

• Sformułowanie globalne : może być formułowanie jako problem optymalizacji

funkcjonału lub jako zasada wariacyjna.

o

Minimalizacja funkcjonału:

1

( )

( , )

( )

2

I u

u u

=

−

B

L u

1

W funkcjonałach energetycznych pierwszy składnik prezentuje energię
wewnętrzną układu, podczas gdy drugi jest równy pracy wykonanej przez siły
zewnętrzne. Nieznana funkcja

może przedstawiać sobą przemieszczenia

, odkształcenia

( )

u P

u

ε , naprężenia σ lub wszystkie z nich. Funkcja

u

realizująca

ekstremum (minimum, punkt stacjonarny) funkcjonału

( )

min ( )

u

I u

jest szukana.

Można rozważać problem optymalizacji funkcjonału bez ograniczeń (w całej
przestrzeni rozwiązań dopuszczalnych) lub z ograniczeniami (ekstremum jest
szukane w podprzestrzeni narzuconych ograniczeń).

o

Zasada wariacyjna

( ,

)

( )

u u

u

dla

u V

∂ =

∂

∂ ∈

B

L

W mechanice powyższe równanie może mieć sens np. zasady prac
wirtualnych. Sformułowanie wariacyjne (tzw. słabe) ma podstawowe
znaczenie przy konstruowaniu rozwiązań przybliżonych. Można go uzyskać ze
sformułowania mocnego w czterech krokach:

ƒ

Przemnożenie równania różniczkowego przez dowolną funkcję (tzw.

funkcja testująca),

ƒ

Przecałkowanie wyniku po rozważanym obszarze Ω ,

ƒ

Całkowanie przez części z wykorzystaniem twierdzenia Greena w celu

zredukowania pochodnych do minimalnego rzędu,

ƒ

Wprowadzenie do funkcjonału warunków brzegowych Neumanna.

Sformułowania globalne wymagają dodatkowego całkowania po obszarze.
Sformułowanie wariacyjne jest ogólniejsze, gdyż możliwe jest w przypadku
wszystkich zagadnień brzegowych, podczas gdy ułożenie funkcjonału możliwe
jest tylko dla niektórych zadań mechaniki, np. dla zadań liniowej sprężystości
(funkcjonał Lagrange’a, Hamiltona, Reissnera, Castigliano, itp.).

• Możliwe są również podejścia mieszane, polegające np. na podziale obszaru Ω na

podobszary, gdzie stosuje się różne sformułowania wraz z odpowiednimi warunkami
ograniczającymi.

Budowa rozwiązania przybliżonego problemu brzegowego zależy przede wszystkim od
wybranej metody dyskretnej. Można wyróżnić dwie główne koncepcje:

• Rozwiązanie dyskretne w postaci kombinacji liniowej współczynników liczbowych

oraz funkcji bazowych:

1 1

2 2

1

( )

( )

( ) ...

( )

( )

n

n n

i i

i

p x

a

x

a

x

a

x

a

x

ϕ

ϕ

ϕ

ϕ

=

=

+

+ +

=

∑

.

Funkcje bazowe (najczęściej: wielomiany, funkcje trygonometryczne, funkcje
specjalne) muszą być liniowo niezależne, odpowiednio ciągłe oraz muszą spełniać
jednorodne warunki brzegowe rozważanego problemu (jednorodne warunki to takie, w
których po prawej stronie stoi 0, (np.

0

0

( ) 0,

( ) 0

I

u x

u x

=

= ). Przy takim zapisie

postaci rozwiązania przybliżonego można szukać budując odpowiednie residua, (czyli
wyrażenia świadczące o spełnieniu przez rozwiązanie przybliżone wyjściowych
równań różniczkowych) odpowiednio w obszarze i na brzegu:

2

( )

( )

,

( )

( )

d

b

x

p x

f

x

p x

g

ε

ε

=

−

=

L

B

−

Funkcjonał ważący powyższe wyrażenia ma postać:

( )

d

d

b

b

I p

w d

w d

ε

ε

Ω

∂Ω

=

Ω +

∫

∫

∂Ω

b

.

Wagi

świadczą o odejściu p(x) od wyniku ścisłego odpowiednio w obszarze i

na jego brzegi. Dla metod residuów ważonych (metoda Bubnowa - Galerkina, metoda
najmniejszych kwadratów, metoda kolokacji) i metod energetycznych (metoda
Rayleigha – Ritza) zakłada się błąd na brzegu

d

w i w

0

b

ε

= (ścisłe spełnienie warunków

brzegowych) i rozważa jedynie

( )

d

d

I p

w d

ε

Ω

=

Ω

∫

. Odmienną koncepcję prezentują

tzw. metody Trefza, w których zakłada się ścisłe spełnienie równania wewnątrz
obszaru a rozwiązań przybliżonych poszukuje na jego brzegu.

• Rozwiązanie dyskretne w wybranych punktach obszaru (lub/i jego brzegu) zwanych

węzłami. W tej koncepcji niezbędna jest dyskretyzacja obszaru (na węzły, elementy
itp.), gdzie zastępuje się wielkości ciągłe wielkościami dyskretnymi. Numeryczne
wyniki dyskretne można aproksymować funkcją ciągłą w ramach tzw. postprocesingu.

Do tych metod należą: metoda różnic skończonych (MRS, zamiana operatorów
różniczkowych na różnicowe, poszukiwanie wartości węzłowych funkcji szukanej,
aproksymacja metodami najmniejszych kwadratów), metoda elementów skończonych
(MES, podział na elementy i aproksymacja funkcjami kształtu) oraz metod elementów
brzegowych (MEB, podział brzegu na odcinki, obliczanie całek brzegowych).

Przykład
Belka swobodnie podparta obciążona obciążeniem ciągłym równomiernie rozłożonym.

Sformułowanie lokalne:

2

2

( )

( )

( )

( )

( )

0

2

d

M x

x

w x

f x

f x

x

l

dx

EJ

=

=

= −

≤

L

≤

1

( )

(2

)

(0) 0

(2 ) 0

2

M x

qx l x

w

w l

=

−

=

=

3

Sformułowanie globalne:

• W postaci funkcjonału:

2

2

0

1

( )

min ( )

( )

[ (

)

] ,

(0)

(2 ) 0

2

l

w

dw

M x

I w

I w

w dx

w

w l

dx

EJ

⇒

=

−

=

∫

=

• W postaci zasady wariacyjnej:

2

2

2

0

( )

[

] ( )

0,

(0)

(2 ) 0

( ) funkcja próbna, odpowiednio ciągla, spelnia warunki brzegowe: v(0)

(2 ) 0

l

dw

M x

v x dx

w

w l

dx

EJ

v x

v l

+

=

=

=

−

=

∫

=

lub po przecałkowaniu przez części (sformułowanie słabe):

2

0

( )

[

( )]

0,

(0)

(2 ) 0,

(0)

(2 ) 0

l

dw dv M x

v x dx

w

w l

v

v l

dx dx

EJ

−

=

=

=

=

∫

=

Rozwiązanie przybliżone dla metod residualnych:

• Funkcje bazowe:

2

1

2

( )

(

2 ),

( )

(

2 )

x

x x

l

x

x x

l

ϕ

ϕ

=

−

=

−

,

• Rozwiązanie próbne:

2

1

2

( )

( )

( )

(

2 )

(

2 )

p x

a

x

b

x

a x x

l

b x x

l

ϕ

ϕ

=

+

=

−

+

−

,

• Residuum w obszarze:

2

2

( )

( )

1

( )

2

(6

4 )

(2

)

2

d p x

M x

x

a b x

l

qx l x

dx

EJ

ε

=

+

=

+

−

−

−

• Dla metody Bubnowa - Galerkina:

2

2

1

0

0

2

2

2

2

0

0

1

( )

( )

0

[2

(6

4 )

(2

)] (

2 )

0

2

,

...

1

( )

( )

0

[2

(6

4 )

(2

)] (

2 )

0

2

l

l

l

l

x

x dx

a b x

l

qx l x x x

l dx

a b

x

x dx

a b x

l

qx l x x x

l dx

ε

ϕ

ε

ϕ

⎧

⎧

⋅

=

+

−

−

−

−

=

⎪

⎪

⎪

⎪

⇒

⇒

⎨

⎨

⎪

⎪

⋅

=

+

−

−

−

−

=

⎪

⎪

⎩

⎩

∫

∫

∫

∫

=

• Dla metody najmniejszych kwadratów:

2

( , )

0

2

2

0

( , )

( ) ( )

min ( , )

( , ) 0

1

( , )

[2

(6

4 )

(2

)]

,

...

2

( , ) 0

l

a b

l

I a b

x

x dx

I a b

I a b

a

I a b

a b x

l

qx l x dx

a b

I a b

b

ε

ε

=

⋅

⇒

∂

⎧

=

⎪⎪∂

=

+

−

−

−

⇒

⇒

=

⎨ ∂

⎪

=

⎪∂

⎩

∫

∫

• Dla metody kolokacji (punkty kolokacji:

1

2

2

,

3

3

l

l

x

x

=

=

):

4

2

1

1

1

1

0

1

2

2

1

2

2

2

0

1

( ) (

)

0

2

(6

4 )

(2

) 0

( ) 0

2

,

...

( ) 0

1

2

(6

4 )

(2

) 0

( ) (

)

0

2

l

l

x

x x dx

a b x

l

qx

l x

x

a b

x

a b x

l

qx

l x

x

x x dx

ε

δ

ε
ε

ε

δ

⎧

⎧

⋅

−

=

⎪

+

−

−

−

=

⎪

=

⎧

⎪

⎪

⇒

⇒

⇒

⎨

⎨

⎨

=

⎩

⎪

⎪

+

−

−

−

=

⋅

−

=

⎪

⎪

⎩

⎩

∫

∫

=

W metodzie różnic skończonych MRS wprowadzono w ramach dyskretyzacji obszaru 3
węzły (patrz: rysunek). Z trzech wartości węzłowych dwie z nich stanowią warunki
brzegowe:

, pozostaje do obliczenia wartość

. Przy sformułowaniu lokalnym

zamianie na operator różnicowy ulega operator różniczkowy na drugą pochodną:

0

2

0

w

w

=

=

1

w

2

0

1

1

1

2

2

II

x l

w

w

w

d w

w

Lw

dx

l

=

−

+

=

≈

=

2

2

. Równania różnicowe generuje się metodą kolokacji

(ścisłe spełnienie równania w węzłach obszaru):

}

}

0

0

2

0

1

2

1

1

1

2

2

1

...

2

w

w

w

ql

Lw

f

w

l

EJ

−

+

=

⇒

=

⇒

=

W sformułowaniu globalnym można ułożyć funkcjonał energii potencjalnej układu. Po jego
dyskretyzacji (całkowanie kwadraturą Newtona-Cotesa między węzłami) otrzymuje się:

2

2

1

0

2

1

0

1

2

0

0

1 1

1 1

2

2

1

( , ,

)

[(

)

(

)

(

)

(

)

]

2

2

w

w

w

w

l

l

I w w w

l

l

M w

M w

M w

M w

l

l

EJ

EJ

−

−

=

+

−

+

−

+

2

Niewiadomą

(oczywiście

) otrzymuje się minimalizując powyższy funkcjonał

względem

:

1

w

0

2

0

w

w

=

=

1

w

1

1

1

( ) 0

...

d

I w

w

dw

=

⇒

=

Przy sformułowaniu wariacyjnym słabym (funkcja testowa: (0)

(2 ) 0

v

v l

=

= ) od razu

otrzymuje się gotowe równanie różnicowe:

1

0

1

0

2

1

2

1

0

0

1 1

1 1

2

2

(

)

(

)

2

2

w

w v

v

w

w v

v

l

l

l

l

M w

M w

M w

M w

l

l

l

l

EJ

EJ

−

−

−

−

+

−

+

−

+

0

= .

Podstawiając

i przyrównując wyrażenie stojące przy

dowolnym do zera otrzymuje się wartość

.

0

2

0

2

0 oraz

0

w

w

v

v

=

=

=

=

1

v

1

w

2. Wybrane metody numeryczne, które mają zastosowanie w komputerowych

modelach mechaniki.

Podczas modelowania komputerowego ważnym etapami są: stworzenie modelu
numerycznego badanego obiektu rzeczywistego, co wiąże się z zamianą wielkości ciągłych na
wielkości dyskretne (punktowe), a następnie z implementacją komputerową algorytmów
działających na tych dyskretnych wielkościach. Niezbędne są tu metody numeryczne, które w
sposób zrozumiały dla maszyny cyfrowej pozwalają na rozwiązywanie podstawowych
problemów algebry i analizy funkcjonalnej.

5

• Rozwiązywanie równań algebraicznych

Równania nieliniowe (z jedną zmienną) pojawiają się w przypadku materiału
nieliniowego (np. w plastyczności) lub mogą mieć pochodzenie geometryczne (np.
uwzględnienie dużych ugięć). Rozwiązuje się je metodami iteracyjnymi, w których na
podstawie znajomości tzw. schematu iteracyjnego

1

( )

n

n

x

f x

+

=

, punktu startowego

0

x

oraz żądanej dokładności rozwiązania otrzymuje się ciąg przybliżeń miejsca zerowego
równania ( ) 0

F x

= . Istnieje kilka metod iteracyjnych, np.:

o

Metoda Newtona (metoda stycznych):

1

0

( )

,

...

'( )

n

n

n

n

F x

x

x

x

F x

+

=

−

= ,

o

Metoda siecznych:

1

1

1

( )

(

),

...,

...

( )

(

)

n

n

n

n

n

n

n

F x

x

x

x

x

x

x

F x

F x

+

−

−

=

−

−

=

=

−

0

1

o

Metoda relaksacji (optymalnie szybka metoda wokół punktu startowego):

0

1

0

0

0

'( )

1

( )

,

...,

1

'( )

1

'( )

n

n

n

f x

x

f x

x

f x

f x

+

=

−

−

−

x

=

Każda z metod ma swoją interpretację geometryczną.

• Rozwiązywanie układów równań algebraicznych.

Konieczność taka występuje właściwie w każdej metodzie dyskretnej, jest zazwyczaj
ostatnim zabiegiem algebraicznym, po rozwiązaniu którego otrzymuje się niewiadome
pierwotne, np. w statyce budowli – przemieszczenia węzłowe układu konstrukcyjnego
(MRS, MES, Metoda sił, Metoda przemieszczeń).
W przypadku układu równań nieliniowych (nieliniowość materiałowa lub
geometryczna) stosuje się najczęściej iteracyjną metodę Newtona - Rhapsona. Dla
rozwiązania układu równań:

( ) 0,

1, 2,...

i

F

i

=

=

x

n

przyjmuje się punkt startowy

obliczeń

(0)

(0)

(0)

0

1

2

(

,

,...,

)

n

x

x

x

=

x

buduje się następujący schemat iteracyjny:

1

( )

( )

( )

n

n

n

n

n

+

⋅

=

⋅

−

J x

x

J x

x

F x

, gdzie:

1

1

1

1

( ) ...

( )

( )

...

...

...

( ) ...

( )

n

n

n

n

d

d

F

F

dx

dx

d

d

F

F

dx

dx

⎡

⎤

⎢

⎥

⎢

⎥

= ⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

⎣

⎦

x

x

J x

x

x

- macierz Jacobiego.

W przypadku układu równań liniowych:

(

) ( 1)

( 1)

(det

0)

n n n

n

×

×

×

=

≠

A x

b

A

metody jego

rozwiązywania można podzielić na eliminacyjne (dają wyniki ścisłe, wymagają
dużego nakładu pracy), iteracyjne (dają wyniki obarczone błędem, łatwe w
implementacji), kombinowane i specjalne (metody macierzy rzadkich, metody analizy
frontalnej). Do metod eliminacyjnych należą m.in. popularna metoda eliminacji
Gaussa (polegająca na doprowadzeniu macierzy A do postaci diagonalnej) i metoda
eliminacji Choleskiego (dla macierzy symetrycznych, dodatnio określonych, polega na

6

rozkładzie macierzy A na macierze samosprzężone: dolno- i górno-trójkątną). Do
metod iteracyjnych zaliczyć można metody iteracji Jacobiego i Gaussa – Seidle’a
będące uogólnieniem metody iteracji prostej dla równań algebraicznych. Powyższe
metody mogą służyć również do odwracania macierzy.

• Rozwiązywanie problemu własnego.

Konieczność taka pojawia się np. w zadaniach stateczności konstrukcji (poszukiwanie
siły krytycznej) lub w zadaniach dynamiki budowli (poszukiwanie częstości drgań
własnych). Algebraicznie problem własny można zapisać następująco:

(

) ( 1)

( 1)

(

) 0,

n n

n

n

λ

λ

λ

×

×

×

⋅

= ⋅

⇒

−

=

∈

A

x

x

x A

I

ℜ

.

Aby powyższy układ miał niezerowe rozwiązanie, (ale niejednoznaczne!), musi być
spełniony warunek:

det(

) 0

λ

−

=

A

I

. Otrzymuje się w ten sposób równanie

wielomianową na niewiadomą

λ , zwane równaniem charakterystycznym. Po jego

rozwiązaniu dostaje się wartości własne ,

1, 2,...,

i

i

n

λ

=

. Wstawiając je kolejno do

wyjściowego równania otrzymuje się n wektorów własnych

,

1, 2,...,

i

i

n

=

x

,

jednoznacznych, co do długości bądź kierunku. W przypadku, gdy wszystkie wartości
własne są równe 0, to zbiorem wektorów własnych jest cała hiperpłaszczyzna.
Problemu własnego nie rozwiązuje się korzystając z definicji, jako iż jest to
czasochłonne. Ponieważ większość problemów własnych z dziedzin mechaniki
dotyczy macierzy symetrycznych, tak więc metody numeryczne są właśnie dla nich
opracowane. Metodą numeryczną można rozwiązywać problem własny kilkoma
wariantami: można szukać wybranej wartości własnej (i odpowiadającego jej wektora
własnego), wartości własnych z wybranego przedziału lub wszystkich wartości
własnych i wektorów własnych (np. metoda Jacobiego). Do pierwszej grupy metod
należy iteracyjna metoda potęgowa i jej modyfikacje. Wykorzystuje ona własności
tzw. ilorazu Rayleigha. Czysta metoda potęgowa prowadzi do znalezienia wartości
własnej największej co do modułu oraz odpowiadającego jej wektora własnego. Jej
algorytm można zapisać następująco:

0

0

0

0

1

1

1

1

1

max

( )

( )

1

1

...

,

0,1, 2,...

,

,

k

k

k

k

k

T

k

k

k

k

k

k

k

k

k

k

v

v

k

λ

λ

λ

λ

λ

ε

ε

λ

+

+

+

+

+

−

−

=

=

⎧

=

=

⎪

⎪

⎪

=

⇒

=

⇒

⎨

⎪

−

⎪

=

=

−

⎪⎩

max

v

x

x

x

x

x

Av

x

v x

v

v

v

Aby znaleźć kolejne wartości własne, można zastosować powyższy algorytm z
odpowiednimi modyfikacjami, np. rozwiązując problem własny metodą potęgową, ale
dla macierzy odwrotnej

-1

A znajdzie się wartość własną najbliższą zeru a przesuwając

dodatkowo widmo (zbiór wartości własnych) macierzy o zadaną liczbę znajdzie się
wartość własną najbliższą danemu przesunięciu .

d%

d%

7

W zagadnieniach mechaniki teorii elementów skończonych ma się też często do
czynienia z tzw. uogólnionym problemem własnym:

(

)

(

)

( 1)

( 1)

,

n n

n n

n

n

λ

λ

×

×

×

×

⋅

= ⋅

⋅

∈

A

x

B

x

ℜ

,

różniącym się od postaci zwykłej obecnością macierzy (z założenia symetrycznej).

B

Najczęściej problem uogólniony sprowadza się do postaci standardowej:

1

T

λ

−

−

⋅ = ⋅

=

⋅ ⋅

y,

L

A

C y

C

L ,

gdzie macierze dolno- i górnotrójkątne (

T

L , L

) pochodzą z rozkładu macierzy

i rozwiązuje dalej jedną z poznanych wcześniej metod Istnieją też

odpowiednie modyfikacje metody potęgowej dla uogólnionego problemu, bez
potrzeby jego liniowej transformacji do postaci zwykłej.

T

=

⋅

B

L L

• Aproksymacja funkcji

Zagadnienie aproksymacji funkcji f(x) (ciągłej – danej wzorem analitycznym bądź
dyskretnej – danej zbiorem punktów) polega na jej zastąpieniu przez inną funkcję p(x)
w danym przedziale

[ ]

,

x

a b

∈

. Poszukiwana funkcja aproksymująca powinna mieć

łatwy wzór, aby jej dalsza analiza była stosunkowo prosta. Zagadnienie aproksymacji
zapisuje się następująco.

0 0

1 1

0

( )

( )

( )

( ) ...

( )

( )

m

m m

i i

i

f x

p x

a

x

a

x

a

x

a

x

ϕ

ϕ

ϕ

ϕ

=

≈

=

+

+ +

=

∑

m oznacza stopień aproksymacji – jest to liczba związana z ilością liniowo
niezależnych funkcji bazowych ( ),

0,1,...,

i

x

i

m

ϕ

=

przyjmowanych z góry

(najczęściej wielomiany, mogą to być także funkcje trygonometryczne lub funkcje
specjalne np. wielomiany Lagrange’a, Czebyszewa, Legendre’a, Bessela, itp.).
Zadanie aproksymacji sprowadza się do znalezienia liczbowych współczynników

, a od strony algebry do rozwiązania układu równań liniowych.

Najczęściej spotykana jest aproksymacja dyskretna, bazująca na danych punktach

,

0,1,...,

i

a

i

=

m

( , ),

1, 2,...,

i

i

x f

i

=

n , do których funkcja p(x) ma się dopasować. Współrzędne

i

x zwane są węzłami, a

i

f - wartościami węzłowymi. Aproksymacja, w której

sumaryczny błąd odejścia funkcji p(x) od zadanych punktów ( , )

i

i

x f podlega

minimalizacji, zwana jest najlepszą aproksymacją. Jeżeli używa się do tego normy

średniokwadratowej (

2

1

n

i

i

x

x

=

=

∑

), to najlepsza aproksymacja zwana jest metodą

najmniejszych kwadratów, co dla jednomianowych funkcji bazowych sprowadza się
do minimalizacji następującego funkcjonału błędu:

0

1

0

1

1

2

0

1

0

1

( , ,...,

)

1

0

...

...

( , ,...,

)

(

)

min

( , ,...,

)

...

...

m

n

m

j

m

j i

i

m

a a

a

i

j

m

a

a

B a a

a

a x

f

B a a

a

a

−

=

=

=

=

=

−

⇒

⇒

=

∑ ∑

8

Można też stosować metodę ważoną polegającą na wprowadzeniu wag

do funkcjonału

,

1, 2,...,

i

w

i

=

n

i

i

1

2

0

1

1

0

( , ,...,

)

(

)

n

m

j

m

j i

i

j

B a a

a

a x

f w

−

=

=

=

−

∑ ∑

regulujących

odejście wielomianu aproksymującego

1

0

( )

m

j

j

i

p x

a x

−

=

=

∑

od oryginalnych wartości

węzłowych

i

f . Im większa waga, tym bliżej będzie przechodzić krzywa danego

punktu.

Szczególnym przypadkiem aproksymacji jest interpolacja, gdzie w ogólności od
krzywej p(x) wymaga się jedynie ścisłego przejścia przez wszystkie dane punkty

( , ),

1, 2,...,

i

i

x f

i

=

n , co można zapisać: ( )

,

1, 2,...,

i

i

p x

f

i

n

=

=

i co implikuje od

razu stopień interpolacji m=n-1. Krzywa interpolacyjna będzie miała wzór:

1

( )

( )

n

i i

i

p x

a

ϕ

=

=

∑

x . Współczynniki ,

1, 2,...,

i

a

i

n

=

znajduje się z następującego

układu równań:

⋅ = F

Φ a

, gdzie:

( ),

,

,

1, 2,...,

ij

j

i

i

i

x

F

f

i j

n

ϕ

Φ =

=

=

. Jeżeli funkcje

bazowe przyjmie się tak, aby spełniały następujący warunek

0,

( )

( ),

( )

1,

i

i

i

j

ij

dla i

j

x

L x

L x

dla i

j

ϕ

δ

≠

⎧

=

=

= ⎨

=

⎩

, to wtedy prawdziwy jest następujący

wzór (tzw. wzór interpolacyjny Lagrange’a):

1

( )

( )

n

i

i

i

p x

f L

=

=

⋅

∑

x

n

n

. Rozwiązywanie

układu równań nie jest wtedy potrzebne. Uogólnieniem interpolacji Lagrange’a jest
tzw. interpolacja l’Hermitte’a, gdzie w węzłach oprócz wartości funkcji mogą być
również dane wartości pochodnych.
Interpolacja wielomianami coraz wyższych stopni (przy coraz większej liczbie danych
punktów) nie jest korzysta ze względu na zachowanie międzywęzłowe funkcji
interpolujących. Dlatego też przy zachowaniu własności interpolacyjnych można
budować tzw. funkcje sklejane (ang. spline) z wielomianów niskich rzędów na
odcinkach wyznaczanych przez kolejne pary sąsiednich węzłów, stosując odpowiednie
kryteria ciągłości dla tych funkcji.
Problem aproksymacji (w tym także interpolacji) funkcji pojawia się zarówno przy
obróbce danych eksperymentalnych jak i numerycznych; często też dokonuje się
aproksymacji funkcji nieznanej np. przy budowie funkcji kształtu w MES lub przy
generacji wzorów różnicowych.

• Numeryczne różniczkowanie i całkowanie

Numeryczne różniczkowanie i całkowanie najczęściej wykonywane przez algorytm
komputerowy numeryczne operacje przy rozwiązywaniu problemów brzegowych.
Numeryczne różniczkowanie polega na znalezieniu na podstawie wartości
dyskretnych funkcji wartości pochodnej w wybranych węzłach. Najprostszym
rozwiązaniem problemu jest więc aproksymacja funkcji i jej zróżniczkowanie, a
następnie obliczenie wartości pochodnej we wskazanym punkcie, ale ze względu na
złożoność takiego rozwiązania stosuje się inne metody generacji tzw. wzorów
różnicowych. Pochodną numeryczną rzędu k przedstawia się jako kombinację liniową
wartości funkcji

i współczynników liczbowych

:

,

1, 2,...,

i

w

i

=

,

1, 2,...,

i

a

i

=

9

( )

1

,

1, 2,...,

(1, 2,..., )

n

k

j

i

i

i

w

a w

k

j

=

=

=

∈

∑

n .

Współczynniki

można znaleźć rozwijając wartości węzłowe w szereg Taylora

wokół węzła centralnego (tego, w którym liczona jest pochodna, czyli (j)), a następnie
przemnożeniu rozwinięć przez odpowiedni współczynnik, dodaniu stronami i

porównaniu z oryginalnym operatorem różniczkowym

i

a

( )

( )

k

i

k

d

w

dx

. Taki sposób noszący

nazwę metody współczynników nieoznaczonych pozwala też na oszacowanie
dokładności wzoru różnicowego poprzez zebranie pierwszych odrzuconych w
rozwinięciach w szereg wyrazów wyższych rzędów. W bezsiatkowej wersji MRS
generacja wzorów różnicowych odbywa się za pomocą techniki metody najmniejszych
ważonych kroczących kwadratów (ang. Moving Weighted Least Squares – MWSL).
Należy również zaznaczyć, iż wzory różnicowe mogą być budowane nie tylko na
samych wartościach funkcji w węzłach, ale też na wartościach pochodnych dowolnych
rzędów. Takie schematy różnicowe noszą nazwę uogólnionych wzorów różnicowych.

Numeryczne całkowanie funkcji może odbywać się na dwa sposoby. Pierwszy,
noszący nazwę wzorów (tzw. kwadratur) Newtona – Cotesa, polega na zastępowaniu
funkcji podcałkowej wielomianami coraz wyższych rzędów. Do tej grupy należą
popularne wzory prostokątów, trapezów i Simpsona:

( )

calka oryginalna

( ) ,

metoda prostokątów

1

[ ( )

( )] ,

metoda trapezów

2

1

[ ( ) 4

(

)

( )] ,

metoda Simpsona

3

2

2

b

a

p

t

S

I

f x dx

I

f a h h b a

I

f a

f b

h h b a

a b

b a

I

f a

f

f b

h h

=

−

⎧

⎪

=

⋅

= −

−

⎪

⎪

=

+

⋅

= −

−

⎨

⎪

+

−

⎪ =

+ ⋅

+

⋅

=

−

⎪⎩

∫

Inną koncepcję prezentują kwadratury Gaussa, które zakładają wartość całki jako
kombinację liniową wartości funkcji podcałkowej f(x) w punktach zwanych punktami
całkowania ,

1, 2,...,

i

x

i

=

N oraz wag ,

1, 2,...,

i

w

i

N

=

:

1

1

2

2

1

( )

( ) ...

( )

( )

N

G

N

N

i

i

i

I

w f x

w f x

w

f x

w f x

=

=

⋅

+

⋅

+ +

⋅

=

⋅

∑

.

Wagi dobierane są według zasady, by wzór całkowania przybliżony był wzorem
ścisłym dla wielomianu możliwe wysokiego stopnia ortogonalnego (z wagą) w
przedziale

[

]

1,1

−

. Dobór odpowiednich wielomianów implikuje rodzaj kwadratur

Gaussa. (np. kwadratury Gaussa - Legendre’a, Gaussa – Laguerra itp.). Miejsca
zerowe i wagi są najczęściej tablicowane. Ilość miejsc zerowych (wag) N świadczy o
punktowości wzoru Gaussa (np. N=2 – wzór 2-punktowy). Przed zastosowaniem
kwadratury Gaussa należy wyjściowy przedział całkowania

[ ]

,

a b

przetransformować

do przedziału

[

]

1,1

−

.

10

W celu podniesienia dokładności wyniku numerycznego całkowania dzieli się
wyjściowy przedział

[ ]

,

a b

na podprzedziały i do każdego z osobna stosuje się

odpowiednio wybrany wzór niskiego rzędu. Następnie wyniki z poszczególnych
podprzedziałów dodaje się do siebie.
Generacja wzorów różnicowych jest podstawą metod różnicowych (np. MRS),
natomiast całkowanie numeryczne występuje niemal w każdej metodzie dyskretnej,
zwłaszcza przy sformułowaniu globalnym, w którym występuje konieczność
całkowania po zadanym obszarze.

• Rozwiązywanie równań różniczkowych

Problemy, w których występują równania różniczkowe (oprócz nich dany jest obszar,
brzeg i odpowiednie warunki) można podzielić na dwie grupy: problemy początkowe
(warunki – początkowe- zadane są w jednym punkcie obszaru na funkcje niewiadomą
i jej kolejne pochodne – klasyczne zagadnienie Cauchy’ego) oraz problemy brzegowe
(warunki – brzegowe – zadane są w kilku punktach obszaru). Rozwiązywanie
problemów brzegowych zostało omówione wyczerpująco w pytaniu 1, tu omówione
zostanie rozwiązywanie numeryczne problemów początkowych. Poniższe metody
opracowane są dla problemów brzegowych rzędu 1-szego, wszystkie wyższe rzędy
można bowiem do niego sprowadzić. Oryginalne zadanie ma postać:

0

0

( , ),

( )

dy

f x y

y x

y

dx

=

=

.

Ogólny zapis poszukiwania wartości dyskretnych funkcji y(x) ma postać wspólną dla
wszystkich metod:

1

1

( , )

n

n

x

n

n

n

x

y

y

f x y dx

y

+

+

=

+

=

+ Δ

∫

n

.

Poszczególne funkcje różnią się między sobą sposobem obliczania przyrostu

n

Δ na

odcinku

[

]

1

,

n

n

x x

+

. I tak: metody jednokrokowe (metoda Eulera, metoda Rungego -

Kutty) bazują przy obliczaniu

n

Δ na jednym punkcie wstecz, podczas gdy metody

wielokrokowe (metoda Adamsa – Bashfortha, metoda Adama – Moultona) bazują na
znajomości kilku punktów wstecz. Najprostsza z nich wszystkich, metoda Eulera,
zakłada stałą postać funkcji podcałkowej na danych odcinku

[

]

1

,

n

n

x x

+

(

):

1

n

n

h x

x

+

=

−

1

( , )

n

n

n

n

y

y

h f x y

+

=

+ ⋅

.

Następne metody aproksymują funkcję podcałkową wielomianami coraz wyższych
rzędów, np. metoda Rungego – Kutty rzędu II:

1

2

1

1

1

( , )

( ,

)

1

(

)

2

n

n

n

n

n

n

K

h f x y

K

h f x y

K

2

y

y

K

K

+

= ⋅

= ⋅

+

=

+

+

lub metoda Adamsa – Bashfortha rzędu II:

1

1

[23

( ) 16

(

) 5

(

)]

12

n

n

n

n

n

h

y

y

f x

f x

f x

+

−

=

+

⋅

− ⋅

+ ⋅

2

−

.

11

Metoda jednokrokowa wysokiego rzędu oraz para metod wielokrokowych (otwarta i
zamknięta, – czyli bazująca na informacjach odpowiednio znanych i nieznanych
a’priori) stanowią aparat do rozwiązywania problemów początkowych zwanych
metodą predyktor - korektor.

3. Idea tworzenia modeli dyskretnych zgodnie z koncepcją MES na przykładzie 1- i

2- wymiarowego układu ciągłego.

Każda metoda dyskretna (w tym także MES) dokonuje dyskretyzacji rozważanego modelu
ciągłego (matematycznego) danego obiektu rzeczywistego. Dyskretyzacji (zamianie problemu
ciągłego na dyskretny – punktowy) podlegają: obszar (pręt, płyta, tarza itp.), funkcja
niewiadoma (poprzez wprowadzenie tzw. wartości węzłowych) oraz warunki brzegowe.
Podstawową cechą MES jest podział obszaru na elementy – elementy skończone. Są to figury
geometryczne, w których skład wchodzą tzw. funkcje kształtu oraz stopnie swobody (w
węzłach – wierzchołkach elementu). Stopnie swobody to najczęściej wartości węzłowe
funkcji niewiadomej, ale mogą one być również dowolną kombinacją jej pochodnych.
Funkcje kształtu to zestaw funkcji liniowo niezależnych przyporządkowanych elementowi w
taki sposób, że ich wartości w węzłach (wierzchołkach) elementu są równe wartościom
węzłowym. Ogólnie funkcje kształtu (zazwyczaj wielomiany Lagrange lub l’Hermitte’a – ale
też np. funkcje trygonometryczne lub funkcje spline) spełniają warunki reproduktywności

(

( )

,

,

1, 2,..., ,

liczba węzlów

i

j

ij

N x

i j

N

N

δ

=

=

=

) oraz konsystentności (

1

( ) 1

N

i

i

N x

=

=

∑

).

Elementy nie mogą na siebie zachodzić, węzły jednego elementu nie mogą znajdować się
wewnątrz drugiego ani na żadnym z jego boków – węzły mogą jedynie pokrywać się w
wierzchołkach. Suma elementów powinna jak najlepiej pokrywać się ze kształtem obszaru.
Topologię układu m MES opisuje się za pomocą tzw. incydencji, czyli zależności między
elementem, jego węzłami i numerami stopni swobody w węzłach. Każdemu elementowi
przypisuje się wiele wielkości wektorowych i macierzowych. Jedną z nich jest macierz
transformacji, informująca o położeniu elementu w stosunku do globalnego układu
współrzędnych.
Podstawowe znaczenie dla dalszej analizy dyskretnej MES ma wybór sformułowania
problemu oraz model numeryczny. Najczęściej stosuje się sformułowanie słabe,
rozpoczynając rozważania od zasady prac wirtualnych lub globalne, gdzie podstawą analizy
jest odpowiedni funkcjonał energetyczny. Model MES może być np. w przypadku zadań
statyki przemieszczeniowy, naprężeniowy lub mieszany - hybrydowy.
W MES postępowanie mające na celu obliczenie globalnych wielkości i uzyskanie
niewiadomych pierwotnych (np. przemieszczeń) i wtórnych (np. naprężeń i odkształceń)
nazywa się syntezą natomiast przeniesienie tych informacji do układów lokalnych nazywa się
analizą.
W metodzie elementów skończonych można wyróżnić następujące etapy analizy konstrukcji:
analiza na poziomie punktu P, przekroju S, ciała C oraz układu U. Synteza zwie się w MES
powrotem do elementu.
Przystępując do rozwiązania zadania za pomocą MES dokonujemy kolejno:

• Dyskretyzacji układu: wprowadzenie węzłów, podział na elementy, wprowadzenie

stopni swobody w węzłach, opis topologii układu (incydencja!),

• Obliczenia wielkości macierzowych i wektorowych na poziomie elementu (w

układach lokalnych): macierzy sztywności, macierzy transformacji (przejścia od
układu lokalnego do globalnego), ew. macierzy mas (dla zadań dynamiki) i macierzy
wstępnych naprężeń (dla zadań stateczności), wektora zastępników elementowych
oraz wektora wstępnych reakcji.

12

• Transformacja wielkości elementowych do ew. układów węzłowych i układu

globalnego,

• Agregacja elementów w układ (złożenie wielkości elementowych zapisanych w

układzie globalnym z pominięciem działania na wielkościach zerowych),

• Budowa układu równań

: K – globalna macierz sztywności, q – wektor

niewiadomych pierwotnych, P – wektor obciążeń węzłowych, R – wektor reakcji
węzłów (niewiadome wtórne),

K q P R

⋅ = +

• Uwzględnienie podstawowych i naturalnych warunków brzegowych,

• Rozwiązanie układu równań, otrzymanie niewiadomych pierwotnych i wtórnych,

• Powrót do elementu, obliczenie wielkości elementowych (reakcji, sił przekrojowych)

na podstawie rozwiązań globalnych (przemieszczeń, reakcji węzłów),

• Ew. proces adaptacyjny związany z obliczeniem poziomu błędu i przyjęcia nowej

siatki elementów.

Przykład 1D – kratownica

Dyskretyzacja (numeracja węzłów, elementów, globalnych stopni swobody, relacje
przylegania elementów do węzłów):

Y,V

X,U

P = 8N

3

1

2

Q

6

Q

5

3

Q

4

Q

3

2

Q

2

Q

1

1

Incydencja:

• element nr 1 - węzły nr 1 i 2 – stopnie swobody 1,2,3,4,

• element nr 2 - węzły nr 2 i 3 – stopnie swobody 3,4,5,6,

• element nr 3 - węzły nr 1 i 3 – stopnie swobody 1,2,5,6.

Dodatkowo określa się położenie elementów w stosunku do globalnego układu
współrzędnych poprzez podanie ich dostaw kierunkowych (sinus i cosinus). Warunki
brzegowe

Przyjęto elementy prętowe - kratowe (2 węzły, 2 stopnie swobody w elemencie, po jednym
na każdy węzeł – związane z poziomym ruchem sztywnym węzłów).

Uwzględnienie siły skupionej w węźle nr 3 odbywa się na poziomie globalnego układu
równań w wektorze obciążeń węzłowych. Uwzględnienie warunków brzegowych implikuje
zerowe przemieszczenia:

i przebudowę układu równań (skreślenie trzech

kolumn i wierszy odpowiadających tym niewiadomym).

1

2

4

0

Q

Q

Q

=

=

=

13

Przykład 2D – tarcza

Tarczę dyskretyzowano elementami tarczowymi (w płaskim stanie naprężenia) CST
(Constant Strain Triangle – element stałego odkształcenia) trójwęzłowymi. Stopniami
swobody w węzłach są ich przemieszczenia poziome i pionowe. Każdy węzeł na 2 stopnie
swobody. Element ma ich 6. Układ ma

2 N

⋅

stopni swobody, gdzie

oznacza liczbę

węzłów.

N

Ponieważ do tarczy przyłożone jest obciążenie ciągłe, będzie ono musiało zostać zamienione
na tzw. zastępniki obciążenia elementowego – siły węzłowe wywołujące ten sam efekt
kinematyczny, co oryginalne obciążenie. Inna nazwa: równoważniki (zastępniki
przeskalowane przez ‘-‘).

q [kN/m]

i

1

2

3

Y,V

X,U

Uwzględnienie warunków brzegowych przebiega podobnie jak w przypadku 1.
Na podstawie obliczonych przemieszczeń węzłów można aproksymować pole przemieszczeń
dla całej tarczy lub poszczególnych elementów za pomocą wprowadzonej wcześniej bazy
wielomianowej funkcji kształtu (funkcje liniowe – opisują przemieszczenia, ich pochodne to
funkcje stałe a mają charakter odkształceń – stąd nazwa: elementy stałego odkształcenia).
Można też stosować elementy wyższego rzędu (np. elementy trójkątne LST – 6-węzłowe
elementy liniowego odkształcenia) lub elementy tarczowe innego kształtu geometrycznego
(najczęściej obok trójkątnych: prostokątne).

4. Zasadnicze podobieństwa i różnice w koncepcjach MES i MRS.

Różnice i podobieństwa zostaną przedstawione na podstawie cech charakterystycznych
obydwu metod dyskretnych. MRS – metoda różnic skończonych jest najstarszą metodą

14

dyskretną do rozwiązywania zadań brzegowych, pochodzi jeszcze sprzed ery komputerów –
była to wersja na siatkach regularnych. Wraz z rozwojem komputeryzacji powstała jej wersja
na siatki dowolnie nieregularne (tzw. bezsiatkowa MRS) należąca do szerokiej grupy metod
bezsiatkowych. MES natomiast jest najbardziej rozpowszechnioną obecnie metodą
komputerową, mającą bardzo szerokie zastosowanie w problemach technicznych.

Sformułowanie
Obydwie metody służą do rozwiązywania problemów brzegowych mechaniki. Mogą one być
różnorako sformułowane. MRS może być równorzędnie użyta w każdym z nich (lokalne,
globalne, słabe), ale najlepiej znana jest jej wersja w sformułowaniu lokalnym (układ równań
różnicowych). Natomiast MES najlepiej pracuje w sformułowaniach wymagających
całkowania po obszarze – przy zasadzie wariacyjnej lub w zagadnieniach minimalizacji
funkcjonału.

Dyskretyzacja
W MRS wszystkie wielkości dyskretne sprowadza się do węzłów. W wyniku działania MRS
otrzymuje się wyniki węzłowe – dopiero ich późniejsza aproksymacja, niezależna od
wcześniejszych zabiegów – daje rozwiązanie w całym obszarze. Węzły w wersji klasycznej
MRS musiały być rozłożone regularnie, w równych odstępach – tworzyły tzw. siatkę
różnicową, natomiast w wersji bezsiatkowej – węzły mogą być rozłożone dowolnie – bez
żadnych dodatkowych warunków typu regularność czy element. Taka cecha pozwala na łatwe
dokładanie węzłów, przesuwanie oraz usuwanie. Topologia siatki w MRS opisywana jest za
pomocą sąsiadów węzłów: mocnych i słabych. Określenia te pochodzą z podziału obszaru na
podobszary przypisywane węzłom (tzw. wielokąty Voronoi). Na tej podstawie dokonywana
też jest triangularyzacja obszaru. W sformułowaniach globalnych potrzebny jest dodatkowy
podział obszaru na podobszary do całkowania, co jest robione na różne sposoby.
Topologię w MES stanowią węzły i elementy razem budujące siatkę MES. To wybór rodzaju
elementu i aproksymacji między jego węzłami decyduje o końcowym wyniku węzłowym.
Elementy skończone stanowią wielką bibliotekę w zależności od potrzeb można wybierać
elementy pod względem kształtu, liczby węzłów, liczby i rodzaju stopni swobody, rzędu i
jakości funkcji kształtu i innych cech. Do generacji siatek w MES stosuje się specjalne
programy zwane generatorami siatek.

Aproksymacja funkcji niewiadomej
W MRS podstawą metody jest zamiana operatorów różniczkowych na operatory różnicowe,
czyli obliczające wartość danej pochodnej na podstawie informacji z otaczających dany węzeł
innych węzłów. Razem ta konfiguracja węzłów wraz z węzłem centralnym nosi nazwę
gwiazdy różnicowej. Są różne kryteria doboru węzłów do gwiazd, np. kryterium
odległościowe – należy przy tym pamiętać, że zły dobór węzłów spowoduje, iż schemat
różnicowy może okazać się osobliwy lub źle określony. W klasycznej MRS (na siatkach
regularnych) schematy różnicowe generuje się najczęściej metodą współczynników
nieoznaczonych. W bezsiatkowej MRS do generacji całego kompletu wzorów różnicowych
od pierwszej pochodnej aż do ostatniego wymaganego rzędu służy metoda zwana Metodą
Najmniejszych Ważonych Kroczących Kwadratów (ang. Moving Weighted Least Squares –
MWSL). Posługuje się ona lokalną aproksymacją funkcji, dla której różnice między
wartościami węzłowymi są regulowane przez funkcje wagowe. Wpływ danego węzła zanika
wraz z jego odległością od węzła centralnego. Wagi osobliwe wprowadzają interpolację.
Ogólnie aproksymację MWLS w gwieździe z węzłem centralnym

0

x opisuje zależność:

15

1

0

0

0

1

( , )

( )

( ) ( )

m

i

t

i

i

u x x

a x

x

a x p x

−

=

=

⋅

=

∑

Współczynniki ,

1, 2,...,

i

a

i

m

=

(m – rząd lokalnej aproksymacji) znajduje się z układu

równań (n – liczba węzłów w gwieździe) :

, gdzie:

2

t

t

P W Pa P W f

=

2

P – macierz wartości interpelantów (jednomiany) w każdym z węzłów,
W – diagonalna macierz wagowa,
f – wektor wartości węzłowych,
a – wektor niewiadomych współczynników,

2

1

1

1

2

2

2

2

2

(

)

3

3

3

2

1

1

1

...

2

!

1

1

1

...

2

!

1

1

1

...

2

!

... ...

...

...

...

1

1

1

...

2

!

m

m

m

n m

m

n

n

n

h

h

h

m

h

h

h

m

P

h

h

h

m

h

h

h

m

×

⎡

⎤

⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

= ⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

⎣

⎦

,

1

1

2

2

0

(

)

( 1)

0

0

0

0

0

0

,

,

0

0

...

0

...

0

0

0

i

i

n n

n

n

n

w

f

w

f

W

f

w

f

×

×

⎡

⎤

⎡ ⎤

⎢

⎥

⎢ ⎥

⎢

⎥

⎢ ⎥ h x x

=

=

=

⎢

⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎣

⎦

⎣ ⎦

−

,

Wtedy współczynniki można wyrazić:

, natomiast aproksymację w

gwieździe:

. Funkcje

noszą nazwę pseudo funkcji

kształtu, dlatego iż nie odtwarzają przy dowolnych wagach wartości węzłowych (

2

1

2

(

)

t

t

a

P W P

P W f

−

=

( )

( )

t

u x

a p N x f

=

=

⋅

%

( ),

1, 2,...,

i

N x

i

m

=

%

( )

i

i

N x

f

≠

%

).

Ich dobór jest podstawą klasyfikacji wszystkich metod bezsiatkowych.

W MES funkcje kształtu zapewniają ciągłą aproksymację w elemencie, w jego wierzchołkach
odtwarzają wartości węzłowe funkcji niewiadomej. Za pomocą funkcji kształtu i ich
pochodnych można opisać dowolne pole fizyczne w elemencie: ( )

( )

u x

N x f

=

⋅ . Dzięki temu

wszystkie rozważania sprowadza się do węzłów, a między nimi (wewnątrz elementu)
wszystko zależy od przyjętej aproksymacji.

Układ równań algebraicznych
W MRS rzadki lub pasmowy, może być symetryczny. Kłopotliwe jest uwzględnianie
warunków brzegowych ze względu na postać funkcji bazowych (ilorazy wielomianów).
Układ równań jest otrzymywany zależnie od sformułowania zagadnienia metodą kolokacji,
Bubnowa – Galerkina albo Rayleigha – Ritza. Dokładne różniczkowanie i całkowanie funkcji
bazowych jest w tej metodzie czasochłonne.
W MES układ jest na ogół pasmowy i symetryczny, uzyskiwany metodą Bubnowa –
Galerkina albo Rayleigha – Ritza. Stosunkowo łatwo oblicza się współczynniki układu
charakterystyczną dla MES metodą agregacji, często połączoną z eliminacją Gaussa (metoda
frontalna). Łatwo wymusza się warunki brzegowe Dirichleta.

Opracowanie wyników (tzw. postprocessing)
W MRS sprowadza się ono do aproksymacji (uciąglenia) danych dyskretnych w węzłach (co
czyni się najczęściej metodami najmniejszych kwadratów – w tym MWLS, gdzie przyjęcie
nieosobliwych wag powoduje dodatkowe wygładzenie wyniku) lub zróżniczkowania
numerycznego wyniku (np. w celu obliczenia odkształceń należy zróżniczkować

16

przemieszczenia węzłowe – tu obserwowana jest tzw. nadzbieżność – lepsza dokładność
pochodnej numerycznej niż samej funkcji wyjściowej).
W MRS gotowy aparat aproksymacji elementowej pozwala na szybkie znajdywanie wartości
funkcji w dowolnym punkcie obszaru.
W razie, gdy konieczna jest dalsza analiza zadania, np. w wyniku adaptacji (polepszenia
dyskretyzacji w miejscu występowania największego błędu) charakterystyczna dla MRS jest
adaptacja typu h, czyli zagęszczenie siatki węzłów, podczas gdy w MES najczęściej podnosi
się stopień lokalnej (elementowej) aproksymacji (adaptacja typu p), gdyż dokładanie nowych
węzłów i elementów wiąże się z kłopotliwą przebudową wyjściowej siatki.
W MES są dobrze rozwinięte metody szacowania błędu. W MRS te techniki są w ciągłym
rozwoju.

5. Klasyfikacja źródeł błędów w metodach komputerowych

Każda metoda komputerowa zakłada budowę ciągu modeli poprzez przyjmowanie różnych
założeń upraszczających pozwalających na zamianę rzeczywistego obiektu na model
komputerowy. Poszczególne etapy modelowania wiążą się więc z budową kolejnych modeli.
W ciągu tego procesu powstają błędy charakterystyczne dla każdego modelu z osobona,
związane z coraz większym upraszczaniem sytuacji wyjściowej dla danego modelu.

Poszczególne etapy modelowania:

• Punktem wyjścia jest OBIEKT RZECZYWISTY ( istniejący bądź projektowany )

będący przedmiotem analizy statycznej lub dynamicznej albo też pewien PROCES
FIZYCZNY odbywający się w określonym środowisku.,

• MODEL MECHANICZNY będący idealizacją rzeczywistego obiektu poprzez

przyjęcie uzasadnionych założeń upraszczających ( dot. materiału, geometrii,
obciążenia, przyszłych warunków eksploatacji), podczas jego tworzenia popełnia się
błąd modelu zwany nieuniknionym; może się wiązać ze złym przyjęciem parametrów
do modelu,

• MODEL MATEMATYCZNY budowany na podstawie modelu mechanicznego,

ujmujący powyższe idealizowane cechy w opisie matematycznym, na ogół w postaci
układu różniczkowych równań cząstkowych lub też funkcjonału albo zasady
wariacyjnej,

• MODEL DYSKRETNY, dyktowany możliwościami obliczeniowymi komputera (

proste operacje algebraiczne i logiczne), polega na przetworzeniu modelu
matematycznego w układ równań algebraicznych; przy jego tworzeniu niezwykle
ważny jest wybór właściwej metody dyskretnej do analizy problemów brzegowych;

• MODEL NUMERYCZNY, wiążący się z właściwym wyborem nowoczesnych

technik numerycznych, służących rozwiązywaniu dużych układów równań
algebraicznych liniowych i nieliniowych, uogólnionych zagadnień własnych, a także
zagadnień programowania liniowego i nieliniowego. Błędy numeryczne to największa
zmora dzisiejszych modeli komputerowych. Najczęściej wiążą się one z tzw. błędami
obcięcia (obcięcie określonej liczby wyrazów szeregu przy rozwinięciu danej funkcji)
i błędami zaokrągleń. Można też mówić o błędzie residuum (kontrola spełnienia przez
dane rozwiązanie wyjściowego równania – zapobiega przed zbieżnością wyniku do

17

liczby niebędącej rozwiązaniem) lub błędzie zbieżności (kontrola zbieżności pozwala
uniknąć rozbieżności procesu – wynik rozbiega się do nieskończoności lub
niestabilności procesu).

• MODEL KOMPUTEROWY, związany z wyborem odpowiedniego sprzętu

komputerowego do obliczeń oraz napisaniu programu w jednym ze znanych języków (
tzw. modelowanie otwarte) w celu kontrolowania implementowanych procedur
numerycznych lub wyborze dostępnego na rynku programu komercyjnego,
traktowanego na ogół przez modelującego jako „czarna skrzynka” ( tzw. modelowanie
zamknięte).

• PRZYGOTOWANIE DANYCH WEJŚCIOWYCH I WYKONANIE OBLICZEŃ

• WERYFIKACJA POSZCZEGÓLNYCH MODELI NA PODSTAWIE

OTRZYMANYCH WYNIKÓW, bardzo ważny etap, który każe patrzeć krytycznie na
otrzymane wyniki z komputera, które mogą być całkowicie błędne ( „śmieci na
wejściu – śmieci na wyjściu”) lub niekiedy nie spełniają wcześniejszych założeń w
poszczególnych modelach. Analiza błędu zakłada istnienie dwóch rodzajów błędów:
a’priori („przed faktem”) i a’posteriori („po fakcie”). Błąd a’priori liczony na
podstawie teorii danej metody pozwala na oszacowanie dokładności wyniku przed
rozwiązaniem zadania (np. w MES – można określić dla zadań liniowych błąd na
podstawie rozmiaru siatki oraz stopnia aproksymacji). Błąd a’posteriori liczony po
rozwiązaniu danego zadania może dotyczyć samego rozwiązania (kontrolna
zbieżności rozwiązania w węzłach wspólnych dla dwóch siatek) lub tzw. residuum,
czyli dokładności spełnienia wyjściowego równania różniczkowego w wybranych
punktach obszaru. Duży poziom błędu residualnego jest najczęściej podstawą do
ponownego rozwiązania tego samego zadania, ale dla innej, lepszej dyskretyzacji.

Należy też zaznaczyć, iż w miarę możliwości obliczenia, zwłaszcza w przypadku
zadań skomplikowanych, należy zawsze weryfikować przez porównanie wyników np.
z różnych metod dyskretnych lub pochodzących z różnych programów
komputerowych.

Przykład : zamodelować ściskany wspornik i przy pomocy MRS obliczyć siłę krytyczną
Eulera.

0

1

f

w

f

w

1

w

0

18

• Etap I : Model mechaniczny wspornika : belka ( ściskanie + zginanie ) zamocowana

jednostronnie, założenia upraszczające :

1. Geometria : pręt o długości L ( ciało 1D ) – rozważania sprowadzone do osi

obojętnej,

2. Materiał : ograniczenie do zakresu sprężystego ( moduł Younga E ),
3. Przekrój : pręt lity, znane wymiary ( charakterystyki : A,I ),
4. Zastosowano liniową teorię sprężystości ( liniowe związki fizyczne i

kinematyczne ), stąd założenie małych przemieszczeń,

5. Obciążenie : rezygnacja z zasady zesztywnienia aby uwzględnić efekt

wyboczenia : warunki równowagi rozpatruje się w konfiguracji odkształconej,
obciążenie jest zachowawcze ( nie zmienia kierunku podczas wyboczenia )

• Etap II : Model matematyczny : zdecydowano się na model matematyczny w

sformułowaniu lokalnym, w postaci równania różniczkowego razem z warunkami
brzegowymi (układ współrzędnych kartezjańskich zaczepiono we wsporniku );
równanie uwzględnia zasadę małych przekrojów i małych przemieszczeń z modelu
mechanicznego :

sformułowanie problemu (stateczność zlinearyzowana) :

l

x

EJ

w

x

M

w

x

f

w

x

f

x

w

dx

d

x

≤

≤

−

=

=

=

0

)

,

(

)

,

(

,

)

,

(

)

(

)

(

2

2

L

gdzie (wzór na moment zginający uwzględnia odkształcenie pręta) :
:

0

)

0

(

'

0

)

0

(

)

(

)

,

(

=

=

−

−

=

w

w

w

f

P

w

x

M

f oznacza strzałkę ugięcia pręta ( dla x = L ).

• Etap III : Model dyskretny : zdecydowano się użyć Metodę Różnic Skończonych,

zdyskretyzowano pręt dwoma węzłami („1” – x = 0, „2” – x = L ), warunki brzegowe i
nieznaną funkcję przemieszczeń punktów pręta zapisano za pomocą przemieszczeń
węzłowych ( w

1

= 0 i w

1

’ = 0 – warunek na pochodną zapisano za pomocą fikcyjnego

węzła „0” – x = -L : w

0

= w

2

, niewiadoma : w

2

= ? ), równanie różniczkowe

zamieniono na różnicowe ( operator na drugą pochodną zastąpiono operatorem
różnicowym centralnym )

w’’ = (1/h^2) * ( w(i-1) – 2*w(i) + w(i+1) )

• Etap IV : Obliczenia : kollokacja w węźle „1” :

( w0 – 2*w1 + w2 ) / (1/L^2) + ( P / EJ )*w1 = P*f / EJ

uwzględnienie warunków podparcia oraz w2 = f prowadzi do równania na siłę
krytyczną :

2*f / (L^2) = P*f / EJ

19

a

stąd :

P = 2*EJ / (L^2)

• Etap V : Weryfikacja wyników :

Analiza błędu : przez porównanie ze znanym wynikiem analitycznym :

Pe= (1/4) * (

π^2) * EJ / (L^2) = 2.467*EJ / (L^2)

E = | Pe – P | / | Pe | = 0.189

Uwaga! Wynik analityczny jest ścisły tylko dla przyjętego wyżej modelu
mechanicznego!

Źródło błędu w modelu mechanicznym : model jest wystarczający dla potrzeb
inżynierskich, wynik dokładniejszy byłby po przedsięwzięciu środków zaradczych :
uwzględnienie dużych ugięć, rozważanie wyboczenia plastycznego ( hipoteza
Karmana ), uwzględnienie rzeczywistej geometrii pręta ( wyboczenie pręta
cienkościennego – hipotezy Własowa ), przeprowadzenie doświadczenia : zbadanie
siły krytycznej.

Źródło błędu w modelu matematycznym : nieuwzględnienie wpływu ściskania
(podwyższenie siły krytycznej) oraz ścinania w okolicach podpory ( hipoteza
Timoszenki ), środki zaradcze : uwzględnienie w równaniu ugięcia również skrócenia
od ściskania ( 1 / ( E*J*L ).

Źródło błędu w modelu dyskretnym : zbyt mała liczba węzłów, niski rząd operatora
różnicowego ( w obszarze i na brzegu ), środki zaradcze : zwiększyć liczbę węzłów (
przy trzech węzłach na belce rozwiązanie wynosi P = 0.95*Pe ), przeprowadzić
analizę adaptacyjną ( typu h, p i mieszane ), zwiększyć rząd operatorów ( np. poprzez
technikę aproksymacji wyższego stopnia ), zmienić metodę ( np. MES – stateczność
zlinearyzowana - porównać wyniki ), odwołać się do danych z eksperymentu.

6. Zachowanie się układów prętowych przy obciążeniach termicznych i

geometrycznych.

Deformację konstrukcji prętowej mogą spowodować różnego rodzaju czynniki, które noszą
nazwę obciążeń zewnętrznych. Wyróżnia się ich trzy zasadnicze rodzaje:

• Obciążenia statyczne (indeks p) czyli wszelkiego rodzaju siły uogólnione,

• Obciążenia termiczne (indeks t) w postaci równomiernego i nierównomiernego

ogrzania elementów konstrukcyjnych,

• Obciążenia geometryczne (indeks g) pochodzące od wszelkiego rodzaju

wymuszeń kinematycznych (np. osiadania podpór, niedokładności montażu
itp.).

Wpływ obciążeń statycznych na pracę konstrukcji zostanie pominięty, gdyż pytanie nie
dotyczy tego problemu bezpośrednio.

20

Obciążenia termiczne
Obciążenia termiczne ( j t

= ) powodują odkształcenia elementarnego wycinka pręta o

długości

, spowodowane zmianą temperatury. Wpływ temperatury może się uwydatnić w

zmianie wartości trzech sił przekrojowych na odcinku

pręta: momentu gnącego M, siły

podłużnej N oraz siły poprzecznej Q. Przyrosty powyższych wielkości wynoszą odpowiednio:

ds

ds

0

,

,

M

N

s

w

t

t

t

t

ds

t ds

h

α

α

Δ − Δ

Δ =

Δ = ⋅ Δ ⋅

Δ = 0.

Q
t

Zmianie mogą ulec, co najwyżej wartości przyrostów momentu zginającego i siły podłużnej
natomiast temperatura nie ma wpływu na wartość zmiany siły poprzecznej – a raczej ten
wpływ jest pomijalny w przyjętym poprzednio modelu prętowym (ciało 1D). Przyrost
momentu zginającego związany jest z nierównomiernym ogrzaniem dwóch stron pręta – inna
temperatura jest po stronie tzw. „spodów” (indeks „s”) a inna po stronie „wierzchów” (indeks
„w”). Rozróżnienie „wierzchów” i „spodów” jest zarówno w mechanice budowli jak i
wytrzymałości materiałów założeniem umownym i związane jest z przyjęciem globalnego
układu współrzędnych – pozwala na rysowanie dodatnich wartości momentu zginającego po
stronie włókien rozciąganych – „spodów”. Podsumowując: nierównomierne ogrzanie
„spodów” i „wierzchów” powoduje powstanie dodatkowego niezerowego momentu na belce i
faktyczne jej ugięcie – w kierunku temperatury wyższej. Natomiast ogrzanie osi obojętnej
pręta powoduje jego wydłużenie (temperatura dodatnia) lub skrócenie (temperatura ujemna).
Efekty te pojawiają się niezależnie od stopnia statycznej niewyznaczalności układu.
Zakłada się, że przyrosty temperatur

t

Δ

zmieniają się liniowo wzdłuż wysokości przekroju

poprzecznego. Ich wartości należy obliczać ze wzorów:

h

0

0

,

,

s

s

m

m

w

m

m

t

t

t

t

t

t

t

t

t

Δ = −

Δ = −

Δ = − ,

w których przez oznaczono tzw. temperaturę montażu konstrukcji, a przez

-

aktualne temperatury włókien odpowiednio dolnych (spody), górnych (wierzchy) i
położonych w osi obojętnej przekroju poprzecznego. Współczynnik

m

t

0

, i

s

w

t t

t

α jest współczynnikiem

rozszerzalności termicznej (ma wymiar

0

1

C

).

Drugim efektem statycznym oprócz powstania naprężeń w przekroju pręta jest wywołanie
przemieszczeń punktów pręta pochodzenia termicznego – nie mylić z odkształceniami.
Przekrój może się odkształcić tylko wtedy, gdy się przemieści, – ale może się przemieścić bez
powstania odkształceń – taka sytuacja ma miejsce w mechanizmach i była omawiana na
przedmiocie „mechanika ogólna” – w postaci ruchu sztywnego. W przypadku działania
obciążeń termicznych tzw. dopełniająca praca wirtualna (zasada prac wirtualnych w
mechanice budowli jest podstawą do wyprowadzenia wszelkiego rodzaju wzorów na
przemieszczenia) uogólnionych sił przekrojowych (

,

,

i

i

i

)

M N Q wynosi dla całego układu:

0

( )

[

( )

(

t

s

w

w

i

u

t

t

L

M s

t

h

δ

α

Δ − Δ

= −

+ Δ

∫

)] .

i

N s ds

Obciążenia geometryczne
Przy działaniu obciążeń geometrycznych ( j g

=

)

w układach statycznie wyznaczalnych nie

powstają żadne odkształcenia, co bynajmniej nie wyklucza występowania przemieszczeń
spowodowanych przez te obciążenia. Wynika to z faktu, że układy takie nie są
przesztywnione i wszelkie wymuszenia kinematyczne sprawiają, że zachowują się one tak jak

21

mechanizmy ruchome, tzn. przemieszczają się bez występowania jakichkolwiek odkształceń
ich ogniw (elementów). Można więc w rozważanym przypadku (dla USW – układów
statycznie wyznaczalnych) zapisać:

0

0,

M

N

Q

g

g

g

Δ =

Δ =

Δ = 0.

Stąd również praca wirtualna odpowiadająca temu typowi obciążeń wynosi:

0.

g

w

L

δ

=

Natomiast w układach statycznie niewyznaczalnych (USN) każda pojedyncza osiadająca
podpora (osiadanie odbywa się po kierunku nieruchomego przemieszczenia uogólnionego)
zwalnia jeden więz kinematyczny, co nie musi koniecznie oznaczać przejścia konstrukcji w
mechanizm. Stąd ogólnie dla USN obciążenia geometryczne wywołują niezerowe
przemieszczenia, naprężenia i odkształcenia.


Praktyczne uwzględnienie obciążeń geometrycznych (dla USW i USN) odbywa się poprzez
końcowy efekt zastosowania zasady prac wirtualnych dla całego układu będący wzorem
pozwalającym obliczyć przemieszczenie dowolnego węzła (tzw. wzór Maxwella – Mohra):

0

1

( )

( )

( )

( )

( )

( )

( )

[

(

n

i

p

i

p

i

p

s

w

i

i

i

k

u

M s M s

N s N s

k Q s Q s

t

t

p

M

EI

EA

GA

h

α

α

=

⋅

⋅

⋅

⋅

Δ − Δ

=

+

+

+

+ Δ

−

∑

∫

)

( )]

ik

k

s

t N s ds

R

Δ

Wzór ten określa wartość poszukiwanego rzeczywistego

i

p o właściwym dla niego wymiarze

fizycznym tylko wtedy, gdy w i-tym stanie pomocniczym (wirtualnym) posłużymy się
odpowiednim wirtualnym, jednostkowym i bezwymiarowym obciążeniem uogólnionym

1

i

P

= . Wzór uwzględnia wpływy od wszystkich trzech rodzajów obciążeń. Ostatni składnik

wzoru – już poza znakiem całki – jest liczbowym uwzględnieniem wielkości osiadań
kinematycznych.

Przykłady

• Rama portalowa z obciążeniem termicznym

A

B

C

D

1

w

t

t

=

2

s

t

t

=

2

t

2

t

Rama z obciążeniem termicznym

(tzw. stan j=t)

A = A’

B

C

D

Rysunek odkształconej ramy

C’

B’

D’

22

Rama, pokazana na rys., poddana jest obciążeniu termicznemu scharakteryzowanemu przez

łókien skrajnych:

C ,

- współczynnik rozszerzalności termicznej:

następujące wartości:
- temperatura montażu

0

5

m

t

C

= −

,

- aktualna temperatura w

0

0

1

2

5

,

15

t

C

t

= +

= +

5

0

1

10

C

α

−

=

.

Wysokość przekroju prostokątnego słupów i rygla ram

0.2

h

m

=

y wynosi

.

amy:

m równym

i przyrostem

.

- pręt BC jest

ie ogrzany – wygięty i w dłużony z następującymi

Rozważane będą obciążenia termiczne każdego z trzech prętów r
- pręty AB oraz DC są równomiernie ogrzane z wydłużeniem termiczny

0

N

t

t

ε

α

= ⋅ Δ

0

0

2

15 ( 5) 20

m

t

t

t

C

Δ = − =

− − =

nierównomiern

y

odkształceniami:

0

oraz

M

N

s

w

t

t

t

t

t

h

κ

α

ε

α

Δ − Δ

= ⋅

= ⋅ Δ ,

związanymi z wartościami przyrostów:

C

= − = − =

− − =

Δ = − = − = − − =
Δ = Δ − Δ

=

−

+

−

=

− −

+

− −

=

+

=

Na rysunku przedstawiono też szkic ramy odkształconej i przemieszczonej.

• Rama portalowa z obciążeniem geometrycznym

y przemieszczonej wydać charakterystyczną cechę statycznie wyznaczalnych

0

2

0

1

0

0

1

2

15 ( 5) 20

5 ( 5) 10

(

) / 2 [(

) (

)]/ 2 {[5 ( 5)] [15 ( 5)]}/ 2 (10 20) / 2 15

s

s

m

m

w

w

m

m

w

s

m

m

t

t

t

t

t

C

t

t

t

t

t

C

t

t

t

t

t

t

t

Δ

Na rysunku ram
układów poddanych obciążeniu geometrycznemu, a mianowicie proste, niezakrzywione i
niewydłużone poszczególne pręty ramy.

A

B

C

D

A

B

C

D

C’

B’

D’

Rama z obciążeniem geometrycznym

(tzw. stan j=g)

Rysunek przemieszczonej ramy

(1 )

A

A

v

Δ

=

(2 )

A

A

ϕ

Δ

=

(3 )

D

D

v

Δ

=

A’

23

7. Drgania własne i wymuszone, rezonans i tłumienie + Klasyfikacja wpływów

zagadnieniach statyki (bądź też stateczności) budowli zakłada się, że obciążenia są

rakter losowy. W

ń dynamicznych należy wymienić przede wszystkim: takie, które

) okresowe,

ednoznaczne opisanie ruchu konstrukcji wymaga znajomości nie tylko jej własności

dynamicznych, charakterystyki dynamiczne konstrukcji, uwzględnienie
wpływów dynamicznych zagadnieniach konstrukcjach budowlanych

W
przykładane do konstrukcji w sposób statyczny, – czyli energia kinetyczna obciążanego
układu jest pomijalnie mała w porównaniu z jego energią potencjalną. W zadaniach, gdzie jej
wpływ nie może być pominięty, należy analizować konstrukcję metodami dynamiki budowli.
Różnorodność obciążeń zmieniających się w czasie, jakie działają na konstrukcje budowlane
w trakcie ich eksploatacji, jest ogromna. Zmianom mogą ulegać wartości tych obciążeń, ich
zwroty kierunki, także punkty przyłożenia. Każda zmiana, jeżeli tylko jest dostatecznie duża i
szybka, wywołuje tzw. efekty dynamiczne związane z ruchem konstrukcji.
Działające na konstrukcję obciążenia mają w większości przypadków cha
dalszych rozważaniach przedstawione zostaną obciążenia deterministyczne o ściśle
określonych parametrach.
Z najważniejszych obciąże
występują w procesach technologicznych i pochodzą od pracy różnych maszyn i urządzeń,
działanie wiatru na budowle wysokie, ruch pojazdów po mostach i dźwigów po torach
podsuwnicowych oraz trzęsienia Ziemi, efekty parasejsmiczne pochodzące od ruchu
drogowego szybkich i ciężkich pojazdów kołowych, odstrzały w kamieniołomach i tzw.
tąpnięcia górotworu na skutek eksploatacji górniczej prowadzonej w kopalniach.
Na wykresie poniżej pokazano typowe obciążenia dynamiczne: a) harmoniczne, b
c) losowe, d) skokowe i e) impulsowe:

t

J
sprężystych, ale również rozkładu jej masy. Konieczna jest także znajomość efektywnego
tłumienia drgań analizowanego układu. O jego wartości decydują zarówno własności tłumiące

t

t

a) b)

c)

t

d)

t

e)

24

samego materiału konstrukcyjnego, jak również parametry charakteryzujące interakcję układu
drgającego z otaczającym go ośrodkiem fizycznym (powietrzem, wodą). Efektywne tłumienie
drgań jest wielkością fizyczną niewyznaczalną na drodze analitycznej. Dlatego też musi ono
być określane eksperymentalnie, za pomocą odpowiedniej aparatury pomiarowej.
Podstawowym problemem dynamiki budowli jest obliczenie tzw. odpowiedzi konstrukcji Q(t)

ajprostszymi modelami obliczeniowymi konstrukcji drgających są modele jednomasowe o

na działanie danego obciążenia zewnętrznego P(t). Przy zadaniach statyki pomijało się efekty
związane ze zmianą energii kinetycznej przy przykładaniu obciążenia. W dynamice to
założenie odrzuca się jednocześnie rozróżniając drgania własne konstrukcji (wywołane
ciężarem konstrukcji) oraz wymuszone (wywołane przez czynniki zewnętrzne – np. drgające
silniki spoczywające na konstrukcji). Elementy powyższych wektorów interpretuje się
odpowiednio jako obciążenie zewnętrzne i przemieszczenie związane z i-tym stopniem
swobody. Macierze M,C,K określają własności modelu odpowiednio: masowe
(bezwładnościowe), tłumiące i sprężyste.

N
prętach nieważkich (tzn. pozbawionych masy). Stanowią one punkt wyjścia do tworzenia
innych modeli, w tym modeli wielomasowych – modeli statycznych uzupełnionych o pewną
liczbę punktowych mas skupionych. Cechą takiego modelu jest tzw. liczba stopni swobody
drgań (LSSD) – określająca w sposób jednoznaczny położenia wszystkich mas układu
drgającego w dowolnej chwili czasu. Dla układu z rysunku powyżej przy założeniu EA

= ∞

jest LSSD = n. Uwzględnienie rzeczywistych własności mechanicznych pręta ( EA

≠ ∞ )

powoduje zmianę modelu na bardziej realny o LSSD = 2n.
Jedną z podstawowych cech układu drgającego jest jego tzw. widmo częstości drgań

ł W. Kolousek. Zaproponował on przyjęcie dla

Q

(t)

P

(t)

własnych. Niskie częstości drgań własnych związane są z dominującym efektem zginania
prętów. Wraz ze wzrostem efektu rozciągania bądź ściskania, częstości rosną i są tym wyższe,
im większa jest sztywność EA. W większości przypadków, z jakimi spotykamy się w praktyce
obliczeniowej najważniejszą rolę odgrywają najniższe częstości drgań własnych. One to,
bowiem najczęściej decydują o występowaniu takich, niekorzystnych dla konstrukcji,
zjawisk, jak rezonans objawiający się nadmiernym wzrostem jej przemieszczeń, a co za tym
idzie – także przyspieszeń i towarzyszących im sił bezwładności powodujących czasem nawet
bardzo znaczne przeciążenie konstrukcji.
Autorem modeli bardziej złożonych by
każdego pręta LSSS =

∞ . Oznaczało to uwzględnienie ciągłego rozkładu masy wzdłuż osi

pręta. W konsekwencji otrzymał dla pojedynczego pręta jego równanie równowagi w postaci
równania różniczkowego cząstkowego. Przejście od analizy pojedynczych prętów do analizy
układu złożonego odbywa się według algorytmu metody przemieszczeń. Natomiast rewolucję
w zakresie budowania modeli dynamicznych wywołała metoda elementów skończonych
MES, gdzie zarówno modele o masach skupionych jak i ciągłe analizowane są na podstawie
równań różniczkowych zwyczajnych.

M,C,K

1

( )

Q t

( )

Q t

( )

Q t

1

m

n

2

EA

= ∞

2

m

m

( )

P t

( )

P t

n

1

2

( )

P t

n

25

Dla dowolnego układu drgającego warunek „równowagi dynamicznej” układu można zapisać

dzie:

•

korzystając z zasady d’Alemberta. Dla układu o jednym stopniu swobody jest ono w postaci:

( )

( )

( )

( ) 0

B t

T t

S t

P t

+

+

+

= ,

g

( )

( )

B t

M Q

••

= −

t

- siła bezwładności d’Alemberta drgającej masy,

•

- si

- siła wym

tąd otrzymuje się następujące równanie różniczkowe ruchu układu:

ła tłumienia występującego w układzie,

( )

( )

T t

C Q t

•

= −

•

- siła reakcji sprężystej układu,

( )

( )

S t

K Q t

= −

•

uszenia (obciążenia) zewnętrznego.

( )

P t

S

M

( )

( )

( )

( )

Q t

C Q t

K Q t

P t

••

•

+

+

=

.

dpowiedź Q(t), wywołana przez dane wymuszenie P(t), zależy, jak wiadomo, od warunków

O

początkowych

( )

Q t

Q

0

0

=

,

( )

Q t

Q

•

•

=

, określających w sposób jednoznaczny stan układu w

chwili początko

0

0

wej

0

t t

= .

Równanie można też

d

prze stawić w postaci:

2

1

( ) 2

( )

( )

( )

Q t

Q t

Q t

P t

M

γ

ω

••

•

+

+

=

,

po wprowadzeniu oznaczeń

2

2

,

C

K

M

M

γ

ω

=

=

a

i korzystając z zadady superpozycji

(t) jako sumę trzech odpowiedzi:

Q a t

Q b t

Q t

•

=

+

+

a trzy niezależne efekty wywołujące drgania:

• Wychylenie początkowe (dla

przedstawić można rozwiązanie Q

0

0

( )

( )

( )

( )

p

Q t

n

0

t t

= ):

0

0

( )

0

Q t

Q

=

≠ ,

• Prędkość początkowa (dla

):

•

0

t t

=

0

0

( )

0

Q t

Q

•

≠

,

=

• Obciążenie zewnętrzne (dla

0

t t

= ): P(t).

yniki można zapisać w postaci:

W

2

2

( )

[cos(

)

sin(

)],

t

d

d

d

d

a t

e

t

t

γ

γ

ω

ω

ω

ω

ω

−

=

+

=

γ

−

,

1

( )

sin(

),

d

d

b t

e

t

γ

ω

ω

−

=

26

0

1

( )

( ) (

)

t

p

Q t

P

b t

d

M

τ

τ τ

=

−

∫

- tzw. całka Duhamela.

Przypadkiem szczególnym jest brak tłumienia (

0

γ

= ). Wtedy otrzymuje się:

0

1

1

( ) cos( ),

( )

sin( ),

( )

( )sin[ (

)]

t

p

a t

t

b t

t

Q t

P

t

d

M

ω

ω

τ

ω

ω

ω

=

=

=

−

∫

τ τ

]

.

W omawianych powyżej wzorach i zależnościach występowały następujące wielkości
fizyczne:

[

/

rad s

ω

- kołowa częstość drgań (własnych),

2

2

1

d

2

ω

ω

γ

ω

ξ

=

−

=

−

- kołowa częstość drgań tłumionych (

/

ξ γ ω

=

),

pozwalające wyznaczyć wielkości z nimi związane:

2

[ ]

T

s

π

ω

=

- okres drgań własnych,

1

[

lub

/ ]

f

Hz

cykl s

T

=

- częstość drgań (własnych).

Powyższe wzory pozwalają obliczyć odpowiedź układu o LSSD = 1 dla dowolnych
warunków początkowych ruchu i dowolnego obciążenia zewnętrznego. Całkę Duhamela
oblicza się zazwyczaj numerycznie, chyba że rozwiązanie udaje się przedstawić w postaci
analitycznej.

W przypadku, gdy do układu znajdującego się w spoczynku (

0

0

( ) 0 , ( ) 0

Q t

Q t

•

=

=

) zostaje, w

chwili , przyłożone obciążenie o stałej wartości

0

p

t

>

P const

=

, otrzymuje się odpowiedź w

postaci:

(

)

( )

sin[

(

)]

p

t

t

d

d t

P

Q t

e

w t

d

M

γ

τ

τ τ

ω

−

−

=

−

∫

,

a dla braku tłumienia (

0

γ

= ):

( )

{1 cos[ (

)]}

st

Q t

Q

t

ω

τ

=

−

−

,

przy czym

2

/

(

st

Q

P K K

M

)

ω

=

=

jest przemieszczeniem analizowanego układu, wywołanym

przez statycznie działanie siły P. Miarą przewyższenia odpowiedzi statycznej przez
dynamiczną jest tzw. współczynnik dynamiczny

μ

wynoszący w rozważanym przypadku:

2

( )

max

2.0

st

t

st

st

Q

Q t

Q

Q

μ

=

=

=

.

Przy występowaniu tłumienia (

0

γ

≠ ) charakterystyczne są dwa stany występujące podczas

odpowiedzi układu w czasie. Pierwszy to zjawisko dążenia odpowiedzi układu do osiągnięcia
pewnego stanu ustalonego. Stanem ustalonym może być stan spoczynkowej równowagi

27

statycznej układu oraz stan równowagi statycznej. Przed osiągnięciem stanu spoczynkowego
odpowiedź Q(t) ma charakter nieustalony, przejściowy (tzw. stan nieustalony). Przechodzenie
od stanu nieustalonego do ustalonego odbywa się dzięki uwzględnianiu rozpraszania energii
drgań do otoczenia na skutek tłumienia obecnego w każdym układzie fizycznym.

Bardzo często występującym wymuszeniem jest wymuszenie harmoniczne ( )

sin( )

P t

P

t

θ

=

.

Wtedy odpowiedzią jest:

(

)

0

1

( )

sin( )

sin[

(

)]

t

t

d

d

Q t

P

e

w t

d

M

γ

τ

θτ

τ

ω

−

−

=

−

∫

τ

,

a dla braku tłumienia (

0

γ

= ):

( )

[sin( )

sin( )]

st

Q t

Q

t

t

μ

θ

ν

ω

=

−

, przy czym:

2

2

1/(1

),

/

/

/(

)

st

m

Q

P K

P M

μ

ν

ν θ ω

ω

=

−

=

=

=

,

Szczególnym przypadkiem jest stan, kiedy różnica częstości drgań własnych i harmonicznych

θ ω

− = Δ

ω maleje, zbliżając się do przypadku granicznego (rezonansowego) θ ω

=

. Taki

charakterystyczny rodzaj drgań nazywa się dudnieniami, które powstają zawsze wtedy, gdy
mamy do czynienia z superpozycją (nakładaniem się) drgań harmonicznych o mało
różniących się częstościach. Wraz z przybliżaniem się do rezonansu mamy

,

1

θ ω ν

≈

≈ , co

oznacza się zbliżanie się amplitudy drgań (maksymalnej odpowiedzi) do sinusoidy o bardzo
długim okresie T

.

Δ

≈ ∞

W przypadku szczególnym

θ ω

=

mamy do czynienia z wymuszeniem ( )

sin( )

P t

P

t

ω

=

,

wywołującym w układzie drgającym tzw. zjawisko rezonansu na jego częstości własnej

ω .

Wtedy cechą charakterystyczną odpowiedzi jest liniowo-asymptotyczny, względem czasu,
wzrost amplitudy do nieskończoności. Nieograniczone narastanie amplitudy nie występuje
jednak nigdy w układach rzeczywistych, a to ze względu n fakt występowania tłumienia.

Drgania własne układu o wielu stopniach swobody

Postać ogólna równania ruchu układu drgającego o wielu (

n) stopniach swobody:

( )

( )

( )

( )

t

t

t

••

•

+

+

=

M Q

C Q

KQ

P t .

Macierze M,C,K są macierzami kwadratowymi o wymiarach (LSSD

× LSSD), natomiast P i

Q

to macierze jednokolumnowe (wektory) o wymiarach (LSSD

×1). Elementy macierzy

sztywności K i bezwładności M oblicza się na podstawie ich następujących definicji
(

):

,

1, 2,...,

i j

LSSD

=

Element

macierzy K jest reakcją i-tego więzu kinematycznego (siłą uogólnioną)

równoważącą działanie obciążenia (wymuszenia) kinematycznego pochodzącego od
jednostkowego przemieszczenia

j-tego więzu;

ij

k

1

j

Δ =

28

Element

macierzy M jest reakcją i-tego więzu kinematycznego (siłą uogólnioną)

równoważącą działanie sił bezwładności d’Alemberta obciążających układ na skuter

jednostkowego przyspieszenia

, wymuszonego za pomocą j-tego więzu kinematycznego.

ij

m

1

j

Q

••

=

Elementy macierzy tłumienia C oblicza się najczęściej ze wzoru

α

β

=

+

C

M

K .

Wyznaczenie stałych ,

α β

odbywa się na drodze eksperymentalnej.

Dla drgań własnych układu, nietłumionych jest:

0, ( ) 0

t

=

≡

C

P

, co prowadzi do równania

ruchu:

( )

( ) 0

t

t

••

+

=

M Q

KQ

.

Podstawiając możliwą odpowiedź układu w postaci ( )

sin(

)

t

t

ω ϕ

=

+

Q

A

otrzymujemy

równanie:

2

( -

)

0

ω

=

K

M A

reprezentujące problem własny dynamiki budowli. Z algebry wiadomo, iż

ma on niezerowe rozwiązanie, gdy

2

Det( -

) 0

ω

=

K

M

. Po rozwinięciu wyznacznika

otrzymuje się równanie wielomianowe mające dokładnie n dodatnich i rzeczywistych
pierwiastków zwanych częstościami (kołowymi) drgań własnych układu. Zbiór tych
częstości, uporządkowany niemalejąco:

1

2

0

...

n

ω ω

ω

<

≤

≤ ≤

nazywa się widmem częstości drgań własnych układu. Z równania problemu własnego można
otrzymać, po podstawieniu

i

ω ω

=

, odpowiednie niezerowe rozwiązanie

i

A , określające i-tą

postać drgań własnych, która jest i-tym wektorem własnym problemu. Pełne rozwiązanie
problemu własnego tworzy zbiór n tzw. par własnych ( , )

i

i

ω

A . Po wprowadzeniu dwóch

dodatkowych macierzy: spektralnej (widmowej)

Ω i modalnej (własnej)

V

:

1

2

n

ω

ω

ω

⎡

⎤

⎢

⎥

⎢

⎥

=

⎢

⎥

⎢

⎥

⎣

⎦

O

Ω

,

,

1

2

[ ,

,...,

]

n

=

V

A A

A

można sformułować tzw. warunek ortogonalności:

T

V MV = M , w którym macierz M jest

macierzą diagonalną o elementach

i

m , nazywaną modalną macierzą mas.

W zapisie problemu własnego zamiast macierzy sztywności K można stosować macierz
podatności D, o następującej własności: KD = DK = I. Wtedy problem własny przybiera
postać:

2

2

( -

)

0 lub (

-

)

0,

1/

ω

λ

λ

=

=

I

DM A

DM

I A

ω

=

.

29

Rozwiązanie problemu własnego stanowi klucz do utworzenia ważnego algorytmu
dekompozycji modalnej przy obliczaniu odpowiedzi Q(t). Algorytm ten nazywa się analizą
modalną układów drgających.

Równania drgań wymuszonych harmonicznie

Zakłada się, na konstrukcję działa obciążenie harmoniczne

0

( )

sin( )

t

t

θ

=

P

P

wywołujące

zmienne w czasie przemieszczenia, a także siły bezwładności. Równanie macierzowe
opisujące drgania wymuszone harmonicznie:

2

0

( -

)

θ

=

K

M Q

P

0

.

Dla danego wektora amplitud obciążenia wymuszającego

0

P należy obliczyć wektor

amplitud przemieszczeń

.

0

Q

Chcąc wyznaczyć amplitudy przemieszczeń

w czasie drgań wymuszonych harmonicznie,

przy danej macierzy podatności D i macierzy mas M, należy rozwiązać równanie
macierzowe:

0

Q

2

0

( -

)

θ

=

0

I

DM Q

DP .

W analizie drgań wymuszonych używane jest także nastające równanie macierzowe:

-1

0

0

2

1

( -

)

0

θ

=

D

M

B + DP

,

zawierające macierz podatności D, macierz mas M oraz – jako poszukiwane wielkości –
amplitudy sił bezwładności

. Równanie powyższe ma także skróconą postać:

0

B

0

0

0

p

+

=

*

D B

Δ

dzięki następującym podstawieniom:

-1

0

0

2

1

,

p

θ

= −

=

*

D

D

M

D

Δ

P .

8. Założenia i podstawy analizy statycznej prętów cienkościennych

Pręt cienkościenny to taki pręt, którego jeden z wymiarów określających przekrój poprzeczny
(grubość) jest nieporównywalnie mały w stosunku do drugiego. Pręty cienkościenne
wymagają innej analizy statycznej niż pręty lite, gdyż w ich przypadku nie obowiązują pewne
założenia upraszczające dla prętów litych: hipoteza Bernoulliego (przekrój płaski przed
przyłożeniem obciążenia pozostaje płaski po obciążeniu i prostopadły do ugiętej osi) i zasada
de Saint Venanta (przyłożone obciążenie można zastąpić innym statycznie równoważnym –
naprężenia, przemieszczenia i odkształcenia nie zmienią się poza niewielkim obszarem
zawierającym obciążone pole). Pręty cienkościenne można podzielić na pręty o profilu
otwartym i zamkniętym.

Analiza statyczna pręta cienkościennego zakłada wprowadzenie tzw. wycinkowych wielkości
geometrycznych na jego przekroju poprzecznym na bazie tzw. teorii grafu. Na profilu pręta
cienkościennego wybiera się tzw. korzeń, od którego odchodzą dendryty (drzewa) do każdego
innego węzła. Definiuje się lokalny krzywoliniowy układ współrzędnych ortogonalnych,

30

którego współrzędną skierowaną wzdłuż krawędzi profilu jest tzw. współrzędna łukowa.
Licząc całkę od bieguna (dowolny punkt płaszczyzny przekroju) aż do dowolnego punktu na
drodze skierowanej otrzyma się tzw. współrzędną wycinkową, stanowiącą podstawę
wyprowadzeń wszystkich innych wielkości geometrycznych profilu (tzw. charakterystyki
wycinkowe i geometryczne): momentu statycznego, momentów odśrodkowych, momentu
bezwładności i innych. Ważnymi pojęciami są także środek zginania (takie położenie
bieguna, dla którego obydwa wycinkowe momenty odśrodkowe zerują się) oraz główny punkt
zerowej współrzędnej wynikowej (główny korzeń).

Dla prętów cienkościennych otwartych poddanych prostemu skręcaniu wprowadzono
następujące założenia upraszczające:

• Przekrój poprzeczny aproksymuje się skończoną ilością prostokątów,

• Prostokąty pracują niezależnie – tak jednak, że jednostkowy kąt skręcenia każdego

nich jest taki sam.

Innymi pojęciami charakterystycznymi dla teorii prętów cienkościennych są bimoment i
bipara. Na ich podstawie powstaje odpowiednik zasady de Saint Venanta w prętach
cienkościennych (o kinematycznej równoważności układów).

U podstaw teorii pręta cienkościennego, oprócz założeń poprzednich, leżą następujące
hipotezy sformułowane przez Własowa:

• Powierzchnia środkowa (powierzchnia równoodległa od powierzchni górnej i dolnej)

deformuje się tak, jakby w płaszczyźnie każdego przekroju rozpostarta była na linii
środkowej sztywna tarcza, idealnie jednak wiotka dla deformacji w kierunku
prostopadłym do tego przekroju,

• Odkształcenia kątowe w punktach powierzchni środkowej są małe (

0

xs

γ

≈ ),

• Naprężenia normalne

n

σ

do powierzchni równoległych do powierzchni środkowej są

pomijalnie małe w stosunku do dwóch pozostałych naprężeń normalnych

,

x

s

σ σ

.

Do wyprowadzenia równania różniczkowego rządzącego (np. równania ruchu) stosuje się
najczęściej w teorii pręta cienkościennego podejście kinematyczne (przypuszczenie funkcji
przemieszczeń).

9. Siły przekrojowe w ustrojach prętowych

Rozważany będzie dowolny pręt określony w układzie współrzędnych globalnych
(

1

2

,

,

3

X X X ). Do pręta przyłożony jest dowolny układ sił zewnętrznych, również określony w

tym układzie poprzez współrzędne sił

1

2

3

( , , )

i

i

i

P

P P P

=

i

)

oraz punkty ich przyłożenia

1

2

3

(

,

,

i

i

i

i

A

X X X

=

1, 2,...

i

=

,

. W dowolnym punkcie Q dokonuje się przekroju poprzecznego

pręta: pręt rozpada się na dwie części: lokalny układ współrzędnych (

1

2

3

, ,

x x x ) powiązany jest

z II-gą częścią. Układ sił z części II zredukowany do bieguna Q składa się z dwóch wektorów
zaczepionych w Q:

• Sumy sił układu S :

i

i II

S

P

∈

=

∑

,

• Momentu układu względem bieguna Q:

,

i

i

i

i II

i

M

r P

r

QA

∈

=

×

=

∑

.

31