Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
1
Wyznaczyć częstości i postacie drgań własnych ramy portalowej. Wykreślić postacie
drgań odpowiadające trzem pierwszym częstościom drgań własnych. Obliczyć
przemieszczenie u(t) wywołane obciążeniem p(t) korzystając z metody różnic centralnych.
Zastosować macierzową metodę przemieszczeń (dyskretyzacja na 9 elementów
belkowych). Dane: L, H, b
1
, h
1
, b
2
, h
2
, E,
,
,
0
0.02
0.04
0.06
0.08
0.1
0
500
1000
1500
2000
p
[
N
]
Rozwiązanie:
Napisać program realizujący zadanie w środowisku MATLAB. Jako bazy użyć gotowy
skrypt o nazwie cwiczenie_5_szablon.m. Program ten wykorzystuje następujące funkcje:
ke_beam, m_beam, agreg_k_beam, agreg_m_beam, mrc (funkcje należy skopiować do
katalogu roboczego).
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
2
Opis stosowanych funkcji:
[Ke]=ke_beam(EJ,L)
Funkcja generuje lokalną macierz sztywności elementu belkowego:
3
2
3
2
2
2
3
2
3
2
2
2
12
6
12
6
6
4
6
2
12
6
12
6
6
2
6
4
e
EI
EI
EI
EI
L
L
L
L
EI
EI
EI
EI
L
L
L
L
EI
EI
EI
EI
L
L
L
L
EI
EI
EI
EI
L
L
L
L
k
%---------------------------------------------------------------
% WEJSCIE:
% EJ = sztywność giętna EJ
% L = długość elementu
%----------------------------------------------------------------
% WYJSCIE:
% Ke = macierz sztywności 4x4 względem przemieszczeń:
% v_a,fi_a,v_b,fi_b
%----------------------------------------------------------------
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
3
[Me]=me_beam(mi,L)
Funkcja generuje lokalną macierz mas elementu belkowego:
2
2
2
2
156
22
54
13
22
4
13
3
54
13
156
22
420
13
3
22
4
e
L
L
L
L
L
L
L
L
L
L
L
L
L
m
%---------------------------------------------------------------
% WEJSCIE:
% mi = rozkład masy na długości elementu
% L = długość elementu
%----------------------------------------------------------------
% WYJSCIE:
% Me = macierz mas 4x4 względem przemieszczeń:
% v_a,fi_a,v_b,fi_b
%----------------------------------------------------------------
Uwaga: funkcje ke_beam oraz me_beam nie są używane bezpośrednio w skrypcie
cwiczenie_5_przyklad.m. Do funkcji ke_beam oraz me_beam odwołują się funkcje
agreg_k_beam oraz agreg_m_beam.
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
4
[KG]=agreg_k_beam(Neq,Ne,CE,ALOK)
Funkcja agregacji globalnej macierzy sztywności dla płaskiego elementu belkowego
1
el
n
i
e
i
K
k
% WEJSCIE:
% Neq - liczba niewiadomych uogólnionych przemieszczeń
% Ne - liczba elementów
% CE - tablica cech elementów
% ALOK - tablica wektorów alokacji
%----------------------------------------------------------------
% WYJSCIE:
% KG - globalna macierz sztywności
[MG]=agreg_m_beam(Neq,Ne,CE,ALOK)
Funkcja agregacji globalnej macierzy mas dla płaskiego elementu belkowego
1
el
n
i
e
i
M
m
% WEJSCIE:
% Neq - liczba niewiadomych uogólnionych przemieszczeń
% Ne - liczba elementów
% CE - tablica cech elementów
% ALOK - tablica wektorów alokacji
%----------------------------------------------------------------
% WYJSCIE:
% MG - globalna macierz mas
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
5
[u,v,a]=mrc(M,C,K,P,t,u0,v0)
Funkcja całkowania równań ruchu metodą różnic centralnych.
%---------------------------------------------------------
% WEJSCIE:
% M - macierz mas (n x n)
% C - macierz tłumienia (n x n)
% K - macierz sztywności (n x n)
% P - wektor obciążeń zewnętrznych (n x nt)
% t - wektor czasu (1 x nt)
% u0 - wektor przemieszczeń początkowych (1 x n)
% v0 - wektor prędkości początkowych (1 x n)
%----------------------------------------------------------
% WYJSCIE:
% u - wektor przemieszczeń (n x nt)
% v - wektor prędkości (n x nt)
% a - wektor przyspieszeń (n x nt)
%----------------------------------------------------------
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
6
Kolejność obliczeń:
a) Dokonać dyskretyzacji układu stosując element belkowy (z pełnym kompletem
więzów).
Wektor przemieszczeń węzłowych elementu belkowego:
T
e
i
i
k
k
v
v
q
.
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
7
b) W programie MATLAB zadeklarować liczbę elementów i liczbę niewiadomych
uogólnionych przemieszczeń:
Ne = 9;
% liczba elementow
Neq = 16;
% liczba niewiadomych uogolnionych przemieszczen
c) Zadeklarować zmienne dotyczące danych geometrycznych i materiałowych:
E=15*10^9;
% modul sprezystosci [N/m2]
H=6.2;
% wysokosc ramy [m]
L=4.8 ;
% rozpietosc ramy [m]
b1=0.1;
% szerokosc przekroju belki [m]
h1=0.2;
% wysokosc przekroju belki [m]
b2=0.1;
% szerokosc przekroju slupa [m]
h2=0.12;
% wysokosc przekroju slupa [m]
ro=2300;
% gestosc [kg/m3]
ksi_1=3/100;
% liczba tlumienia ksi 1
ksi_2=5/100;
% liczba tlumienia ksi 2
Ib=(b1*h1^3)/12;
% moment bezwladnosci belki [m4]
Ic=(b2*h2^3)/12;
% moment bezwladnosci slupa [m4]
mi_b=ro*b1*h1;
% masa rozlozona belki [kg/m]
mi_c=ro*b2*h2;
% masa rozlozona slupa [kg/m]
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
8
d) Uzupełnić tablicę cech elementów (tablica CE).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TABLICA CECH ELEMENTOW
%-------------------------------------------
% Ne EI L ro
%-------------------------------------------
CE = [ 1 E*Ic H/3 mi_c
2 E*Ic H/3 mi_c
3 E*Ic H/3 mi_c
4 E*Ib L/3 mi_b
5 E*Ib L/3 mi_b
6 E*Ib L/3 mi_b
7 E*Ic H/3 mi_c
8 E*Ic H/3 mi_c
9 E*Ic H/3 mi_c];
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
9
e) Uzupełnić tablicę wektorów alokacji (tablica ALOK).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TABLICA WEKTOROW ALOKACJI
%-------------------------------------------
% Ne | Vi | Fi | Vk | Fk
%-------------------------------------------
ALOK=[ 1 0 0 1 2
2 1 2 3 4
3 3 4 10 5
4 0 5 6 7
5 6 7 8 9
6 8 9 0 11
7 12 13 10 11
8 14 15 12 13
9 0 16 14 15];
f) Dokonać agregacji globalnych macierzy sztywności K i mas M :
[K]=agreg_k_beam(Neq,Ne,CE,ALOK)
[M]=agreg_m_beam(Neq,Ne,CE,ALOK)
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
10
g) Rozwiązać problem własny:
[mode,vale] = eig(K,M)
mode
– macierz modalna
11
12
1
21
22
2
1
2
N
N
jn
N
N
NN
Φ
vale
– macierz widmowa
2
1
2
2
2
2
0
0
0
0
0
0
N
Ω
omega=sqrt(diag(vale));
– częstości kołowe drgań własnych (ekstrakcja diagonali
macierzy widmowej)
f = omega/2/pi;
– częstotliwości drgań własnych
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
11
h) Sortowanie i normalizacja częstości i postaci drgań własnych.
[mode,f]=sort_norm(mode,f)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% funkcja sortujaca i normalizujaca czestosci
% i postaci drgan wlasnych
%----------------------------------------------------------
% WEJSCIE:
% mode - macierz modalna
% f - wektor czestotliwosci
%----------------------------------------------------------
% WYJSCIE:
% mode_s - posortowana macierz modalna
% f_s - posortowany wektor czestotliwosci
%----------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
12
i) Wykreślić postacie drgań odpowiadające trzem pierwszym częstościom drgań
własnych.
mode =
0.2591 -0.4811 -1.0000 -0.8532 0.3807 0.9551 -0.0377 -0.1408 -0.0942 0.0296 0.1079 0.1159 -0.0019 0.0267 -0.0337 0.0118
0.2128 -0.3061 -0.4350 -0.2603 -0.0710 -0.5533 0.1031 0.5251 1.0000 0.4361 0.9081 0.8072 -0.0106 0.1337 -0.1564 0.0516
0.7287 -0.7810 -0.5941 -0.1618 -0.3382 -1.0000 0.0391 0.0382 -0.1825 -0.1109 -0.1555 -0.0683 -0.0004 0.0394 -0.0688 0.0333
0.2079 0.0639 0.7837 0.7205 -0.2885 -0.0034 -0.0996 -0.6169 -0.8840 0.0536 0.7161 1.0000 -0.0206 0.3664 -0.5136 0.2014
0.0323 0.1869 0.3825 -0.6099 -0.1325 0.2758 0.7307 0.3748 0.0556 -0.2315 -0.6771 0.1770 0.0804 1.0000 -0.9958 0.8635
0.0249 0.2311 0.3956 -0.9640 -0.2587 0.0218 0.5723 0.1678 -0.1058 -0.0107 0.0401 0.0369 -0.0106 -0.2067 0.1199 0.0123
0.0018 0.0878 0.0966 -0.4183 -0.0995 -0.1657 -0.3562 -0.2306 -0.0106 0.2031 0.5860 -0.3402 -0.0951 -0.4932 -0.6805 0.9988
0.0152 0.2481 0.3188 -1.0000 -0.1896 -0.2034 -0.5111 -0.2524 0.0678 0.0420 -0.0403 -0.0086 0.0258 0.2034 0.1213 -0.0120
-0.0108 -0.0719 -0.1678 0.3851 0.1611 -0.0287 -0.4537 -0.0672 0.1123 -0.1908 -0.5799 0.3561 0.0524 -0.5169 0.6717 1.0000
1.0000 -0.3908 0.9954 0.3535 -0.3599 0.7718 0.3627 -0.0834 0.2207 0.1204 -0.1861 -0.0976 0.0516 0.1980 -0.0013 0.2037
-0.0052 -0.2343 -0.1725 0.7064 -0.0379 0.2341 0.8612 0.1372 -0.1387 0.2970 0.4880 -0.3175 0.0786 0.9844 1.0000 0.8664
0.8966 0.5858 0.2416 -0.1157 0.8791 -0.3112 -0.1854 0.2943 -0.1603 0.0824 0.0706 -0.0155 0.0711 0.0275 0.0638 0.0324
0.1123 -0.5023 0.7629 -0.1443 -0.1391 0.4316 -0.8617 0.6940 -0.0066 -0.9462 0.9330 -0.0789 -0.4217 0.4118 0.5447 0.2077
0.5372 1.0000 -0.9283 0.2921 -0.5475 0.1826 0.0979 -0.1553 0.1139 -0.2488 0.1945 -0.0169 0.0424 0.0109 0.0225 0.0089
0.2287 0.1796 0.0939 -0.0832 0.7717 -0.5302 0.9664 -0.9158 0.3331 -0.1693 -0.1823 0.0471 -0.8429 0.2785 0.2438 0.0708
0.2760 0.6509 -0.7657 0.2798 -1.0000 0.5780 -1.0000 1.0000 -0.4654 1.0000 -1.0000 0.1145 -1.0000 0.2397 0.1661 0.0409
f =
1.4602
5.7043
9.4556
12.2021
20.1370
25.5384
41.3254
45.5184
56.9852
82.4466
100.4504
114.2185
148.3424
187.8369
257.4507
424.6326
postac_1_2_3=mode(:,1:3);
postac_1_2_3=
0.2591 -0.4811 -1.0000
0.2128 -0.3061 -0.4350
0.7287 -0.7810 -0.5941
0.2079 0.0639 0.7837
0.0323 0.1869 0.3825
0.0249 0.2311 0.3956
0.0018 0.0878 0.0966
0.0152 0.2481 0.3188
-0.0108 -0.0719 -0.1678
1.0000 -0.3908 0.9954
-0.0052 -0.2343 -0.1725
0.8966 0.5858 0.2416
0.1123 -0.5023 0.7629
0.5372 1.0000 -0.9283
0.2287 0.1796 0.0939
0.2760 0.6509 -0.7657
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
13
Rysunek pierwszej postaci drgań własnych:
postac_1_2_3=mode(:,1:3);
postac_1_2_3=
0.2591 -0.4811 -1.0000
0.2128 -0.3061 -0.4350
0.7287 -0.7810 -0.5941
0.2079 0.0639 0.7837
0.0323 0.1869 0.3825
0.0249 0.2311 0.3956
0.0018 0.0878 0.0966
0.0152 0.2481 0.3188
-0.0108 -0.0719 -0.1678
1.0000 -0.3908 0.9954
-0.0052 -0.2343 -0.1725
0.8966 0.5858 0.2416
0.1123 -0.5023 0.7629
0.5372 1.0000 -0.9283
0.2287 0.1796 0.0939
0.2760 0.6509 -0.7657
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
14
j) Utworzyć macierz tłumienia C
Macierz tłumienia C wyznaczamy jako macierz tłumienia proporcjonalnego (macierz
Rayleigha):
0
1
a
a
C
M
K .
Współczynniki a
0
i a
1
macierzy tłumienia wyznaczamy z równania:
0
1
1
1
1
2
i
i
i
j
j
j
a
a
,
1
1
1
0
1
1
2
2
2
1
2
1
a
a
.
jeżeli znane są dwie wartości liczby tłumienia
i
i
j
, odpowiadające częstościom drgań
i
oraz
j
.
omega_1=f(1)*2*pi;
omega_2=f(2)*2*pi;
wspl=2*inv([1/omega_1 omega_1; 1/omega_2 omega_2])*[ksi_1;ksi_2];
a0=wspl(1);
a1=wspl(2);
C=a0*M + a1*K;
k) Wyznaczyć krytyczny krok całkowania dt z warunku stabilności metody różnic
centralnych:
dt_mrc=2/(f(end)*2*pi);
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
15
l) Przyjąć
kr
dt
dt
, zdefiniować wektor czasu:
t = [ 0:dt:tk];
Czas końcowy
tk
dobrać tak, by zaobserwować przejście układu w stan spoczynku.
m) Zdefiniować wektor obciążenia p(t), przyłożyć go wzdłuż wybranego stopnia
swobody.
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
0
500
1000
1500
2000
tc=0.1
% czas trwania impulsu [s]
po=2000;
% amplituda impulsu [N]
t1=[0:dt:tc];
p= po*sin(pi/tc*t1);
n) Utworzyć macierz P zawierającą obciążenia zewnętrzne wzdłuż wszystkich stopni
swobody:
P = zeros(Neq,length(t));
Na zadanym stopniu swobody (np. nr 5) w macierzy P umieścić wektor p
P(5,1:(length(p)))=p;
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
16
o) Zdefiniować wektory przemieszczeń i prędkości początkowych u
0
oraz v
0
.
u0=zeros(1,Neq);
% wektor przemieszczenia poczatkowego
v0=zeros(1,Neq);
% wektor predkosci poczatkowej
p) Dokonać całkowania równań ruchu
( )
t
Mu Cu Ku
p
metodą różnic centralnych:
[u,v,a]=mrc(M,C,K,P,t,u0,v0);
q) Wykreślić odpowiedź układu w czasie wzdłuż wybranych stopni swobody (np. nr 10).
figure(1)
plot(t,u(12,:),
'k'
);grid
on
;ylabel(
'u_1_0 [m]'
);xlabel(
't [s]'
)
figure(2)
plot(t,v(12,:),
'k'
);grid
on
;ylabel(
'v_1_0 [m/s]'
);xlabel(
't [s]'
)
figure(3)
plot(t,a(12,:),
'k'
);grid
on
;ylabel(
'a_1_0 [m/s^2]'
);xlabel(
't [s]'
)
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
17
0
1
2
3
4
5
6
7
8
0
1000
2000
p
5
[
N
]
0
1
2
3
4
5
6
7
8
-5
0
5
x 10
-3
u
1
2
[
m
]
0
1
2
3
4
5
6
7
8
-0.1
0
0.1
v
1
2
[
m
/s
]
0
1
2
3
4
5
6
7
8
-5
0
5
a
1
2
[
m
/s
2
]
t [s]
Dynamika Budowli
– laboratorium
Ćwiczenie 5
Magdalena Rucka
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów
18
Zadanie do samodzielnego rozwiązania: