DYNAMIKA BUDOWLI
wszystko płynie
wszystko drga
Wszystkie procesy twórcze w przyrodzie dzieją się w stanach dalekich
od równowagi.
Michał Heller
Szczęście w przestrzeniach Banacha, 1997
Gdańsk, 2007
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
1. Drgania swobodne nietłumione
Drgania konstrukcji o jednym stopniu swobody:
1. Drgania swobodne:
a) Nietłumione: mu + ku = 0
b) Tłumione: mu + cu + ku = 0
2. Drgania wymuszone: mu + cu + ku = p(t)
a) Siłą harmoniczną p(t) = p0 sin(t)
b) Impulsem
c) Siłą dowolną
Drgania swobodne nietłumione
Siły sprężystości fs i bezwładności fi :
fs (t) = ku(t), (1.1)
fi (t) = mu(t). (1.2)
3 2
sila bezwladosci fi sila bezwladosci fi
sila sprezystosci fs sila sprezystosci fs
2
1
1
0 0
-1
-1
-2
-3 -2
0 0.5 1 1.5 2 -0.04 -0.02 0 0.02 0.04
u [m]
t [s]
________________________________________________________________________________________________
2007-02-26 2
sila [N]
sila [N]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Drgania swobodne nietłumione opisane są następującym równaniem:
mu + ku = 0 , (1.3)
gdzie m oznacza masę, k sztywność. Drgania swobodne zapoczątkowane są poprzez wytrącenie
układu z pozycji równowagi poprzez warunki początkowe, tzn. przyłożenie wychylenia początko-
wegou0 lub/i prędkości początkowej u0 w czasie t = 0 . Warunki początkowe mają następującą po-
stać:
u(0) = u0, u(0) = u0 . (1.4)
Równanie (1.3) wraz z warunkami początkowymi (1.4) tworzy zagadnienie własne. Rozwiązanie
powyższego równania stanowią w dynamice konstrukcji drgania własne nierzeczywistego układu
bez tłumienia.
Równanie charakterystyczne ma postać:
ms2 + k = 0. (1.5)
Jego rozwiązaniem są dwa pierwiastki zespolone:
s1 = in s2 =-in , (1.6)
gdzie n oznacza częstość kołową drgań (naturalną) mierzoną w rad/s wyrażoną wzorem:
k
n = . (1.7)
m
Rozwiązaniem ogólnym równania (1.3) jest:
1 2 nn
u(t) = Aes t + Aes t = Aei t + Ae-i t . (1.8)
1 2 1 2
Korzystając z zależności między funkcjami wykładniczymi a trygonometrycznymi:
n
ei t = cosnt + i sinnt,
(1.9)
n
e-i t = cosnt - i sinnt,
równanie (1.8) może być zapisane jako:
u(t) = Acosnt + Bsinnt . (1.10)
Stałe A oraz B zostaną wyznaczone z warunków brzegowych (1.4):
u(0)
A = u(0) B = . (1.11)
n
________________________________________________________________________________________________
2007-02-26 3
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
0.5
Tn
.
u(0)
C
u(0)
0
f
-0.5
0 2 4 6 8 10 12 14 16 18
t [s]
Amplituda drgań zależna od warunków początkowych:
C = A2 + B2 (1.12)
pozwala zapisać wzór (1.10) w postaci zwiniętej:
u(t) = C sin(nt +) , (1.13)
gdzie:
A
= arctg . (1.14)
B
Czas potrzebny do wykonania jednego cyklu drgań nazywany jest okresem drgań:
2Ą
Tn = . (1.15)
n
Częstotliwość mierzona w Hertzach wyraża się następującym wzorem:
1 n
fn = = . (1.16)
Tn 2Ą
Naturalna częstość kołowa drgań może być wyrażona w alternatywnej formie:
g
n = , (1.17)
st
gdzie ugięcie statyczne st wyraża się wzorem:
mg
st = . (1.18)
k
Stąd wzór (1.17) zapisać można jako:
k 1
n = = . (1.19)
m mst
Drgania swobodne nietłumione zaprezentowane są w pliku cw1_01.m
________________________________________________________________________________________________
2007-02-26 4
u [m]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 1.1
Dany jest wspornik o długości L, z masą m skupioną na jego końcu, EI = const. Obliczyć częstość
drgań swobodnych nietłumionych.
L
2
M1 L3
11 = ds =
+"
EI 3EI
0
1 3EI
n ==
m11 mL3
________________________________________________________________________________________________
2007-02-26 5
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 1.2
Dany jest wspornik o długości L = 2 m wykonany ze stali o module sprężystości E = 200GPa i gę-
stości = 7850kg/m3. Belka ma przekrój dwuteowy o następujących wymiarach: b = 10 cm,
h = 12 cm, g = 2 cm, t = 1 cm. Obliczyć częstotliwości drgań w obu kierunkach. Z jakiej długości
wspornika należy skupić masę, aby otrzymać poprawne wyniki?
Charakterystyki geometryczne przekroju poprzecznego:
Pole powierzchni A = 52 cm2
Ix = 2117.333 cm4
Momenty bezwładności:
Iy = 334.333 cm4
Masa skupiona z ź długości belki M = 20.41 kg
3EIx
y
n == 278.93 rad/s ! fny = 44.39 Hz
mL3
3EIy
x
n ==110.84 rad/s ! fnx =17.64 Hz
mL3
Porównując z analitycznym wzorem na pierwszą częstość drgań wspornika o masie rozłożonej ź :
3.515 EI
1 =
L2 ź
masa rozłożona po długości belki ź = 40.82 kg/m
3.515 EIx
y
n == 283.03 rad/s ! fny = 45.05 Hz
L2 ź
EIy
3.515
x
n ==112.47 rad/s ! fnx =17.90 Hz
L2 ź
Wniosek: skupienie masy z ź długości wspornika i wymodelowanie przy użyciu jednego stopnia
swobody dobrze odzwierciedla pierwszą częstość drgań.
________________________________________________________________________________________________
2007-02-26 6
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 1.3
Pierwsza częstość drgań belki swobodnie podpartej przy uwzględnieniu rzeczywistego ciągłego
rozkładu masy dana jest wzorem:
2
Ą EI
= .
L2 ź
Z jakiej długości belki swobodnie podpartej należy skupić masę, aby otrzymać poprawną częstość
drgań swobodnych. Przyjąć belkę długości L = 1200 mm, o wysokości przekroju h = 20 i szeroko-
ści b = 60 mm wykonaną z pleksiglasu; gęstość pleksiglasu wynosi r = 1190 kg/m3 a moduł sprę-
żystości E = 3300 MPa.
Ix = 4"10-8 m4
ź = 1.428 kg/m
Rzeczywista częstość drgań wynosi:
2
Ą EI
== 65.90 rad/s
L2 ź
Masa skupiona do środka belki z L:
M = 0.8568 kg
L3
11 =
48EI
1 48EI
n === 65.42 rad/s
M11 ML3
Masa skupiona do środka belki z L:
M =1.7136 kg
L3
11 =
48EI
1 48EI
n === 46.26 rad/s
M11 ML3
Wniosek: w przypadku modelowania belki swobodnie podpartej za pomocą jednego stopnia swobo-
dy masę należy zbierać z długości belki.
________________________________________________________________________________________________
2007-02-26 7
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 1.4
Dana jest rama jak na rysunku. Obliczyć częstość drgań swobodnych, jeżeli EI = const.
Jest to układ statycznie niewyznaczalny, rozważyć musimy dwa stany: obciążenie obciążeniem ze-
wnętrznym, oraz nadliczbową reakcją:
10 + 11X1 = 0
L
MM0 -L3
1
10 = dx =
+"
EI 32EI
0
L
MM1 L3
1
11 = dx =
+"
EI 2EI
0
1
X1 =
16
________________________________________________________________________________________________
2007-02-26 8
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Aby obliczyć przemieszczenie od siły P = 1 w układzie statycznie niewyznaczalnym możemy sko-
rzystać z twierdzenia redukcyjnego:
L
MM0 13 L3
11 = dx =
+"
EI 1536 EI
0
1 1536 EI
Częstość drgań naturalnych: n ==
m11 13 mL3
________________________________________________________________________________________________
2007-02-26 9
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 1.5
Wyznaczyć boczną sztywność ramy oraz odpowiadającą jej częstość drgań własnych.
a) EIb = "; b) EIb = 0 ; c) EIb = const.
EIc EIc 96 EIc
Odp. a) k = 24 ; b) k = 6 ; c= k =
h3 h3 7 h3
________________________________________________________________________________________________
2007-02-26 10
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
2. Więzi sprężyście odkształcalne
P = ku
P
k =
u
1 u
= =
k P
a) Połączenie szeregowe
Sumaryczne przemieszczenie układu sprężyn:
u = u1 + u2 +& + un . (2.1)
W każdej sprężynie działa stała siła:
PP P
u1 = , u2 = , & , un = . (2.2)
k1 k2 kn
Po podstawieniu otrzymujemy:
#ś#
1 1 1 P
u = P + +& + = , (2.3)
ś#
k1 k2 kn ź# k
# #
gdzie k oznacza sztywność zastępczą układu szeregowego:
n
1 1 1 1 1
= + +& + = . (2.4)
"
k k1 k2 kn i=1 ki
1 1 1
= +
k k1 k2
________________________________________________________________________________________________
2007-02-26 11
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
b) Połączenie równoległe
Siła działająca na układ jest sumą sił występujących we wszystkich sprężynach:
P = P1 + P2 +& + Pn . (2.5)
Przemieszczenie jest jednakowe dla wszystkich sprężyn:
P1 = k1u P2 = k2u & Pn = knu . (2.6)
Po podstawieniu otrzymujemy:
P = u(k1 + k2 +& + kn ) = ku , (2.7)
gdzie k oznacza sztywność zastępczą połączenia równoległego:
n
k = k1 + k2 +& + kn = . (2.8)
"k
i
i=1
k = k1 + k2
________________________________________________________________________________________________
2007-02-26 12
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 2.1
9EI
Wyznaczyć sztywność układu i częstość drgań swobodnych, jeżeli EI = const, kB = .
L3
Schemat rozpatrzymy jako superpozycję dwóch stanów:
a) ugięcie punktu 1 przy założeniu, że w punkcie B istnieje podpora stała
L
M12 L3
1 = dx =
+"
EI 8EI
0
b) ugięcie punktu 1 przy założeniu, że nieskończenie sztywna belka opiera się na podporze spręży-
stej
1 L3
B = RB =
kB 6EI
3 L3
2 = B =
2 4EI
3L3
= 1 +2 =
8EI
1 8EI
k = =
3L3
3L3
=
8mEI
________________________________________________________________________________________________
2007-02-26 13
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 2.2
Wyznaczyć sztywność i częstość drgań układu jak na rysunku:
96EI
Sztywność sprężyny łączącej belki wynosi ks = , EI = const.
L3
Układ składa się z trzech części o ustalonej sztywności:
- belka górna
L3 48EI
g = kg =
48EI L3
- belka dolna
L3 192EI
d = kd =
192EI L3
________________________________________________________________________________________________
2007-02-26 14
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
96EI
- sprężyna ks =
L3
Belka dolna i sprężyna łączą się szeregowo w układ, który połączony jest z belką górną równolegle:
kdks 64EI
k ' ==
kd + ks L3
112EI
Zastępcza sztywność układu: k = kg + k ' =
L3
112EI
=
mL3
________________________________________________________________________________________________
2007-02-26 15
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 2.3
Obliczyć częstość drgań własnych maszyny o ciężarze Q = 3 kN znajdującej się na końcu belki
wspornikowej długości 3 m, o przekroju jak na rysunku (b = 18 cm, h = 20 cm, grubość ścianki
d = 1 cm, E = 205 GPa), jeżeli pomiędzy belką a maszyną znajduje się podkładka o współczynniku
sprężystości k = 1 MN/m.
Ix = 2574.713 cm4
1 3EI
sztywność wspornika: kw = = = 586462.406 N/m
11 L3
sztywność podkładki: k = 1 MN/m
1 1 1
sztywność układu: = + = 0.000002705 ! kc =369666.753 N/m
kc kw k
k kg
= = = 34.77 rad/s
m Q
________________________________________________________________________________________________
2007-02-26 16
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 2.4
Znalezć częstość drgań własnych, okres drgań oraz ilość drgań na minutę fundamentu o ciężarze
Q = 2000kN. Pole powierzchni podstawy wynosi A = 10 m2, a współczynnik sprężystości podłoża
kz = 25 MN/m3.
________________________________________________________________________________________________
2007-02-26 17
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
3. Drgania swobodne tłumione
Siły sprężystości fs , bezwładności fi dane są wzorami (1.1) i (1.2) natomiast siła tłumienia fd
zdefiniowana jest następująco:
fd (t) = cu(t) . (3.1)
2 2
sila bezwladosci fi sila bezwladosci fi
1.5 1.5
sila sprezystosci fs
sila sprezystosci fs
sila tlumienia fd
sila tlumienia fd
1 1
0.5 0.5
0 0
-0.5 -0.5
-1 -1
-1.5 -1.5
-2 -2
0 5 10 15 20 -0.4 -0.2 0 0.2 0.4 0.6
u [m]
t [s]
Drgania swobodne tłumione opisane są następującym równaniem:
mu + cu + ku = 0 (3.2)
gdzie m oznacza masę, k sztywność, c tłumienie. Warunki początkowe, czyli wychylenie początko-
we u0 i/lub prędkość początkowa u0 , mają następującą postać:
u(0) = u0 u(0) = u0 . (3.3)
Równanie charakterystyczne jest postaci:
ms2 + cs + k = 0 . (3.4)
Jeżeli zdefiniujemy liczbę tłumienia jako:
c
= , (3.5)
ckr
________________________________________________________________________________________________
2007-02-26 18
i
s
d
i
s
d
sily f f f [N]
sily f f f [N]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
gdzie ckr oznacza tłumienie krytyczne:
ckr = 2mn , (3.6)
równanie (3.4) można zapisać w formie:
2
s2 + 2ns +n = 0 . (3.7)
Jego rozwiązaniem są dwa pierwiastki:
2
s1,2 =-n ą n -1 . (3.8)
Pierwiastki (3.8) mogą być rzeczywiste lub urojone, w zależności od wartości liczby . Rozważy-
my 3 przypadki:
>1 - tłumienie nadkrytyczne
Pierwiastki (3.8) równania (3.4) są rzeczywiste, a odpowiedz układu opisuje następujące równanie:
1 2
u(t) = C1es t + C2es t (3.9)
czyli:
22
n nn
u(t) = e- t (C1e- t -1 + C2e t -1) . (3.10)
Z warunków brzegowych mamy:
u(t = 0) = u0 = C1 + C2 ! C1 = u0 - C2 (3.11)
22
n nn
u(t) = e- t (-n)(C1e- t -1 + C2e t -1) +
(3.12)
22
22
n nn
+ e- t (-n -1C1e- t -1 + n -1C2e t -1)
22
u(t = 0) = u0 = C1(-n -n -1) + C2(-n + n -1) (3.13)
22 2
u0 = u0(-n -n -1) + C2(n +n -1 -n + n -1) (3.14)
2
u0 - u0(-n -n -1)
C2 = (3.15)
2
2n -1
2
u0(n -1 -n ) - u0
C1 = (3.16)
2
2n -1
________________________________________________________________________________________________
2007-02-26 19
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Po uwzględnieniu warunków brzegowych:
2 2
2 2
u0(n -1 -n) - u0 t -1 u0 - u0(-n -n -1)
n n n
u(t) = e- t ( e- + e t -1) =
22
2n -1 2n -1
(3.17)
22 22
nn nn
e- t -1 + e t -1 u0n + u0 e- t -1 - e t -1
n
= e- t[u0( ) + )]
(
2
22
n -1
Ostatecznie:
#ś#
u(0) + u(0)n
22
n
u(t) = e- t ś#u(0) cosh nt -1 + sinh nt -1 ź# (3.18)
( ) ( )
2
ś#ź#
n -1
# #
=1 - tłumienie krytyczne
Pierwiastki są rzeczywiste:
s1 = s2 =-n . (3.19)
Odpowiedz układu:
n
u(t) = e- t C1 + C2t . (3.20)
( )
Z warunków brzegowych:
u(t = 0) = u0 = C1 (3.21)
nn
u(t) = e- t (-n)(C1 + C2t) + C2e- t(3.22)
u(t = 0) = u0 = -nC1 + C2 ! C2 = u0 + nu0 (3.23)
Po uwzględnieniu warunków brzegowych:
n
u(t) = e- t[u0 + (u0 +nu0)t] (3.24)
n
u(t) = e- t u(0)(1+ nt) + u(0)t (3.25)
( )
<1 - tłumienie podkrytyczne
Pierwiastki są zespolone:
2
s1,2 =-n ą in 1- . (3.26)
Częstość drgań kołowych tłumionych zdefiniujmy jako:
2
D = n 1- . (3.27)
Odpowiedz układu:
1 2
u(t) = C1es t + C2es t , (3.28)
________________________________________________________________________________________________
2007-02-26 20
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
czyli:
n DD
u(t) = e- t C1ei t + C2e-i t . (3.29)
( )
Zamieniając na postać trygonometryczną:
n
u(t) = e- t AcosDt + B sinDt (3.30)
( )
Uwzględniając warunki brzegowe:
u(t = 0) = u0 = A (3.31)
nn
u =-ne- t AcosDt + B sinDt + e- t AsinDt + DB cosDt (3.32)
() (-D
)
u0 +nu0
u(t = 0) = u0 = -n A +DB ! B = (3.33)
D
#ś#
u0 + nu0
n
u(t) = e- t ś#u0 cosDt + sinDtź# (3.34)
D
# #
gdzie stałe A i B są następujące:
u(0) +nu(0)
A = u(0) B = (3.35)
D
0.6
tlumienie krytyczne
=1
0.4
tlumienie nadkrytyczne
>1
0.2
0
-0.2
tlumienie podkrytyczne
<1
czas [s]
-0.4
0 2 4 6 8 10 12 14 16 18 20
Drgania swobodne tłumione zaprezentowane są w pliku cw2_01.m
________________________________________________________________________________________________
2007-02-26 21
przemieszczenie [m]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Logarytmiczny dekrement tłumienia
Liczbę tłumienia można wyznaczyć eksperymentalnie poprzez pomiar drgań swobodnych.
W tym celu określić należy iloraz dwóch kolejnych maksymalnych wychyleń:
#ś#
u
u(t) 2Ą
j
= = expś#ź# (3.36)
2
ś#ź#
u(t + TD ) u
1-
j+1
# #
Logarytmicznym dekrementem tłumienia określa się wielkość będącą logarytmem naturalnym z
ilorazu dwóch kolejnych wychyleń:
uj 2Ą
= ln = (3.37)
uj+1 1- 2
Jeżeli tłumienie jest małe ( < 0.2 ), co ma miejsce w rzeczywistych konstrukcjach, logarytmicz-
ny dekrement tłumienia może być uproszczony do postaci:
2Ą (3.38)
Aby określić liczbę tłumienia, niekoniecznie trzeba posługiwać się ilorazem dwóch kolejnych am-
plitud. Po j cyklach amplituda zmniejsza wartość z u1 do uj+1 :
u1 u1 u2 u3 uj j
= = e (3.39)
uj+1 u2 u3 u4 uj+1
a logarytmiczny dekrement tłumienia przyjmuje postać:
1 u1
= ln 2Ą (3.40)
j u
j+1
Stąd łatwo wyznaczyć liczbę cykli potrzebną do zmniejszenia amplitudy o zadaną wartość. Na
przykład liczba cykli j potrzebna do zmniejszenia amplitudy o 50 % wynosi:
0.11
j50% (3.41)
10
2Ą
=
8
2
1-
6
4
= 2Ą
2
0
0 0.2 0.4 0.6 0.8 1
Liczba tłumienia
________________________________________________________________________________________________
2007-02-26 22
Logarytmiczny dekrement tłumienia
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 3.1
Dany jest układ o znanej liczbie tłumienia . Znalezć liczbę cykli drgań swobodnych potrzebną do
redukcji amplitudy drgań do 10% w stosunku do amplitudy początkowej. Prędkość początkowa jest
równa zeru.
0.366
j10%
ZADANIE 3.2
Dany jest zbiornik na wodę. Za pomocą kabla przymocowanego do
górnej części zbiornika przyłożono boczną siłę Q = 10 MN
powodując wychylenie zbiornika o 2 m. Kabel przecięto
wprowadzając konstrukcję w drgania. Po upływie 2 sekund, zbiornik
wykonał 4 pełne cykle drgań, a amplituda zmniejszyła się do 1 m.
Obliczyć:
a) liczbę tłumienia
b) okres drgań nietłumionych
c) sztywność
d) masę
e) tłumienie
f) liczbę cykli potrzebną do zmniejszenia amplitudy drgań do 0.2m.
a) = 2.76%
b) Tn = 0.4998 s
c) k = 5 MN/m
d) m = 31.63955 t
e) c = 21955.29 kg/s
f) j 13.28 cykli
________________________________________________________________________________________________
2007-02-26 23
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 3.3
Znalezć częstość naturalną drgań i liczbę tłumienia belki wspornikowej na podstawie eksperymen-
talnie pomierzonego przyspieszenia swobodnego końca belki. Porównać wyznaczoną naturalną
częstość drgań z częstością wyznaczoną analitycznie poprzez skupienie masy belki do jednego
punktu. Wykonana z pleksiglasu belka ma długość L = 480 mm, a jej przekrój jest prostokątny wy-
sokości h = 20 i szerokości b = 60 mm. Gęstość pleksiglasu wynosi r = 1190kg/m3 a moduł Youn-
ga E = 3300 MPa.
80
..
..
u11=1.4088 m/s2
u5=5.3293 m/s2
60
t11=0.6376 s
t5=0.3794 s
40
20
0
-20
-40
-60
-80
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
czas [s]
1 ui
= ln
2Ą j u
j+i
1 u5 1 5.3293
= ln = ln = 0.035 = 3.5%
2Ą 6 u11 12Ą 1.4088
Ponieważ w czasie t = 0.6376 s - 0.3794 s = 0.2582 s belka wykonuje 6 cykli drgań to okres drgań
0.2582 s
tłumionych TD jest równy TD == 0.0403 s
6
2
Okres drgań naturalnych Tn = TD 1- = 0.0403 1- (0.035)2 = 0.0402 s
2Ą 2Ą n
Naturalna częstość drgań n = = =156.298 rad/s ! fn = = 24.875 Hz
Tn 0.0402 2Ą
Teraz obliczymy częstość drgań skupiając masę z ź długości wspornika: M = 0.17136 kg
Ix = 4"10-8 m4
3EIx n
n == 144.554 rad/s ! fn = = 23.006 Hz
ML3 2Ą
________________________________________________________________________________________________
2007-02-26 24
2
p r z y s pi es z eni e [ m/ s ]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 3.4
Maszyna o ciężarze Q = 250 kN zamocowana jest do podłoża za pomocą 4 sprężyn i 4 tłumików.
Pionowe przemieszczenie pod wpływem ciężaru własnego maszyny wynosi 0.8 cm. Tłumiki zapro-
jektowano tak, aby redukowały amplitudę pionowych drgań do 1/8 początkowej amplitudy po 2
cyklach drgań. Obliczyć:
a) częstość drgań nietłumionych;
b) liczbę tłumienia;
c) częstość drgań tłumionych.
Odp. n = 35.02 rad/s , = 0.165 , d = 34.54 rad/s
________________________________________________________________________________________________
2007-02-26 25
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
4. Drgania wymuszone siłą harmoniczną
Drgania wymuszone siłą harmoniczną nietłumione
Drgania wymuszone sinusoidalną siłą harmoniczną nietłumione opisane są następującym równa-
niem:
mu + ku = p0 sint , (4.1)
gdzie m oznacza masę, k - sztywność, - częstość siły wymuszającej, p0 amplitudę siły wymu-
szającej. Warunki początkowe, czyli wychylenie początkowe u0 i/lub prędkość początkowa u0 ,
mają następującą postać:
u(0) = u0 u(0) = u0 . (4.2)
Równanie (4.1) jest niejednorodne. Jego rozwiązanie jest sumą całki ogólnej i całki szczególnej:
u(t) = u(t)c + u(t)p . (4.3)
Całka ogólna równania jednorodnego:
mu + ku = 0 . (4.4)
Równanie charakterystyczne ma postać:
s2 + n2 = 0 , (4.5)
s1,2 = ąin . (4.6)
Całka ogólna ma postać:
uc (t) = Acosnt + Bsinnt . (4.7)
Całka szczególna równania niejednorodnego (4.1) ma postać:
up (t) = Dsint + E cost . (4.8)
Wyznaczenie stałych C i D:
up = D cost - E sint , (4.9)
up =-D2 sint - E2 cost . (4.10)
Podstawiając(4.9) i (4.10) do równania (4.1) zapisanego w postaci:
po
u +n2u = sint , (4.11)
m
________________________________________________________________________________________________
2007-02-26 26
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
po
-D2 sint - E2 cost + n2Dsint +n2E cost = sint (4.12)
m
po
cost En2 - E2 + sint Dn2 - D2 = sint (4.13)
() ()
m
ż# - E2 = 0 ! E = 0
En2
#
(4.14)
#ś#
# po po 1
2
ś#
n
#D - D2 = ! D =
mm n2 -2 ź#
# #
#
po 1
D = (4.15)
2
m
n2 /n
( )
(1- )
po 1
D = (4.16)
2
k
1-( )
/n
Całka ogólna równania niejednorodnego (4.1):
u(t) = Acosnt + Bsinnt + Dsint . (4.17)
Wyznaczenie stałych A i B z warunków brzegowych (4.2):
u(0) = A ! A = uo (4.18)
u(t) =-An sinnt + Bn cosnt + D cost (4.19)
uo - D
u(0) = Bn + D ! B = (4.20)
n
uo po /n
B = - (4.21)
2
n k
( )
(1- /n )
Po podstawieniu stałych A, B, D do równania (4.17) otrzymujemy ostateczną odpowiedz układu:
Ą#ń# ń#
Ą#
u(0) p0 /n p0 1
u(t) = u(0)cosnt + - sinnt + sint (4.22)
ó# ó#
n k 1- ( /n )2 Ą# k 1- ( /n)2 Ą#
Ł#Ś# Ś#
Ł#
drgania "zanikające" drgania ustalone
0.03
odpowiedz calkowita
drgania ustalone
0.02
0.01
0
-0.01
-0.02
-0.03
0 1 2 3 4 5 6 7
czas [s]
________________________________________________________________________________________________
2007-02-26 27
przemieszczenie[m]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Współczynnik dynamiczny, zjawisko rezonansu:
Rozważmy drgania ustalone o częstości siły wymuszającej :
p0 1
u(t) = sint (4.23)
k 1- ( /n)2
Wielkość p0 / k może być interpretowana jako ugięcie statyczne układu wywołane przez siłę p0
przyłożoną w sposób statyczny:
p0
(ust )0 = (4.24)
k
a równanie (4.23) przyjmie postać:
1
u(t) = (ust )0 sint (4.25)
1- ( /n)2
odpowiedz rezonansowa ukladu nietlumionego
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
0 1 2 3 4 5 6 7
t [s]
Współczynnikiem dynamicznym Rd (zwielokrotnienia amplitudy drgań) nazywamy iloraz amplitu-
dy drgań do amplitudy ugięcia statycznego:
1
Rd = (4.26)
1- ( /n)2
u(t) = (ust )0 Rd sint (4.27)
max u(t) = u0 = (ust )0 Rd (4.28)
Funkcja dynamiczności obciążenia R(t) - iloraz odpowiedzi konstrukcji u(t) spowodowanej siłą
zmienną w czasie do przemieszczenia statycznego, tzn. przemieszczenia powstałego pod działaniem
siły statycznej równej maksymalnej wartości siły zmiennej w czasie:
u(t)
R(t) = = Rd sint (4.29)
(ust )0
________________________________________________________________________________________________
2007-02-26 28
u [ m]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Częstotliwość rezonansowa zdefiniowana jest jako częstotliwość, dla której współczynnik Rd osią-
ga maksymalną wartość, czyli dla n .
60
60
50
40
40
20
30
0
20
-20
10
-40
-60 0
0 1 2 3 4 0 1 2 3 4
/ [-]
/ [-] n
n
Drgania wymuszone siłą harmoniczną tłumione
Drgania wymuszone sinusoidalną siłą harmoniczną tłumione opisane są następującym równaniem:
mu + cu + ku = p0 sint , (4.30)
gdzie m oznacza masę, k - sztywność, c - tłumienie, - częstość siły wymuszającej, p0 amplitudę
siły wymuszającej. Warunki początkowe, czyli wychylenie początkowe u0 i/lub prędkość począt-
kowa u0 , mają następującą postać:
u(0) = u0 u(0) = u0 . (4.31)
Równanie (4.30) jest niejednorodne. Jego rozwiązanie jest sumą całki ogólnej i całki szczególnej:
u(t) = u(t)c + u(t)p (4.32)
Całka ogólna równania (4.30) dla tłumienia podkrytycznego ma postać (3.30):
n
u(t) = e- t AcosDt + B sinDt , (4.33)
( )
natomiast całka szczególna jest postaci (4.8):
up (t) = Dsint + E cost . (4.34)
________________________________________________________________________________________________
2007-02-26 29
2
n
d
R [-]
1/ ( 1- ( /
) ) [ -]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Wyznaczenie stałych D i E:
up = D cost - E sint (4.35)
up =-D2 sint - E2 cost (4.36)
Podstawiając (4.35) i (4.36) do równania (4.30) zapisanego w postaci:
po
u + 2nu + n2u = sint (4.37)
m
otrzymujemy związek:
-D2 sint - E2 cost + 2n D cost - E sint +n2 Dsint + E cost =
( ) ( )
(4.38)
po
= sint
m
po
sint -D2 - 2nE +n2D + cost -E2 + 2nD +n2E = sint (4.39)
() ()
m
ż#
E 2 -n2
( )
#-E2 + 2nD + n2E = 0 ! D =
#
2n
(4.40)
#
# po
2
#-D - 2nE + n2D =
# m
E 2 -n2 E 2 -n2 po
() ( )
-2 - 2nE +n2 = (4.41)
2n 2n m
2
#ś#
- n2 -2
()
po
ś#ź#
E - 2n = (4.42)
ś#ź#
2n m
# #
2
#
n2 -2 + (2n)2 ś# po
()
ś#ź#
E = (4.43)
ś#ź#
-2n m
# #
pon
-2
E = (4.44)
2
m
n2 -2 + (2n)2
()
p0 -2 ( /n )
E = (4.45)
2
k
Ą#Ś# 2
[]
Ł#1- ( /n)2 ń# + 2 ( /n)
Na podstawie zależności (4.40) otrzymujemy stałą D :
E 2 -n2 p0 2 -n2
( ) ( )
-2 ( /n)
D == (4.46)
2
2 k
Ą#Ś# 2 2
n
[]
Ł#1- ( /n )2 ń# + 2 ( /n) n
________________________________________________________________________________________________
2007-02-26 30
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
p0 1- ( /n)2
D = (4.47)
2
k
Ą#Ś# 2
[]
Ł#1- ( /n)2 ń# + 2 ( /n )
Wyznaczenie stałych A i B. Całka ogólna równania niejednorodnego (4.33):
n
u(t) = e- t AcosDt + B sinDt + Dsint + E cost (4.48)
()
u(t = 0) = A + E = uo ! A = uo - E (4.49)
p0 -2 ( /n)
A = uo + (4.50)
2
k
Ą#Ś# 2
[]
Ł#1- ( /n)2 ń# + 2 ( /n)
n
u(t) = e- t Ą#-n AcosDt + BsinDt AD sinDt + BD cosDtń# +
()-
Ł# Ś#
(4.51)
+ D cost - E sint
u(t = 0) = -n A + BD + D = uo (4.52)
uo - D + An uo n
B == - D + A (4.53)
D D D D
uo p0 1- ( /n)2
B = - +
2
D D k
Ą#Ś# 2
[]
Ł#1- ( /n )2 ń# + 2 ( /n )
(4.54)
ś#
n # p0
ś#uo + -2 ( /n ) ź#
+
2
ś#ź#
D k
Ą#Ś# 2
[]
Ł#1- ( /n )2 ń# + 2 ( /n )
# #
Ostatecznie odpowiedz układu jest wyrażona następująco:
n
u(t) = e- t AcosDt + BsinDt + Dsin E cos (4.55)
()
t t
+
drgania ustalone
drgania zanikające
0.02
odpowiedz calkowita
drgania ustalone
0.015
0.01
0.005
0
-0.005
-0.01
-0.015
0 1 2 3 4 5 6 7 8 9 10
czas [s]
________________________________________________________________________________________________
2007-02-26 31
pr z emi es z c z eni e[ m]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Dla częstości wymuszenia = n odpowiedz układu przyjmuje postać:
Ą#ń#
#ś#
p0 1
n
u(t) = ó#Ą# (4.56)
e- t ś#cosDt + sinDt ź# - cosnt
2
ś#ź#
k 2
ó#Ą#
1-
# #
Ł#Ś#
Odpowiedz układu nie przekracza wartości u0 :
(ust )0
u0 = (4.57)
2
odpowiedz rezonansowa ukladu tlumionego
ust 0
( )
0.05
2
0
ust 0
( )
-0.05
-
0 2 4 6 8 10 12 14
2
t [s]
Współczynnik dynamiczny:
1
Rd = (4.58)
2
Ą#Ś# 2
[]
Ł#1- ( /n)2 ń# + 2 ( /n)
10
=0.01
=0.1
5
=0.2
1
0
=1.0
=0.7
-5
0 0.5 1 1.5 2 2.5 3
/ [ - ]
n
________________________________________________________________________________________________
2007-02-26 32
u [ m]
d
R [ - ]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 4.1
Na belce o długości L = 4 m umieszczony jest silnik o ciężarze Q = 40 kN wywołujący drgania si-
nusoidalne o amplitudzie p0 = 2 kN i częstości wymuszenia 500 obrotów na minutę. Obliczyć
amplitudę drgań nietłumionych oraz tłumionych dla = 5% . Zbadać przypadek rezonansu dla
drgań tłumionych. EI = 24150 kNm2.
L
M12 L3
11 = ds =
+"
EI 48EI
0
QL3
Ugięcie statyczne od ciężaru silnika: Q = = 0.2208 cm
48EI
1 48EI
Sztywność układu: k = = =1811250 N/m
11 L3
1 48EI
Częstość drgań: n == = 66.65 rad/s
m11 mL3
Ugięcie statyczne układu wywołane przez siłę p0 przyłożoną w sposób statyczny:
p0
(ust )0 = = 0.01104 cm
k
Częstość wymuszenia: = 52.36 rad/s
1
Współczynnik dynamiczny bez uwzględnienia tłumienia: Rd == 2.61
1- ( /n)2
Współczynnik dynamiczny z uwzględnieniem tłumienia:
1
t
Rd == 2.557
2
Ą#Ś# 2
[]
Ł#1- ( /n)2 ń# + 2 ( /n)
Amplituda drgań nietłumionych: u0 = (ust )0 Rd = 0.0288 cm
t
Amplituda drgań tłumionych: u0 = (ust )0 Rd = 0.0282 cm
W przypadku rezonansu: = n = 66.65 rad/s
Współczynnik dynamiczny:
1
t
Rd == 10
2
Ą#Ś# 2
[]
Ł#1- ( /n )2 ń# + 2 ( /n)
t
Amplituda drgań tłumionych: u0 = (ust )0 Rd = 0.1104 cm
Odpowiedz układu nie przekracza wartości u0 :
(ust )0
u0 = = 0.1104 cm
2
________________________________________________________________________________________________
2007-02-26 33
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 4.2
Jaka powinna być stała sprężyny ks , którą należy wstawić pod silnik o wadze Q = 2 kN, aby
współczynnik dynamiczny Rd drgań harmonicznych o częstości = 10 rad/s spełniał warunek
Rd < 0.5? Belka ma długość L = 2 m, I = 328 cm4, E = 200 GPa.
Sztywność belki:
L3
11 =
3EI
3EI
kb =
L3
Sztywność szeregowego połączenia układu belka-sprężyna:
kbks
k =
kb + ks
3EI
ks
k kbks
2
L3
Częstość naturalna drgań: n = = =
3EI
m m(kb + ks )
m# + ks ś#
ś#
L3 ź#
# #
1 2
Rd =< 0.5 czyli > 3
2
2
n
1-( )
/n
3EI
2
m23EI
L3
ks <=
9EI
#
-2 ś# 9EI - m2L3
ś#
mL3 ź#
# #
ks < 6988.85 N/m
Przyjęto: ks =
n = rad/s
________________________________________________________________________________________________
2007-02-26 34
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 4.3
Masa m, sztywność k oraz naturalna częstość n nietłumionego układu o jednym stopniu swobody
są nieznane. Do wyznaczenia tych wielkości zastosowano test wzbudzenia harmonicznego. Przy
częstotliwości wzbudzenia 4 Hz wystąpił rezonans. Następnie masę zwiększono o 5 kg i wówczas
rezonans wystąpił przy częstotliwości wymuszenia 3 Hz. Wyznaczyć masę i sztywność.
m = 6.43 kg , k = 4061 N/m
ZADANIE 4.4
Układ o jednym stopniu swobody poddano wzbudzeniu siłą sinusoidalną. Przy rezonansie amplitu-
da przemieszczenia została pomierzona i wynosiła 2 m. Przy wzbudzeniu częstością równą 1/10
naturalnej częstości n amplitudą wyniosła 0.2 m. Ile wynosi liczba tłumienia?
= 0.0495
________________________________________________________________________________________________
2007-02-26 35
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
5. Projektowanie konstrukcji o jednym stopniu swobody obciążonych dy-
namicznie
Dany jest układ o jednym stopniu swobody (m, k, c) obciążony siłą p(t).
Równanie ruchu:
mu + cu + ku = p(t) (5.1)
Przenosząc wyrazy związane z masą i tłumieniem na prawą stronę otrzymujemy:
ku = p(t) mu cu (5.2)
- -
pz
gdzie pz oznacza zastępczą siłę statyczną.
pz = ku (5.3)
Siła zastępcza pz osiąga wartość minimalna i maksymalną dla minimalnego i maksymalnego prze-
mieszczenia u:
pz min = kumin
(5.4)
pz max = kumax
________________________________________________________________________________________________
2007-02-26 36
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Obwiednia momentów jest sumą momentów od obciążenia statycznego oraz dynamicznego:
Mdyn max
obw M = MQ ą (5.5)
Mdyn min
M1 pz max
obw M = M1Q ą (5.6)
M1 pz min
Przypadki szczególne:
Drgania swobodne bez tłumienia:
mu + ku = 0
u(t) = C sin(nt +)
2
pz =-mu = mnC = kC pz min = -kC pz max = kC
Drgania wymuszone harmonicznie:
mu + cu + ku = p0 sint
umin,max = "(ust )0 Rd
pz min =-kumin =- p0Rd
pz max = kumax = p0Rd
________________________________________________________________________________________________
2007-02-26 37
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 5.1
Wyznaczyć obwiednię momentów zginających odciężaru masy i dynamicznych momentów wywo-
łanych warunkiem początkowym u(t = 0) = 0.4 m/s . Belka ma długość L = 4.8 m, wykonana jest z
dwóch stalowych dwuteowników I 180 (Ix-x = 1450 cm4). Moduł sprężystości stali E = 200 GPa.
Ciężar Q = 20 kN.
Moment bezwładności przekroju Ix = 2900 cm4
M12
Przemieszczenie od siły jednostkowej: 11 = ds = 0.496"10-6 m/N
+"
EI
1
Częstość drgań własnych: n == 31.447 rad/s
m11
2
Ą# ń#
2
Amplituda drgań masy: C = A2 + B2 = u(0) + = 0.01272 m
[ ]u(0)
ó# Ą#
n
Ł# Ś#
C
Zastępcza siła statyczna pz = kC = = 25645.161 N
11
MQ = MQ
1
Mdyn =ąM1 pz
________________________________________________________________________________________________
2007-02-26 38
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 5.2
Na belce (I220, Ix-x = 3060 cm4) znajduje się maszyna o ciężarze Q = 20 kN wywołująca wymusza-
jącą siłę harmoniczną amplitudzie p0 = 2 kN . Liczba obrotów wirnika maszyny wynosi 400 obro-
tów/minutę. Dobrać podkładkę pod maszynę tak, by maksymalne naprężenia przy zginaniu nie
przekraczały dop = 90 MPa . L = 5m.
M12
Przemieszczenie od siły jednostkowej: 11 = ds = 0.373"10-6 m/N
+"
EI
Częstość naturalna drgań n = 36.263 rad/s
Częstość wymuszenia: = 41.888 rad/s
1
Współczynnik dynamiczny: Rd == 2.991
1- ( /n)2
Zastępcza siła statyczna: pz = Rd p0 = 5.982 kN
Moment maksymalny: Mmax = (Q + pz )M1 = 31.178 kNm
Mmax
Maksymalne naprężenia: max = y =112.077 MPa > 90 MPa = dop
Ix
Naprężenia maksymalne przekraczają wartość dopuszczalną, trzeba dobrać podkładkę zmieniającą
częstość drgań własnych układu. Dopuszczalna zastępcza siła statyczna musi wynosić:
pz d" 0.8636 kN
pz = Rd p0 = 0.8636 kN
1
Rd == 0.432
1- ( /n)2
________________________________________________________________________________________________
2007-02-26 39
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
2
# ś#
= 3.315
ś# ź#
n
# #
2
n = 529.292 rad/s
1 1 1
= +
k kb ks
1 1 1
= +
2 2
2
n n s
Szukana sztywność podkładki wynosi:
ks =1806 kN/m
________________________________________________________________________________________________
2007-02-26 40
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 5.3
Na belkę ciągłą przegubową działa obciążenie harmoniczne. Obliczyć maksymalne wychylenie
masy oraz amplitudy dynamicznych momentów zginających oraz narysować obwiednię momentów,
jeżeli L = 4 m, EI = 40000 kNm2, p0 = 6 kN, = 35 rad/s, m = 4.8 Mg.
M12 1
Ugięcie od siły jednostkowej: 11 = ds = m/kN
+"
EI 12000
1
Sztywność: k = = 12000 kN/m
11
Częstość drgań naturalnych: n = 50 rad/s
1
Współczynnik dynamiczny: Rd ==1.961
1- ( /n)2
Zastępcza siła statyczna pz = p0Rd = 11.8 kN
Ugięcie maksymalne: umax = uQ + ust 0 Rd = 11(Q + pz ) = 0.49 cm
( )
Momenty zginające wywołane statycznym działaniem ciężaru masy:
Amplitudy dynamicznych momentów zginających:
Obwiednia momentów:
________________________________________________________________________________________________
2007-02-26 41
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 5.4
Rygiel ramy obciążony jest w środku rozpiętości masą skupioną m = 3Mg i siłą harmoniczną o am-
plitudzie p0 = 6 kN i częstości wymuszenia = 24 rad/s. Narysować obwiednię momentów zgina-
jących, jeżeli L = 4 m, EI = 12500 kNm2.
L
MM0 20
Ugięcie od siły jednostkowej: 11 = dx =
+"
EI 3EI
0
1 3 EI
Częstość drgań naturalnych: n == = 25 rad/s
m11 20 m
________________________________________________________________________________________________
2007-02-26 42
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
1
Współczynnik dynamiczny: Rd ==12.755
1- ( /n )2
Zastępcza siła statyczna pz = Rd p0 = 76.531 kN
MQ = MQ = 29.43 M
Mdyn =ąMpz =ą76.531 M
obwM = (29.43 ą 76.531)M
________________________________________________________________________________________________
2007-02-26 43
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
6. Drgania wywołane dowolnym obciążeniem wymuszającym
6.1. Działanie impulsu jednostkowego
impuls jednostkowy gdy 0
p = mu
t2
pdt = m(u2 - u1) = m"u ! impuls
+"
t1
Drgania swobodne układu o jednym stopniu swobody (bez tłumienia) opisane są równaniem:
u(0)
u(t) = u(0)cosn (t - ) + sinn(t - ). (6.1)
n
Wstawiając warunki brzegowe:
u(0) = 0 i u(0) ="u = 1/ m , (6.2)
do równania (6.1) otrzymujemy odpowiedz układu:
1
h(t - ) a" u(t) = sin n(t - ) t e" . (6.3)
[]
mn
Dla układu tłumionego odpowiedz jest następująca:
1
n
h(t - ) a" u(t) = e- (t- ) sin d (t - ) t e" . (6.4)
[]
md
________________________________________________________________________________________________
2007-02-26 44
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
6.2. Impuls prostokątny
Równanie ruchu:
p0 t d" td
ż#
mu + ku = p(t) = (6.5)
#
0 t e" td
#
z warunkami brzegowymi u(0) = 0 i u(0) = 0 . Ruch układu przebiega w 2 fazach:
a) faza działania impulsu t d" td , podczas której układ poddany jest działaniu siły p(t) = p0 ( step
force ). Odpowiedz układu:
p0
u(t) = Acosnt + Bsinnt + up , up = (6.6)
k
- p0
Z warunków brzegowych : u(0) = 0 i u(0) = 0 otrzymujemy A = i B = 0
k
p0
u(t) = (1- cosnt), t d" td (6.7)
k
u(t) = ust 0 (1- cosnt), t d" td (6.8)
( )
b) faza drgań swobodnych t e" td
u(td )
u(t) = u(td )cosn t - td + sinn t - td (6.9)
( ) ( )
n
Drgania swobodne zapoczątkowane są prędkością i przemieszczeniem masy w czasie t = td , wy-
znaczonymi z równania (6.7):
u(td ) = ust 0 (1- cosntd ), u(td ) = ust 0 n sinntd . (6.10)
( ) ( )
Podstawiając (6.10) do równania (6.9) otrzymujemy:
ust 0 n sinntd
u(t) = ust 0 (1- cosntd )cosn t - td + sinn t - td , t e" td (6.11)
( ) ( )( ) ( )
n
u(t) = ust 0 Ą#(1- cosntd )cosn t - td + sinntd sinn t - td ń# , t e" td (6.12)
( ) ( ) ( )Ś#
Ł#
u(t) = ust 0 Ą#cosn t - td - cosntd ń# , t e" td (6.13)
( ) ( )
Ł#Ś#
________________________________________________________________________________________________
2007-02-26 45
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Odpowiedz układu nietłumionego o jednym stopniu swobody na impuls prostokątny:
12 12
Tn/td =2 Tn/td =1.75
10 10
8 8
6 6
4 4
2 2
0 0
-2 -2
-4 -4
0 1 2 3 4 0 1 2 3 4
t [s] t [s]
12 12
Tn/td =1.5 Tn/td =1.25
10 10
8 8
6 6
4 4
2 2
0 0
-2 -2
-4 -4
0 1 2 3 4 0 1 2 3 4
t [s] t [s]
12 12
Tn/td =1 Tn/td =0.5
10 10
8 8
6 6
4 4
2 2
0 0
-2 -2
-4 -4
0 1 2 3 4 0 1 2 3 4
t [s] t [s]
12 12
Tn/td =0.25 Tn/td =0.1
10 10
8 8
6 6
4 4
2 2
0 0
-2 -2
-4 -4
0 1 2 3 4 0 1 2 3 4
t [s] t [s]
________________________________________________________________________________________________
2007-02-26 46
u [ m]
u [ m]
u [ m]
u [ m]
u [ m]
u [ m]
u [ m]
u [ m]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
6.3. Przybliżona odpowiedz układu o jednym stopniu swobody dla krótkich impulsów:
Impuls krótki = impuls spełniający warunek td /Tn < 0.5
Jeżeli czas trwania impulsu td jest krótszy od Tn / 2 , wówczas maksymalne przemieszczenie nie
występuje w trakcie trwania impulsu, lecz w fazie drgań swobodnych, a impuls może być traktowa-
ny jako czysty impuls amplitudzie J:
td
J = p(t)dt
+"
0
Odpowiedz układu na impuls J o czasie trwania td spełniającym warunek td /Tn < 0.5 jest odpo-
wiedzią układu na impuls jednostkowy zgodnie z równaniem (6.3):
#ś#
1
u(t) = J sinnt . (6.14)
ś#ź#
mn
# #
ZADANIE 6.1
Na belce o rozpiętości L = 7 m znajduje się masa m = 1800 kg. Do masy przyłożono nagle impuls
prostokątny o amplitudzie P = 9 kN działający przez czas 0.01 s. Wyznaczyć amplitudę drgań belki
jeżeli E = 210 GPa, Ix = 4000 cm4.
48 1
11 =
7 EIx
7EIx
n == 26.087 rad/s
48m
2Ą
Tn = = 0.241
n
td /Tn = 0.041
impuls :
J = P "td = 90 Ns
amplituda:
J
umax = =1.917 "10-3 m
mn
________________________________________________________________________________________________
2007-02-26 47
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
6.4. Działanie dowolnej siły wymuszającej wybrane sposoby rozwiązania równania ruchu o
jednym stopniu swobody:
a) całka Duhamela
function [u]=duhamel(m,wn,wd,ksi,p,t)
20
15
10
5
0
0 0.5 1 1.5 2 2.5 3 3.5 4
t [s]
0.4
0.2
0
-0.2
-0.4
0 0.5 1 1.5 2 2.5 3 3.5 4
t [s]
Całka Duhamela dla układu nietłumionego:
t
p( )
u(t) = sin n (t - ) d . (6.15)
[]
+"
mn
0
Całka Duhamela dla układu tłumionego:
t
p( )
n
u(t) = e- (1- ) sin d (t - ) d . (6.16)
[]
+"
md
0
Całka Duhamela z uwzględnieniem warunków brzegowych dla układu nietłumionego:
t
#ś#
u(0) p( )
u(t) = cosnt + sindt + sin n (t - ) d . (6.17)
[]
ś#u(0) ź#
+"
n mn
# # 0
Całka Duhamela z uwzględnieniem warunków brzegowych dla układu tłumionego:
t
#ś#
u(0) + nu(0) p( )
-nt
n
u(t) = cosdt + sindt + e- (t- ) sin d (t - ) d .(6.18)
[]
ś#u(0) ź#e md
+"
d
# # 0
________________________________________________________________________________________________
2007-02-26 48
p(t) [N]
u(t) [m]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
b) równanie ruchu w przestrzeni stanów
Równanie ruchu układu tłumionego:
mu + cu + ku = p(t) (6.19)
gdzie: p(t) - dowolna siła wymuszająca
Proces dynamiczny może byś opisany zmiennymi xi definiującymi stan zjawiska. Zmienne xi nazy-
wamy zmiennymi przestrzeni stanu (lub prościej stanami). Do jednoznacznego opisu procesu dy-
namicznego o jednej zmiennej niezbędne są dwa stany. Na przykład przemieszczenie u i przyspie-
szenie u .
Załóżmy, że:
x1 = u - stan nr 1
x2 = u - stan nr 2
Przekształcając równanie (6.19) otrzymujemy:
c k 1
u + u + u = p(t)
m m m
c k 1
oraz u =- u - u + p(t)
m m m
zapisując równanie poprzez stany uzyskujemy:
c k 1
x2 =- x2 - x1 + p(t) (6.20)
m m m
dodatkowo zależność między stanami można zapisać jako:
x1 = x2 (6.21)
Równania (6.20) i (6.21) zapisane wektorowo przyjmują postać:
x=Ax+Bp(t) (state space equation of motion)
x1
Ą# ń#
gdzie: x= , wektor stanów,
ó#x Ą#
Ł# 2 Ś#
0 1
Ą#ń#
A=ó# - macierz układu
Ł#-k / m -c / mĄ#
Ś#
0
Ą# ń#
B=
ó#1/ mĄ#
Ł# Ś#
Drugim równaniem uogólniającym procesy dynamiczne jest równanie wyjścia:
y =Cx +Dp(t) (output equation)
________________________________________________________________________________________________
2007-02-26 49
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Jeżeli jako decydujące dla opisu dynamiki wybierzemy stan związany z przemieszczeniem, xi, to
macierze C i D przyjmą postać:
C= 1 0 D= 0
[ ] [ ]
Macierze A, B, C i D definiują układ dynamiczny w przestrzeni stanów.
Wyjaśnienie użycia funkcji lsim w programie Matlab.
Format funkcji:
[Y, X] = lsim(A,B,C,D, p, t, X0);
Wynik całkowania funkcją lsim
Y - wektor odpowiedzi wybranych stanów
X wektory odpowiedzi układu wszystkich stanów
Dane wejściowe do funkcji lsim
A, B, C, D macierze stanu układu dynamicznego
p wektory wymuszeń zewnętrznych
t wektor czasu
X0 wektor warunków początkowych (przemieszczenia i prędkości),
np. X0=[0 0]
20
m=1 %kg
15
k=100 %N/m
10
c=0.6 %Ns/m
5
tk = 4; % czas końcowy [s]
0
0 0.5 1 1.5 2 2.5 3 3.5 4
fs = 1000; % częstotliwość próbkowania [Hz]
t [s]
0.4
dt=1/fs ; % krok czasowy
N = tk*fs; % liczba próbek 0.2
t = linspace(0, tk, N); % wektor czasu
0
-0.2
p = [ utworzyć żądany wektor obciążenia];
-0.4
0 0.5 1 1.5 2 2.5 3 3.5 4
t [s]
XO = [0 0]; % warunki początkowe X0= [uo, vo]
A = [ 0 1
-inv(m)*k -inv(m)*c ];
B= [ 0 ; inv(m)];
C = eye(2);
D= [0;0];
[Y,X] = lsim(A,B,C,D,p,t,XO);
przem=Y(:,1);
predk=Y(:,2);
przysp = -inv(m)*c*predk-inv(m)*k*przem+inv(m)*p';
________________________________________________________________________________________________
2007-02-26 50
p(t) [N]
u(t) [m]
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
c) metoda różnic centralnych
Siła wymuszająca jest dana w postaci zbioru dyskretnych wartości pi = p(ti ) w punktach czaso-
wych ti . Przedział czasu "ti = ti+1 - ti zwykle przyjmowany jest jako stały. Odpowiedz układu okre-
ślona jest w dyskretnych chwilach czasu ti : Równanie ruchu w chwili ti jest następującej postaci:
mui + cui + kui = pi (6.22)
Podstawą metody jest zastąpienie pochodnych (prędkości i przyspieszeń) centralnymi ilorazami
różnicowymi. Dla stałego kroku czasowego "t różnice centralne przyjmują postać:
ui+1 - ui-1 ui+1 - 2ui + ui-1
ui = , ui = . (6.23)
2
2"t
"t
( )
Podstawiając przybliżone wartości (6.23) do równania (6.22) otrzymujemy:
ui+1 - 2ui + ui-1 ui+1 - ui-1
m + c + kui = pi (6.24)
2
2"t
"t
( )
Równanie (6.24) można przekształcić do postaci:
do obliczenia znane znane znane
Ą#ń# Ą#ń# Ą# ń#
m c m c 2m
ó#Ą# ui-1 k ui (6.25)
+ ui+1 = pi - ó# - Ą# - ó# - Ą#
22 2
2"t 2"t
"t "t "t
ó#( ) Ą# ó#( ) Ą# ó# ( ) Ą#
Ł# Ł#Ś# Ł# Ś#
Ś#
k pi
lub:
kui+1 = pi (6.26)
Nieznana wartość przemieszczenia w chwili ti+1 można obliczyć ze wzoru:
pi
ui+1 = (6.27)
k
Rozwiązanie ui+1 w chwili ti+1 jest określone z równania ruchu w chwili ti bez korzystania z rów-
nania ruchu w chwili ti+1 ). Stąd metodę różnic centralnych nazywamy metodą bezpośrednią (jaw-
ną). Do określenia ui+1 wymagane są wartości ui oraz ui-1.
________________________________________________________________________________________________
2007-02-26 51
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Do określenia przemieszczenia u1 wymagane są wartości u0 oraz u-1. Wartość u0 znana jest z wa-
runków początkowych. Aby wyznaczyć wartość u-1 korzystamy z zależności (6.23) dla czasu i=0:
u1 - u-1 u1 - 2u0 + u-1
u0 = , u0 = , (6.28)
2
2"t
"t
( )
Z równań (6.28) wyznaczamy:
2
"t
u-1 = u0 - "t u0 + u0 . (6.29)
( )( )
2
Przyspieszenie u0 można określić z równania ruchu (6.22) dla czasu t0 :
mu0 + cu0 + ku0 = p0 (6.30)
p0 - cu0 - ku0
u0 = (6.31)
m
Metoda różnic skończonych wymaga, aby krok czasowy "t spełniał warunek:
"t 1
< (6.32)
Tn Ą
W praktyce, aby wynik obliczeń były bardziej dokładne przyjmuje się "t /Tn d" 0.1.
Algorytm metody różnic skończonych:
1. Obliczenia początkowe
p0 - cu0 - ku0
1.1. u0 =
m
2
"t
1.2. u-1 = u0 - "t u0 + u0
( )( )
2
m c
1.3. k = +
2
2"t
"t
( )
m c
1.4. a = -
2
2"t
"t
( )
2m
1.5. b = k -
2
"t
( )
2. Obliczenia dla kroku czasowego i
2.1. pi = pi - aui-1 - bui
pi
2.2. ui+1 =
k
ui+1 - ui-1 ui+1 - 2ui + ui-1
2.3. jeżeli wymagane: ui = , ui =
2
2"t
"t
( )
3. Powtórzenie dla następnego kroku czasowego zamień i na i+1, powtórz 2.1, 2.2, 2.3 dla
następnego kroku czasowego (i = 0,1,2,& )
________________________________________________________________________________________________
2007-02-26 52
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 6.2
Układ o jednym stopniu swobody (m = 0.2533 kg, k = 10 N/m, Tn = 1 s., = 0.05 ) poddany jest
działaniu siły wymuszającej jak na rysunku. Obliczyć metodą różnic centralnych przemieszczenie
układu. Krok czasowy "t = 0.1 s., czas końcowy 1 s, u(0) = 0, u(0) = 0 .
1. Obliczenia początkowe
p0 - cu0 - ku0
1.1. u0 =
m
2
"t
1.2. u-1 = u0 - "t u0 + u0
( )( )
2
m c
1.3. k = +
2
2"t
"t
( )
m c
1.4. a = -
2
2"t
"t
( )
2m
1.5. b = k -
2
"t
( )
2. Obliczenia dla kroku czasowego i
2.1. pi = pi - aui-1 - bui
pi
2.2. ui+1 =
k
teoretycznie
ti pi ui-1 ui ui+1 (2.2)
pi (2.1)
ui+1
0.0328
0.2332
0.6487
1.1605
1.5241
1.4814
0.9245
0.0599
-0.7751
-1.2718
-1.2674
________________________________________________________________________________________________
2007-02-26 53
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
d) metoda Newmarka
W 1959r. Newmark opublikował zbiór metod bazujących na następujących równaniach:
ui+1 = ui + Ą# 1+ ł "tń# ui + ł"t ui+1
( ) ( )
Ł#Ś#
(6.33)
2
Ą# ń# Ą# 2 ń#
ui+1 = ui + "t ui + 0.5 - "t ui + "t ui+1
( ) ( )( ) ( )
Ł# Ś# Ł# Ś#
Parametry i ł określają zmianę przyspieszenia w jednym kroku czasowym oraz charakteryzują
dokładność i stabilność metody. Następujący wybór parametrów: " 1/ 6,1/ 4 , ł =1/ 2 jest sa-
tysfakcjonujący ze wszystkich punktów widzenia metod numerycznych, włączając dokładność obli-
czeń.
Algorytm metody Newmarka dla układów liniowych:
Przypadki specjalne: (a) metoda średniego przyspieszenia ( =1/ 4, ł = 1/ 2 )
(b) metoda liniowego przyspieszenia ( =1/ 6, ł =1/ 2 )
1. Obliczenia początkowe
p0 - cu0 - ku0
1.1. u0 =
m
1 ł
1.2. k = k + m + c
2
"t
"t
( )
m ł c 1 # ł ś#
1.3. a = + , b = m + "t -1ź#c
ś#
"t 2 2
# #
2. Obliczenia dla kroku czasowego i
2.1. " = "pi + aui + bui
pi
"
pi
2.2. "ui =
k
łł# ł ś#
2.3. "ui = "ui - ui + "t
ś#1- ź#u i
"t 2
# #
11 1
2.4. "ui = "ui - ui - ui
2
"t 2
"t
( )
2.5. ui+1 = ui + "ui, ui+1 = ui + "ui, ui+1 = ui + "ui
3. Powtórzenie dla następnego kroku czasowego zamień i na i+1, powtórz 2.1 do 2.5 dla
następnego kroku czasowego
Metoda Newmarka jest stabilna, jeżeli:
"t 1 1
d" . (6.34)
Tn
Ą 2 ł - 2
Dla parametrów = 1/ 4, ł =1/ 2 warunek ma postać: "t /Tn d" "
natomiast dla parametrów = 1/ 6, ł = 1/ 2 : "t /Tn d" 0.551
Aby obliczenia były dokładne, należy przyjąć przedział czasowy znacznie mniejszy niż wynika to z
warunku (6.34).
________________________________________________________________________________________________
2007-02-26 54
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 6.3
Układ o jednym stopniu swobody (m = 0.2533 kg, k = 10 N/m, Tn = 1 s., = 0.05 ) poddany jest
działaniu siły wymuszającej jak na rysunku. Obliczyć metodą Newmarka (średniego przyspiesze-
nia) przemieszczenie układu. Krok czasowy "t = 0.1 s., czas końcowy 1 s., u(0) = 0, u(0) = 0 .
1. Obliczenia początkowe
p0 - cu0 - ku0
1.1. u0 =
m
1 ł
1.2. k = k + m + c
2
"t
"t
( )
m ł c 1 # ł ś#
1.3. a = + b = m + "t -1ź#c
ś#
"t 2 2
# #
2. Obliczenia dla kroku czasowego i
2.1. " = "pi + aui + bui
pi
"
pi
2.2. "ui =
k
łł# ł ś#
2.3. "ui = "ui - ui + "t
ś#1- ź#u i
"t 2
# #
11 1
2.4. "ui = "ui - ui - ui
2
"t 2
"t
( )
2.5. ui+1 = ui + "ui, ui+1 = ui + "ui, ui+1 = ui + "ui
3. Powtórzenie dla następnego kroku czasowego
teore-
ti pi ui "pi ui
"
pi "ui "ui "ui ui
tycznie
(2.5) (2.2) (2.3) (2.4) (2.5) (2.5)
(2.1)
ui+1
0.0328
0.2332
0.6487
1.1605
1.5241
1.4814
0.9245
0.0599
-0.7751
-1.2718
-1.2674
________________________________________________________________________________________________
2007-02-26 55
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
7. Drgania swobodne układów dyskretnych o n stopniach swobody.
Równanie macierzowe ruchu:
Mu + Cu + Ku = P(t) (7.1)
Rozpatrywać będziemy drgania swobodne bez tłumienia:
Mu + Ku = 0 (7.2)
Rozwiązanie u(t) będzie w formie:
u(t) = qn (t)Ćn (7.3)
gdzie:
qn(t) = An cosnt + Bn sinnt (7.4)
Wektor Ćn odpowiada za kształt odpowiedzi a qn (t) są współrzędnymi modalnymi. Po podstawie-
niu (7.3) i (7.4) do (7.2) otrzymujemy:
2
K -nM Ćn = 0 (7.5)
( )
Równanie (7.5) posiada nietrywialnie rozwiązanie, jeżeli
2
det K -nM = 0 (7.6)
( )
2
N pierwiastków n równania (7.6) określa N naturalnych częstości układu n (n = 1,2,& N). Każ-
dej częstości n odpowiada wektor własny Ćn obliczony z dokładnością do stałego mnożnika.
Wektor dany przez relatywne wartości N przemieszczeń Ćjn (j = 1,2,& N) określa postać drgań.
N wartości własnych i N wektorów własnych (postaci drgań) można złożyć w macierze. Macierz,
której kolumny składają się z N wektorów własnych nazywa się macierzą modalną:
Ć11 Ć12 Ć1N
Ą# ń#
ó#Ć Ć22 Ć2N Ą#
Ą#
Ą# ń#ó# 12
Ś = = (7.7)
Ł#Ćjn Ś#ó#
Ą#
ó# Ą#
ĆN ĆNN Ś#
Ł#ĆN1 2
2
N wartości własnych n może być złożonych w diagonalną macierz widmową układu:
2
Ą# ń#
1 0 0
ó# Ą#
2
0 2 0
ó# Ą#
&!2 = (7.8)
ó# Ą#
ó# Ą#
2
0 0 N Ś#
ó# Ą#
Ł#
________________________________________________________________________________________________
2007-02-26 56
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
2
Każda wartość własna n i każdy wektor własny Ćn spełnia równanie:
2
KĆn = MĆnn (7.9)
Użycie macierzy modalnej i macierzy widmowej pozwala na zapisanie N równań (7.9) w postaci
jednego równania macierzowego:
KŚ = MŚ&!2 (7.10)
Wektory własne obliczone z dokładnością do stałego mnożnika można znormalizować tak, aby
największa współrzędna wektora była 1:
Ćn
Ćn = (7.11)
max Ćjn
2
Jeżeli wszystkie wartości własne n ą rzeczywiste to wektory własne odpowiadające różnym
częstościom n `" k są ortogonalne:
T
Ćn KĆr = 0 (7.12)
T
Ćn MĆr = 0 (7.13)
Ortogonalność postaci drgań zapewnia, że następujące macierze są diagonalne:
K = ŚTKŚ (7.14)
M=ŚTMŚ (7.15)
Elementy leżące na diagonali:
Kn = ĆnTKĆn (7.16)
Mn = ĆnTMĆn (7.17)
Możliwe jest również takie znormalizowanie wektorów własnych:
Ćn
Ćn = (7.18)
TMĆn
Ćn
że macierz modalna będzie diagonalizować macierz mas do macierzy jednostkowej, a macierz
sztywności do macierzy widmowej:
ŚTKŚ = &!2 (7.19)
ŚTMŚ = I (7.20)
________________________________________________________________________________________________
2007-02-26 57
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Wektor przemieszczenia może być wyrażony we współrzędnych modalnych:
N
u = qr = Śq (7.21)
"Ć
r
r=1
T
gdzie q = [q1 q2 ... qn]T oznacza współrzędne modalne. Mnożąc (7.21) przez Ćn M :
N
TT
Ćn Mu = MĆrqr (7.22)
"Ć
n
r=1
Wykorzystując właściwość ortogonalności:
TT
Ćn Mu = Ćn MĆnqn (7.23)
TT
Ćn Mu Ćn Mu
qn == (7.24)
T
Ćn MĆn Mn
W przypadku obliczeń ręcznych niejednokrotnie łatwiej jest posłużyć się macierzą podatności F,
gdzie
F = K-1 (7.25)
Przekształcając równanie (7.2) poprzez lewostronne pomnożenie przez macierz K-1 otrzymujemy:
FMu + u = 0 (7.26)
Podstawiając (7.4) do równania (7.26) otrzymamy
2
I -nFM Ćn = 0 (7.27)
( )
Równanie (7.27) ma nietrywialnie rozwiązanie, jeżeli:
2
det I -nFM = 0 (7.28)
( )
2
N pierwiastków n równania (7.28) określa N naturalnych częstości układu n (n = 1,2,& N). Każ-
dej częstości n odpowiada wektor własny Ćn .
________________________________________________________________________________________________
2007-02-26 58
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 7.1
Obliczyć częstości i postacie drgań swobodnych układu belkowego z dwiema masami skupionymi.
Sprawdzić ortogonalność postaci drgań.
4 L3
11 =
9 EI
4 L3
22 =
9 EI
7 L3
12 =
18 EI
7 L3
21 =
18 EI
4 7
Ą#ń#
ó#Ą#
Ą# ń#
L3 1 L3 8 7
9 18
Macierz podatności F ==
ó#Ą#
ó# Ą#
7 4 EI 18 EI 7 8
ó#Ą# Ł# Ś#
ó#18 9 Ą#
Ł#Ś#
m 01 0
Ą#ń# Ą# ń#
Macierz mas M ó#Ą# ó# Ą#
== m
0 2 m 0 2
Ł#Ś# Ł# Ś#
I-n2FM Ćn = 0
()
# 1 0
Ą# ń# m L3 8 14 ś# n2mL3
Ą# ń#
=
ś#
n
ó#0 1Ą# -n2 ó#7 16Ą# ź#Ć = 0
18 EI 18EI
Ł# Ś# Ł# Ś#
# #
1-8 -14
Ą#ń#
det = 0
ó#
-7 1-16Ą#
Ł#Ś#
1 = 0.0441
2 = 0.7559
________________________________________________________________________________________________
2007-02-26 59
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
EI
12 = 0.7938
mL3
EI
22 =13.6062
mL3
# 1 0 8 14 ś# 0
Ą# ń# Ą# ń# Ą# ń#
ś#
n
ó#0 1Ą# - n ó#7 16Ą# ź#Ć = ó#0Ą#
Ł# Ś# Ł# Ś# Ł# Ś#
# #
1 = 0.0441 2 = 0.7559
1- 8"0.0441 -14"0.0441 Ć11 0 1- 8"0.7559 -14"0.7559 Ć12 0
Ą#ń# Ą# ń# Ą# ń# Ą# ń# Ą# ń# Ą# ń#
= =
ó# Ą# ó#
-7 "0.0441 1-16"0.0441Ą# ó#Ć21Ś# ó#0Ą# -7 "0.7559 1-16"0.7559Ą# ó#Ć22 Ą# ó#0Ą#
Ł#Ś# Ł# Ł# Ś# Ł# Ś# Ł# Ś# Ł# Ś#
Ć21 = 1 Ć12 =1
Ć11 = 0.954 Ć22 = -0.477
0.954 1
Ą# ń# Ą# ń#
Ć1 = Ć2 =
ó# Ą# ó# Ą#
1
Ł# Ś# Ł#-0.477Ś#
EI
1 = 0.891
mL3
EI
2 = 3.689
mL3
________________________________________________________________________________________________
2007-02-26 60
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 7.2
Obliczyć częstości oraz postacie drgań własnych układu dwóch prętów kratowych obciążonych
masą skupioną m. E = 210 GPa, Q = 10kN. Sprawdzić ortogonalność postaci drgań.
n
zik zim
( )l
mk =
"
EAi i
i=1
m
11 =1.33862"10-7
N
m
12 = 0.50791"10-7
N
m
22 = 0.38095"10-7
N
1.33862 0.50794 1019.368 0
Ą#ń# Ą# ń#
F = =
ó#0.50794 0.38095Ą# "10-7 M ó# 0 1019.368Ą#
Ł#Ś# Ł# Ś#
I Ćn = 0
-n2FM
()
#
1 0 Ą#ń#
Ą# ń# 1364.54638"10-7 517.77778"10-7 ś#
ś#ź#Ćn = 0
ó#
ó#0 1Ą# -n2
ś#
517.77778"10-7 388.32823"10-7 Ą# ź#
Ł# Ś#
Ł#Ś#
# #
1 = 6297.1656
2 = 60658.0517
________________________________________________________________________________________________
2007-02-26 61
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
rad
1 = 79.35
s
rad
2 = 246.29
s
1 = 6297.1656 2 = 60658.0517
1- 0.85928 -0.32605 Ć11 0 1- 8.27707 -3.14074 Ć12 0
Ą#ń# Ą# ń# Ą# ń# Ą# ń# Ą# ń# Ą# ń#
= =
ó# Ą# ó#
-0.32605 1- 0.24454Ą# ó#Ć21Ś# ó#0Ą# -3.14074 1- 2.35552Ą# ó#Ć22 Ą# ó#0Ą#
Ł#Ś# Ł# Ł# Ś# Ł# Ś# Ł# Ś# Ł# Ś#
Ć11 =1 Ć12 = -0.43
Ć21 = 0.43 Ć22 =1
1
Ą# ń# Ą#-0.43
ń#
Ć1 = Ć2 =
ó#0.43Ą# ó# Ą#
1
Ł# Ś# Ł# Ś#
rad rad
1 = 79.35 2 = 246.29
s s
________________________________________________________________________________________________
2007-02-26 62
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 7.3
Dla 2-kondygnacyjnej ramy obliczyć częstotliwości i postacie drgań własnych.
E = 200 GPa
m2
u2
EIb= "
h = 3 m
Ic = 2500"10-8 m4
EIc
h
EIc
m1
u1
EIb= "
m1 = 2000 kg
m2 = 1000 kg
EIc
h
EIc
2h
24EIc
k1 =
h3
24EIc
k2 =
h3
k1 + k2 -k2
Ą#ń#
K=ó#
-k2 k2 Ą#
Ł#Ś#
m1 0
Ą#ń#
M ó#
=
0 m2 Ą#
Ł#Ś#
[mode,vale]=eig(K,M);
val=diag(vale);
omega=sqrt(val);
f=omega/2/pi
f =
5.7423
13.8631
mode =
-0.7071 -0.7071
-1.0000 1.0000
________________________________________________________________________________________________
2007-02-26 63
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 7.4
Obliczyć częstości drgań wspornika; m = 40 kg, L = 3 m, E = 205 GPa, I = 250 cm4.
192EI -96EI 24EI
Ą#ń#
0
ó#
L3 L3 L2 Ą#
ó#Ą#
ó#-96EI 96EI -24EI -24EI Ą#
ó#
L3 L3 L2 L2 Ą#
K =
ó#Ą#
-24EI 16 4
ó#Ą#
0
ó# L2 LL Ą#
ó#Ą#
24EI -24EI 4 8
ó#Ą#
Ł# L2 L2 LL Ś#
mL
Ą#
0 0 0ń#
ó#Ą#
2
ó#Ą#
mL
ó#
00 0Ą#
M =
ó#Ą#
4
ó#
0 0 0 0Ą#
ó#Ą#
ó#Ą#
0 0 0 0Ś#
Ł#
________________________________________________________________________________________________
2007-02-26 64
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Sposoby postępowania:
1) uzupełnienie bezmasowych stopni swobody małymi masami
przyjęto M(3,3)=M(4,4)=m*10-6
[mode,vale]=eig(K,M)
val=diag(vale)
omega=sqrt(val)
mode =
0.3274 1.0000 0.0000 0.0000
1.0000 -0.6547 -0.0000 -0.0000
0.3792 0.0987 -0.4142 1.0000
0.4830 -1.7041 1.0000 0.4142
omega =
39.6957
204.4760
164592.101
274607.936
2) kondensacja względem bezmasowych stopni swobody:
768 -240 mL
Ą#ń# Ą# ń#
0
ó#Ą# ó# Ą#
EI
7 7 2
K ' = M =
ó#Ą# ó# Ą#
-240 96 L3 mL
ó#Ą# ó# Ą#
0
ó#Ą# ó# Ą#
Ł# 7 7 Ś# Ł# 4 Ś#
mode =
-0.3274 1.0000
-1.0000 -0.6547
omega =
39.6957
204.4761
________________________________________________________________________________________________
2007-02-26 65
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 7.5
Dla wspornika z zadania 7.4 przedstawić na rysunku rozkład modalny wektora przemieszczeń
u = [1 1]T .
mL
Ą#ń#
0
ó#Ą#
0.3274 1
Ą# ń# Ą#ń#
2
Ć1 = Ć2 = M =
ó#Ą#
ó# Ą# ó#Ą#
1 mL
Ł# Ś# Ł#-0.6547Ś# ó#Ą#
0
ó#Ą#
Ł# 4 Ś#
mL
Ą#ń#
0
ó#Ą#
1
Ą# ń#
2
0.3274 1
[]
ó#Ą#
ó#1Ą#
mL
Ł# Ś#
ó#Ą#
0
ó#Ą#
Ł# 4 Ś#
q1 ==
mL
Ą#ń#
0
ó#Ą#
0.3274
Ą# ń#
2
0.3274 1
[]
ó#Ą#
ó# Ą#
mL 1
Ł# Ś#
ó#Ą#
0
ó#Ą#
Ł# 4 Ś#
mL
Ą#ń#
0
ó#Ą#
1
Ą# ń#
2
1
[ -0.6547
]
ó#Ą#
ó#1Ą#
mL
Ł# Ś#
ó#Ą#
0
ó#Ą#
Ł# 4 Ś#
q2 ==
mL
Ą#ń#
0
ó#Ą#ń#
1
Ą#
2
1
[ -0.6547
]
ó#Ą#Ą#
ó#
mL
Ł#-0.6547Ś#
ó#Ą#
0
ó#Ą#
Ł# 4 Ś#
u = q1Ć1 + q2Ć2
________________________________________________________________________________________________
2007-02-26 66
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 7.6
Unormować wektory własne tak, by diagonalizowały macierz mas do macierzy jednostkowej.
3 0 5 2
Ą# ń# Ą# ń# Ą# ń#
M = Ć2 =
ó#0 4Ą# Ć1 = ó#5Ą# ó# Ą#
Ł# Ś# Ł# Ś# Ł#-1.5Ś#
Sprawdzenie ortogonalności wektorów własnych:
3 0 2
5 5 = 0
[ ]Ą# ń# Ą# ń#
ó#0 4Ą# ó# Ą#
Ł# Ś# Ł#-1.5Ś#
Ć1TMĆ1 =175
Ć2TMĆ2 = 21
5 2
Ą#ń#
ó#Ą#
175 21
ó#Ą#
Ś =
5 -1.5Ą#
ó#
ó#Ą#
175 21
Ł#Ś#
5 2 5 2
Ą#ń# Ą#ń#
ó#Ą# ó#Ą#
3 0
Ą# ń#
175 21 175 21
ó#Ą# ó#Ą#
ŚTMŚ == &
5 -1.5Ą# ó#0 4Ą# ó# 5 -1.5Ą#
ó#
Ł# Ś#
ó#Ą# ó#Ą#
175 21 175 21
Ł#Ś# Ł#Ś#
1 0
Ą# ń#
& =
ó#0 1Ą#
Ł# Ś#
________________________________________________________________________________________________
2007-02-26 67
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
8. Drgania wymuszone harmonicznie układów dyskretnych o n stopniach
swobody.
Drgania wymuszone harmonicznie układów dyskretnych o n stopniach swobody.
Równanie macierzowe ruchu:
Mu + Cu + Ku = P(t) (8.1)
Rozpatrywać będziemy drgania bez tłumienia:
Mu + Ku = P0 sint (8.2)
Rozwiązanie u(t) będzie w formie:
u(t) = Ć sint (8.3)
Po podstawieniu (7.3) do (7.2) otrzymujemy:
K -2M Ć = P0 (8.4)
( )
Kd
KĆ = P0 + 2MĆ (8.5)
Pz
gdzie Pz jest zastępczą siłą statyczną, a Kd dynamiczną macierzą sztywności.
Amplitudy drgań obliczymy z równania:
Ć = Kd -1P0 (8.6)
________________________________________________________________________________________________
2007-02-26 68
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 8.1
Obliczyć amplitudy drgań Ć dla = 39 rad/s oraz dla = 40 rad/s , E = 205GPa, I = 250cm4,
m = 100kg, L = 3m, p0 = 1 kN.
L3
11 =
24EI
L3
22 =
3EI
5L3
12 =
48EI
macierz podatności : macierz mas :
2 5 4 0 1
L3 Ą#ń# m Ą# ń# Ą# ń#
F= M P = " p0
=
ó#Ą# ó# Ą# ó#0Ą#
0
48 EI 5 16 4 0 1
Ł#Ś# Ł# Ś# Ł# Ś#
macierz sztywności:
16
ń#
48 EI Ą# - 5
K
=
Ą#
7 L3 ó# - 5 2
Ł#Ś#
-1
K Ć = P0 ! Ć = K P0 = K-1P0
-2M -2M
() ()
d
16
Ą#ń#Ą# ń#
48EI -5 4 0
2m
Kd =-
Ą#ó#0
7L3 ó# 2 4 1Ą#
Ł#-5 Ś#Ł# Ś#
________________________________________________________________________________________________
2007-02-26 69
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
rad
= 39
s
1.9304 -0.6508
Ą#ń# N
Kd ="106
ó#Ą#
m
Ł#-0.6508 0.2223 Ś#
0.0398
Ą# ń#
Ć = K-1P0 = m
d ó#0.1164Ą#
Ł# Ś#
rad
= 40
s
1.9225 -0.6508
Ą#ń# N
Kd ="106
ó#Ą#
m
Ł#-0.6508 0.2203 Ś#
6.0058
Ą#ń#
Ć = K-1P0 = m
d ó#17.7404Ą#
Ł#Ś#
________________________________________________________________________________________________
2007-02-26 70
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 8.2
Wyznaczyć amplitudy drgań Ć oraz narysować wykres momentów dynamicznych.
EI = 16000kNm2, m1 = 4Mg, m2 = 1.2Mg, = 31 rad/s, L = 1 m.
8L3
11 =
3EI
26L3
22 =
3EI
4L3
12 =-
3EI
8 4
Ą#ń#
-
L3 ó# 3 3 Ą#
F
=
ó#Ą#
EI 4 26
ó#Ą#
-
ó#Ą#
3 3
Ł#Ś#
26 4
Ą#ń#
ó#Ą#0.40625 0.0625ń#
EI 3 EI Ą#
3 3
K ==
ó#Ą#
Ą#
L3 64 4 8 L3 ó# 0.0625 0.125
ó#Ą# Ł#Ś#
ó#Ą#
3 3
Ł#Ś#
4000 0
Ą#ń# 1600
M ó#Ą# Ą# ń#
= P =
0 ó# Ą#
0 1200 0
Ł# Ś#
Ł#Ś#
K Ć = P0
-2M
()
2656000 1000000
Ą#ń#
K =
-2M ó#Ą#
Ł#1000000 846800 Ś#
________________________________________________________________________________________________
2007-02-26 71
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
1.0847
Ą#ń#
Ć ="10-3 m
ó#Ą#
Ł#-1.2809Ś#
Zastępcze siły statyczne Pz = P0 +2MĆ
5.7695
Ą#ń#
Pz = kN
ó#Ą#
Ł#-1.4772Ś#
Mdyn = Pz1 " M1 + Pz 2 " M2
________________________________________________________________________________________________
2007-02-26 72
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
9. Drgania swobodne tłumione układów dyskretnych o n stopniach swobo-
dy.
Tłumienie:
o Klasyczne (proporcjonalne) - np. tłumienie Rayleigh a
C = a0M + a1K (9.1)
gdzie a0, a1 są stałymi.
o Nieklasyczne (nieproporcjonalne)
Drgania swobodne tłumione w układach n-stopni swobody.
Równanie macierzowe ruchu:
Mu + Cu + Ku = 0 (9.2)
Wektor przemieszczenia wyrażony we współrzędnych modalnych:
N
u = qr = Śq (9.3)
"Ć
r
r=1
Jeżeli macierz C jest macierzą tłumienia proporcjonalnego to macierz modalna Ś (wyznaczona dla
układu bez tłumienia) diagonalizuje ją:
C=ŚTCŚ (9.4)
Równanie (9.2) we współrzędnych modalnych:
Mq + Cq + Kq = 0 (9.5)
Otrzymujemy N równań różniczkowych
Mnqn + Cnqn + Knqn = 0 (9.6)
Liczba tłumienia każdej postaci:
Cn
n = (9.7)
2Mnn
Równanie (9.6) można przekształcić do postaci:
2
qn + 2nnqn +n qn = 0 (9.8)
Równanie (9.8) jest takie samo jak rówanie dla jednego stopnia swobody. Rozwiązanie jest nastę-
pujące:
________________________________________________________________________________________________
2007-02-26 73
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
#ś#
qn (0) + qn (0)nn
n
qn(t) = e- nt ś# qn (0)cosnDt + sinn Dt (9.9)
ź#
nD
# #
Przemieszczenie:
NN
#ś#
qn (0) + qn (0)nn
n
u = qn = e- nt ś# qn(0)cosnDt + sinn Dt (9.10)
"Ć "Ć #ź#
n n
nD
n=1 n=1
#
Rozwiązanie problemu własnego z uwzględnieniem tłumienia.
Rozwiązanie problemu własnego z tłumieniem nastąpi po zamianie równanie różniczkowego dru-
giego stopnia:
u + M-1Cu + M-1Ku = 0 (9.11)
na równanie pierwszego stopnia:
x=Ax (9.12)
gdzie:
u
Ą# ń#
x = (9.13)
ó#uĄ#
Ł# Ś#
Rówanie (9.12) można przedstawić jako:
u0 I u
Ą# ń# Ą# ń# Ą# ń#
= (9.14)
ó#uĄ# ó# -1 ó#uĄ#
Ł# Ś# Ł#-M K -M-1CĄ# Ł# Ś#
Ś#
________________________________________________________________________________________________
2007-02-26 74
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 9.1
Wyprowadzić macierze tłumienia proporcjanalnego:
1
1 =12.57 1 = 5%
s
400 0 0
Ą#ń#
1 1
ó#Ą#
M = 0 400 0 2 = 39.33 1 = 5%
ó#Ą#
386 s
ó#Ą#
0 0 200Ś#
Ł#
1
3 = 46.89
s
2
Ą# -1 0
ń#
ó#
K = 610 2 -1Ą#
ó#-1 Ą#
ó# -1 1
0
Ł#Ą#
Ś#
0.401 0.803 0.401
ż# # ż# # ż# #
# # # # # #
Ć1 = Ć2 = 0 Ć3 =
#0.695Ź# # Ź# #-0.695Ź#
#0.803# # # # #
0.803
# # #-0.803# # #
Rozwiązanie:
C = a0M C = a1K
2
Cn = a0Mn Cn = a1Kn = a1n Mn
Cn Cn
n = n =
2Mnn 2Mnn
a0 1 a1
n = " n = n
2 n 2
a0 1 a1
n = " + n
2 n 2
Współczynniki a0 i a1 określimy dla dwóch znanych liczb tłumienia i i odpowiadających po-
j
staciom drgań i oraz j
a0 1 a1
n = " + i
2 i 2
a0 1 a0
n = " + j
2 j 2
________________________________________________________________________________________________
2007-02-26 75
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
1
Ą#
i ń#
ó#
i Ą# Ą# ń# i
a0 ż# #
1
ó#Ą#
=
ó#a Ą## Ź#
ó# 1 Ą#
2
j
Ł# 1 Ś#
# #
ó# j Ą#
j
Ł#Ś#
1
Ą#
ó#12.57 12.57ń#
Ą#
a0 0.05
Ą# ń# ż# #
=
ó#Ą#
ó#a Ą##0.05Ź#
1
ó# Ł# 1 Ś# # #
ó#39.33 39.33Ą#
Ł#Ą#
Ś#
a0 = 0.9198
a1 = 0.0021
3.35
Ą# -1.3 0
ń#
ó#
C = a0M + a1K =
ó#-1.3 3.55 -1.3Ą#
Ą#
ó# -1.3 1.78Ś#
0
Ł#Ą#
a0 1 a1
3 = + 3 = 0.0593
2 3 2
________________________________________________________________________________________________
2007-02-26 76
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 9.2
Wyznaczyć postacie i częstości drgań z uwzględnieniem tłumienia proporcjonalnego i niepropor-
cjonalnego
2
Ą# -1 2 0 0.4 0.4
ń# Ą# ń# Ą# -0.1
ń# Ą# -0.1
ń#
K = M =
Cn =
ó# Ą# ó#0 3Ą# Cp = ó# Ą# ó# Ą#
Ł#-1 1 Ś# Ł# Ś# Ł#-0.1 0.4 Ś# Ł#-0.1 0.1 Ś#
algorytm:
[mode,ome]=eig(K,M)
omega_nietlum=sqrt(ome)
A=zeros(4,4);
O=zeros(2,2);
I=eye(2,2);
A=[O I;-inv(M)*K -inv(M)*Cp];
[mode_t,ome_t]=eig(A)
ome_t=diag(ome_t)
czestosc_kolowa_nietlum=abs(ome_t)
czestosc_kolowa_tlum=imag(ome_t)
wyniki
mode =
-0.3031 -0.6388
-0.5216 0.2475
omega_nietlum =
0.3737 0
0 1.0926
mode_t =
0.0632 + 0.6264i 0.0632 - 0.6264i -0.4707 -0.4707
-0.0245 - 0.2427i -0.0245 + 0.2427i -0.8099 -0.8099
-0.6878 -0.6878 0.0268 - 0.1738i 0.0268 + 0.1738i
0.2665 - 0.0000i 0.2665 + 0.0000i 0.0461 - 0.2991i 0.0461 + 0.2991i
ome_t =
-0.1097 + 1.0871i
-0.1097 - 1.0871i
-0.0570 + 0.3693i
-0.0570 - 0.3693i
czestosc_kolowa_nietlum =
1.0926
1.0926
0.3737
0.3737
czestosc_kolowa_tlum =
1.0871
-1.0871
0.3693
-0.3693
________________________________________________________________________________________________
2007-02-26 77
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
algorytm:
[mode,ome]=eig(K,M)
omega_nietlum=sqrt(ome)
A=zeros(4,4);
O=zeros(2,2);
I=eye(2,2);
A=[O I;-inv(M)*K -inv(M)*Cn];
[mode_t,ome_t]=eig(A)
ome_t=diag(ome_t)
czestosc_kolowa_nietlum=abs(ome_t)
czestosc_kolowa_tlum=imag(ome_t)
wyniki
mode =
-0.3031 -0.6388
-0.5216 0.2475
omega_nietlum =
0.3737 0
0 1.0926
mode_t =
0.0580 + 0.6276i 0.0580 - 0.6276i 0.4703 - 0.0167i 0.4703 + 0.0167i
-0.0473 - 0.2382i -0.0473 + 0.2382i 0.8099 0.8099
-0.6881 -0.6881 -0.0014 + 0.1759i -0.0014 - 0.1759i
0.2638 - 0.0275i 0.2638 + 0.0275i -0.0131 + 0.3026i -0.0131 - 0.3026i
ome_t =
-0.1005 + 1.0872i
-0.1005 - 1.0872i
-0.0162 + 0.3736i
-0.0162 - 0.3736i
czestosc_kolowa_nietlum =
1.0918
1.0918
0.3739
0.3739
czestosc_kolowa_tlum =
1.0872
-1.0872
0.3736
-0.3736
________________________________________________________________________________________________
2007-02-26 78
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
10. Układy dyskretne o n stopniach swobody dowolne wymusznie.
Równanie ruchu w przestrzeni stanów dla układu o n stopniach swobody
Równanie ruchu układu tłumionego:
Mu + Cu + Ku = p
T
gdzie: u = u1 u2 un - wektor n przemieszczeń opisujących ruch
[ ]
M- macierz mas układu
C - macierz tłumienia
K macierz sztywności
p wektor obciążeń zewnętrznych
Równanie ruchu (4) w przestrzenie stanu przyjmuje postać:
x=Ax+Bp (state space equation of motion)
u
Ą# ń#
gdzie: x= , wektor stanów,
ó#uĄ#
Ł# Ś#
0I
Ą#ń#
A=ó#-M-1 - macierz układu
K -M-1CĄ#
Ł#Ś#
I macierz jednostkowa (n x n)
0 macierz zerowa (n x n)
0
Ą# ń#
B=
ó#M-1 Ą#
Ł# Ś#
Równanie wyjścia:
y =Cx+Dp (output equation)
Macierze A, B, C i D opisują proces dynamiczny.
Uwagi:
1. Wektory własne (także dla układów o proporcjonalnym tłumienie) wyznaczone przy pomocy
macierzy stanu są wektorami zespolonymi.
2. Sformułowanie w przestrzeni stanów jest uogólnieniem klasycznego sformułowanie równania
ruchu.
3. Rozwiązanie układu w czasie można uzyskać stosując w Matlabie funkcję lsim.
4. Dla proporcjonalnych macierzy tłumień wektory własne układu są rzeczywiste.
5. Wyznaczenie wektorów własnych układu wymaga użycia funkcji eig(A) w programie Matlab.
________________________________________________________________________________________________
2007-02-26 79
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Macierz tłumienia C jest proporcjonalna, jeżeli wyznaczamy ją według wzoru:
C = ąK + M
gdzie: ą i to współczynniki proporcjonalności względem K i M
Wszystkie programy komercyjne (dla konstrukcji budowlanych) rozwiązują wyłącznie problemy z
tłumieniem proporcjonalnym.
Wyjaśnienie użycia funkcji lsim w programie Matlab.
Format funkcji
[Yres, Xres] = lsim(A,B,C,D, U, t, X0);
Wynik całkowania funkcją lsim
Yres - wektor odpowiedzi wybranych stanów
Xres wektory odpowiedzi układu wszystkich stanów
Dane wejściowe do funkcji lsim
A, B, C, D macierze stanu układu dynamicznego
U wektory wymuszeń zewnętrznych
t wektor czasu
X0 wektor warunków początkowych (przemieszczenia i prędkości)
Przykład użycia jest podane w programie st_swobo3_w10.m. Wspornik ustawiony pionowo mode-
lowany jest 3 stopniami swobody. Pierwsza masa skupiona znajduje się przy utwierdzeniu. Prze-
mieszczenie 2 opisuje ruch środka wspornika, zaś współrzędna x3 znajduje się na końcu wspornika.
Animację wykonuje podprogram dof3_ani
________________________________________________________________________________________________
2007-02-26 80
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
11. Transformata Fouriera.
Jean Baptiste Joseph Fourier (1768-1830), francuski matematyk w 1807 roku udowodnił, że każda
okresowa funkcja może zostać rozłożona na sumę sinusów i cosinusów.
2
4
2 1.5
0
1
-2
0.5
-4
0 1 2 3 4 5 0
0 1 2 3 4
t [s]
f [Hz]
2 2
2
1 1 1
0 0 0
-1 -1 -1
-2 -2 -2
0 2 4
0 2 4 0 2 4
t [s]
t [s] t [s]
2
0
-2
5
3
4
2 3
2
1
1
0 0 t [s]
f [Hz]
Szeregi Fouriera
Funkcję p(t) nazywamy funkcją okresową z okresem T jeżeli spełnia ona następujący warunek:
p(t + nT) = p(t), n = - ",...,-2,-1,0,1, 2,..., " (11.1)
gdzie n jest liczbą całkowitą. Przykładowe funkcje okresowe pokazane są na poniższym rysunku.
________________________________________________________________________________________________
2007-02-26 81
u( t )
u( t )
u( t )
u( t )
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
10
0
-10
0 2 4 6 8 10
1
0
-1
0 2 4 6 8 10
5
0
-5
0 2 4 6 8 10
Funkcje okresowe
Okresową funkcję można rozłożyć w nieskończony szereg Fouriera będący sumą sinusów i cosinu-
sów różnej częstotliwości:
" "
p(t) = a0 + cos(2Ą fnt) + sin(2Ą fnt) (11.2)
"a "b
n n
n=1 n=1
gdzie p(t) jest funkcją w dziedzinie czasu, a0, an, bn są współczynnikami szeregu Fouriera, fn ozna-
cza częstotliwość fali. Zależność pomiędzy częstotliwością fn mierzoną w Hertzach a częstością
kołową n mierzoną w rad/s opisana jest następująco:
n = 2Ą fn (11.3)
Za pomocą szeregu Fouriera można rozłożyć sygnał na składniki o częstotliwościach będących cał-
kowitymi wielokrotnościami częstotliwości podstawowej (fundamentalnej) f (lub ):
fn = nf , n = n (11.4)
Częstotliwość będąca całkowitą wielokrotnością częstotliwości podstawowej nazywana jest często-
tliwością harmoniczną.
Ostatecznie współczynniki szeregu Fouriera wyrażone są następująco:
T
1
a0 = p(t)dt
+"
T
0
T
2
an = p(t)cos(2Ą fnt)dt (11.5)
+"
T
0
T
2
bn = p(t)sin(2Ą fnt)dt
+"
T
0
________________________________________________________________________________________________
2007-02-26 82
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 11.1
1 dla -Ą d" t < 0
ż#
Daną funkcję p(t) =
#3 dla - 0 d" t < Ą rozwinąć w szereg Fouriera.
#
4
3
2
1
0
-5 0 5
Rozwinięcie okresowe funkcji p(t).
0 Ą
Ą 0 Ą
#ś#
Ą#ń#
11 1
a0 = p(t)dt = dt + 3 = t + 3t = 2
ó#Ą#
+"+" +"dtŚ# 2Ą ś#ź#
ś#ź#
2Ą 2Ą
-Ą Ł#-Ą 0
-Ą 0
# #
0 Ą
ĄĄ
0
Ą#ń#
Ą#ń#
11 1
an = p(t)cos(nt)dt = cos(nt) + 3
ó#Ą#
+"+" +"cos(nt)Ś# = ó#sin(nt) + 3sin(nt) Ą# = 0
ĄĄ nĄ
ó#Ą#
-Ą Ł#-Ą 0
-Ą 0
Ł#Ś#
0 Ą
ĄĄ
0
Ą# ń#
Ą#ń#
11 -1
bn = p(t)sin(nt)dt = sin(nt) + 3
ó#Ą#
+"+" +"sin(nt)Ś# = ó#cos(nt) + 3cos(nt) Ą# =
ĄĄ nĄ
ó# Ą#
-Ą Ł#-Ą 0
-Ą 0
Ł# Ś#
22
= 1- cos(nĄ ) = 1- (-1)n
()
()
nĄ nĄ
Ostatecznie szereg Fouriera przyjmuje postać:
"
2
p(t) = 2 + (1- (-1)n )sin(nt)
"
nĄ
n=1
Przeanalizujmy, ile wyrazów nieskończonego szeregu Fouriera potrzebnych jest do opisania zada-
nej funkcji. W tym celu obliczymy sumę szeregu dla różnych wartości n.
N
2
SN (t) = 2 + (1- (-1)n)sin(nt)
"
nĄ
n=1
________________________________________________________________________________________________
2007-02-26 83
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
4 4
N=2 N=3
3 3
2 2
1 1
0 0
-5 0 5 -5 0 5
4 4
N=10 N=20
3 3
2 2
1 1
0 0
-5 0 5 -5 0 5
Sumy szeregu dla N=2,3,10,20.
Ciągła transformata Fouriera.
Ciągłe przekształcenie Fouriera opisują zależności:
"
P() = p(t)e-itdt (11.6)
+"
- "
"
1
p(t) = P()eitd (11.7)
+"
2Ą
- "
Współczynnik Fouriera P() otrzymuje się poprzez korelację rozpatrywanej funkcji p(t) z falą sinu-
soidalną eit. Jeżeli sygnał zawiera składnik częstotliwościowy wówczas osiąga duże wartości, w
przeciwnym razie równa się zeru. Informacją, jakiej dostarcza nam transformata Fouriera jest, jakie
składniki częstotliwościowe zawarte są w sygnale. Nic więcej, nic mniej.
________________________________________________________________________________________________
2007-02-26 84
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 11.2
Obliczyć transformatę Fouriera funkcji p(t):
1 dla -Ą d" t d" Ą
ż#
p(t) =
#0 dla -Ą > t > Ą
#
2
1.5
1
0.5
0
-0.5
-1
-4 -2 0 2 4
Impuls prostokątny.
Wykorzystując zależność oraz wzór obliczamy całkę Fouriera:
" " ĄĄ
P() = p(t)eitdt = p(t) cos(t) - i sin(t) dt = cos(t)dt - i sin(t)dt =
()
+"+" +" +"
- " - " -Ą -Ą
Ą
sin(x) 2sin(Ą )
==
-Ą
Ponieważ transformata Fouriera mierzy częstotliwości, jakie zawarte są w funkcji, a rozważana jest
funkcja stała, stąd największa wartość P() występuje przy zerowej częstotliwości.
8
6
4
2
0
-2
-10 -5 0 5 10
Całka Fouriera impulsu prostokątnego.
________________________________________________________________________________________________
2007-02-26 85
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
ZADANIE 11.3
Obliczyć transformatę Fouriera funkcji p(t):
cos(4t) dla -Ą d" x d" Ą
ż#
p(t) =
#0
dla -Ą > x > Ą
#
2
1
0
-1
-2
-4 -2 0 2 4
Wykres funkcji p(t).
Wykorzystując zależność oraz wzór obliczamy całkę Fouriera:
" "
P() = p(t)e-itdt = p(t) cos(t) - i sin(t) dt =
()
+"+"
- " - "
ĄĄ
-2 sin(Ą )
= cos(4t) cos(t)dt - i cos(4t)sin(t)dt =
+"+"
16 -2
-Ą -Ą
Największe wartości przekształcenia Fouriera występują dla = 4 i = - 4 czyli dla częstości
funkcji p(t).
4
3
2
1
0
-1
-10 -8 -6 -4 -2 0 2 4 6 8 10
Całka Fouriera funkcji p(t)
________________________________________________________________________________________________
2007-02-26 86
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Dyskretna transformata Fouriera i jej realizacja w postaci FFT
Szeregi Fouriera mają zastosowanie w analizie sygnałów okresowych, zaś transformata
(całka)Fouriera stosowana jest w analizie sygnałów nieokresowych. Dyskretne przekształcenie Fo-
uriera DFT (Discrete Fourier Transform) rozwijano w latach 40. i 50. XX wieku. Szeregi Fouriera
zastąpiono postacią dyskretną, aby możliwe było użycie cyfrowych technik obliczeniowych.
Ciągły przebieg p(t) można spróbkować za pomocą przetwornika analogowo-cyfrowego co
ts sekund otrzymując N próbek. Częstotliwość próbkowania fs mierzy jak często poszukujemy war-
tości funkcji dyskretnej:
1
fs = (11.8)
ts
Zakładamy, że sygnał jest okresowy z okresem T, który otrzymujemy mnożąc liczbę próbek N
przez przedział czasowy pomiędzy poszczególnymi próbkami ts:
T = Nts (11.9)
Zakładamy również, że częstotliwość podstawowa wyrażona jest wzorem:
1 1 fs
f0 = = = (11.10)
T N"t N
Widać zatem, że od liczby przyjętych próbek zależy rozdzielczość analizy.
Dyskretna transformata Fouriera DFT wywodzi się z całki Fouriera:
"
X ( f ) = x(t)e-12Ą ftdt (11.11)
+"
- "
i ma postać:
-i2Ą mn
N -1
N
X (m) = x(n)e (11.12)
"
n=0
gdzie:
n indeks próbek wejściowych w dziedzinie czasu, n = 0,1,...,N-1
m indeks próbek wyjściowych w dziedzinie częstotliwości, m = 0,1,...,N-1
x(n) ciąg próbek wejściowych, x(0), x(1),...
X(m) ciąg próbek wyjściowych DFT, X(0), X(1),...
N liczba próbek ciągu wejściowego i wyjściowego
Realizacja DWT używając algorytmu FFT (Fast Fourier Transform)
t % wektor czasu
a % rozpatrywany sygnał
fs=1/dt % czestotliwosc probkowania
N=length(t)
fo=fs/N % czestotliwosc podstawowa
fa=fft(a); % szybka transformata
base=fo*(0:N/2-1); % wyznaczenie osi czestotliwosci
mag=abs(fa(1:N/2)); % wyznaczenie widma
A=2*mag./N; % normalizacja odpowiedzi
________________________________________________________________________________________________
2007-02-26 87
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Przykłady:
Drgania silosu podczas wypływu materiału sypkiego:
40
30
20
10
0
-10
-20
-30
-40
0 5 10 15 20 25 30
time [s]
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 50 100 150 200 250
f [Hz]
________________________________________________________________________________________________
2007-02-26 88
2
v
a [m/s ]
normalized ampitude FFT
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Drgania płyty z pleksi wywołane uderzeniem za pomocą młotka modalnego:
40
30
20
10
0
-10
-20
-30
-40
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
time [s]
0.2
0.18
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
0 100 200 300 400 500 600 700 800 900
f [Hz]
________________________________________________________________________________________________
2007-02-26 89
2
v
a [m/s ]
normalized ampitude FFT
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Drgania płyty stalowej wywołane uderzeniem za pomocą młotka modalnego:
40
30
20
10
0
-10
-20
-30
-40
0 0.5 1 1.5 2 2.5 3 3.5 4
time [s]
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 100 200 300 400 500 600 700 800 900
f [Hz]
________________________________________________________________________________________________
2007-02-26 90
2
v
a [m/s ]
normalized ampitude FFT
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Siła pochodząca z młotka modalnego (uderzenie w płytę stalową):
16
15
14
10
12
5
10
0
8 0.198 0.2 0.202 0.204 0.206 0.208
time [s]
6
4
2
0
-2
0 1 2 3 4 5 6 7
time [s]
0.045
0.04
0.035
0.03
0.025
0.02
0.015
0.01
0.005
0
0 50 100 150 200 250 300
frequency [Hz]
________________________________________________________________________________________________
2007-02-26 91
force [N]
force [N]
nor mal i z ed ampl i t ude FFT
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Załącznik A całkowanie graficzne
________________________________________________________________________________________________
2007-02-26 92
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Załącznik B wyjściowe siły przywęzłowe
________________________________________________________________________________________________
2007-02-26 93
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Załącznik C element belkowy
Kondensacja eliminuje nieistotne, niezerowe przemieszczenia, którym odpowiadają zerowe war-
tości sił przywęzłowych.
Kondensacja układu Kq =P względem podwektora q0 zawartego w wektorze q .
K11 K12 q ' R1
Ą# ń# Ą# ń# Ą# ń#
Kq == (1)
ó#K K22 Ą# ó#q Ą# ó#R Ą#
Ł# 21 Ś# Ł# 0 Ś# Ł# 2 Ś#
Rozpisując (1) na 2 równania otrzymujemy:
K11q '+ K12q0 = R1
(2)
K21q '+ K22q0 = R2 ! q0 = K-1(R2 - K21q ')
22
A następnie:
K11q '+ K12K-1R2 - K12K-1K21q ' = R1
22 22
(K11 - K12K-1K21 ' = R1 - K12K-1R2 (3)
22 22
)q
K ' P'
Macierz skondensowana K ' ostatecznie przyjmuje postać:
K ' = K11 - K12K-1K21 (4)
22
Przykłady:
a)Kondensacja względem fa
va a vb b
va
K21
K22 a
K12
vb
K11
b
Kondensacja względem fa oraz fb
vb a b
K11 K12
vb
a
K22
K21
b
________________________________________________________________________________________________
2007-02-26 94
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Macierz sztywności elementu belkowego:
f
k
f
i
v
k
v
i
12EI 6EI -12EI 6EI
Ą#ń#
ó#
l3 l2 l3 l2 Ą#
ó#Ą#
6EI 4EI -6EI 2EI
ó#Ą#
ó#Ą#
l2 l l2 l
K = (5)
ó#Ą#
ó#-12EI -6EI 12EI -6EI Ą#
ó# l3 l2 l3 l2 Ą#
ó#Ą#
6EI 2EI -6EI 4EI
ó#Ą#
Ł# l2 l l2 l Ś#
Macierz skondensowana elementu belkowego:
f
k
v
k
v
i
3EI -3EI 3EI
Ą# ń#
ó#
l3 l3 l2 Ą#
ó# Ą#
-3EI 3EI -3EI
ó# Ą#
K = (6)
ó#
l3 l3 l2 Ą#
ó# Ą#
3EI -3EI 3EI
ó# Ą#
ó# Ą#
l2 l2 l
Ł# Ś#
Macierz skondensowana elementu belkowego:
f
i
v
v
i
k
3EI 3EI 3EI
Ą# ń#
-
ó#
l3 l2 l3 Ą#
ó# Ą#
3EI 3EI -3EI
ó# Ą#
K = (7)
ó#
l2 l l2 Ą#
ó# Ą#
-3EI -3EI 3EI
ó# Ą#
ó# Ą#
l3 l2 l3 Ś#
Ł#
________________________________________________________________________________________________
2007-02-26 95
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Kondensacja względem i :
f
k
f
i
v
k
v
i
vi i vk k
12EI 6EI -12EI 6EI
Ą#ń#
ó#
l3 l2 l3 l2 Ą#
ó#Ą#
6EI 4EI -6EI 2EI
ó#Ą#
ó#Ą#
l2 l l2 l
K =
ó#Ą#
ó#-12EI -6EI 12EI -6EI Ą#
ó# l3 l2 l3 l2 Ą#
ó#Ą#
6EI 2EI -6EI 4EI
ó#Ą#
Ł# l2 l l2 l Ś#
K ' = K11 - K12K-1K21
22
12EI -12EI 6EI 6EI
Ą#ń# Ą# ń#
ó# ó# Ą#
l3 l3 l2 Ą# l2
ó#Ą# ó# Ą#
-12EI 12EI -6EI -6EI 6EI 4EI
Ą# -6EI 2EI
ń# Ą# ń#
ó#Ą# ó# Ą#
K11 = K12 = K21 = K22 =
ó# Ą# ó# Ą#
ó# ó# Ą#
l3 l3 l2 Ą# l2 l2 l l
l2
Ł# Ś# Ł# Ś#
ó#Ą# ó# Ą#
6EI -6EI 4EI 2EI
ó#Ą# ó# Ą#
ó#Ą# ó# Ą#
l2 l2 l l
Ł#Ś# Ł# Ś#
12EI -12EI 6EI 6EI 3EI -3EI 3EI
Ą#ń# Ą# ń# Ą# ń#
ó#
l3 l3 l2 Ą# ó# l2 Ą# ó# l3 l3 l2 Ą#
ó#Ą# ó# Ą# ó# Ą#
l 6EI
Ą# ń# Ą# -6EI 2EI
ń#
ó#-12EI 12EI -6EI Ą# ó#-6EI Ą# ó#-3EI 3EI -3EI Ą#
K ' =- =
Ą# ó# Ą#
ó# ó#
l3 l3 l2 Ą# ó# l2 Ą# ó# 4EI l2 l2 ll3 l3 l2 Ą#
Ł# Ś# Ł# Ś#
ó#Ą# ó# Ą# ó# Ą#
6EI -6EI 4EI 2EI 3EI -3EI 3EI
ó#Ą# ó# Ą# ó# Ą#
ó#Ą# ó# Ą# ó# Ą#
l2 l2 l l l2 l2 l
Ł#Ś# Ł# Ś# Ł# Ś#
f
k
v
k
v
i
________________________________________________________________________________________________
2007-02-26 96
Dynamika Budowli Wydział Inżynierii Lądowej i Środowiska
Magdalena Rucka & Krzysztof Wilde
Macierz bezwładności elementu belkowego:
f
k
f
i
v
k
v
i
156 22L 54 -13L
Ą# ń#
ó#
22L 4L2 13L -3L2 Ą#
źL
ó# Ą#
M = (8)
ó# Ą#
420 54 13L 156 -22L
ó#
-22L 4L2 Ą#
Ł#-13L -3L2 Ś#
gdzie ź jest masą elementu na jednostkę długości.
________________________________________________________________________________________________
2007-02-26 97
Wyszukiwarka
Podobne podstrony:
INSTALACJE ELEKTRYCZNE skrypt PG 2004Dynamika Budowli wyklad 4 11 12Prawo budowlane skrypt wersja 1Dynamika Budowli wyklad 3 11 12Skrypt z prawa budowlanegoMechanika Techniczna I Skrypt 1 7 1 Przedmiot dynamikiEntropia w Układach Dynamicznych II Downarowicz skrypt p16dynamika cieplna przegrod budowlanych8 37 Skrypty w Visual Studio (2)Rys budowlany 7MATLAB cw Skryptysyst oper skrypty 2Skrypt Latexwięcej podobnych podstron