Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
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, b1, h1, b2, h2, E, r�, x�1�, x�2�.� �
2000
1500
1000
500
0
0 0.02 0.04 0.06 0.08 0.1
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).
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 1
p [N]
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
Opis stosowanych funkcji:
[Ke]=ke_beam(EJ,L)
Funkcja generuje lokalną macierz sztywności elementu belkowego:
12EI 6EI -�12EI 6EI
��ł�
ę�
L3 L2 L3 L2 ś�
��
6EI 4EI -�6EI 2EI
��
��
L2 L L2 L
ke =�
��
ę�-�123EI -�6EI 12EI -�6EI ś�
L L2 L3 L2 ś�
ę�
��
6EI 2EI -�6EI 4EI
��
�� L2 L L2 L ��
%---------------------------------------------------------------
% 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
%----------------------------------------------------------------
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 2
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
[Me]=me_beam(mi,L)
Funkcja generuje lokalną macierz mas elementu belkowego:
156 22L 54 -�13L
��ł�
ę�
22L 4L2 13L -�3L2 ś�
m�L
��
me =�
420 ę� 54 13L 156 -�22Lś�
��
��-�13L -�3L2 -�22L 4L2 ��
%---------------------------------------------------------------
% 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.
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 3
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
[KG]=agreg_k_beam(Neq,Ne,CE,ALOK)
nel
i
Funkcja agregacji globalnej macierzy sztywności dla płaskiego elementu belkowego K =�
��k
e
i=�1
% 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)
nel
i
Funkcja agregacji globalnej macierzy mas dla płaskiego elementu belkowego M =�
��m
e
i=�1
% 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
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 4
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
[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)
%----------------------------------------------------------
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 5
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
Kolejność obliczeń:
a) Dokonać dyskretyzacji układu stosując element belkowy (z pełnym kompletem
więzów).
T
Wektor przemieszczeń węzłowych elementu belkowego: qe =� vi j�i vk j�k .
[� ]�
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 6
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
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]
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 7
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
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];
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 8
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
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)
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 9
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
g) Rozwiązać problem własny:
[mode,vale] = eig(K,M)
mode  macierz modalna
f�11 f�12 f�1N
��ł�
ę�f� f�22
f�2 N ś�
21
�� ł�ę�ś�
Ś =�=�
��f�jn ����
ę�f� f�N
f�NN ś�
�� N12 ��
vale  macierz widmowa
2
��ł�
w�1 00
��
2
0 w�2 0
��
�2 =�
��
��
2
0 0 w�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
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 10
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
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
%----------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 11
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
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
postac_1_2_3=mode(:,1:3);
f =
postac_1_2_3=
1.4602
0.2591 -0.4811 -1.0000
5.7043
0.2128 -0.3061 -0.4350
9.4556
0.7287 -0.7810 -0.5941
12.2021
0.2079 0.0639 0.7837
20.1370
0.0323 0.1869 0.3825
25.5384
0.0249 0.2311 0.3956
41.3254
0.0018 0.0878 0.0966
45.5184
0.0152 0.2481 0.3188
56.9852
-0.0108 -0.0719 -0.1678
82.4466
1.0000 -0.3908 0.9954
100.4504
-0.0052 -0.2343 -0.1725
114.2185
0.8966 0.5858 0.2416
148.3424
0.1123 -0.5023 0.7629
187.8369
0.5372 1.0000 -0.9283
257.4507
0.2287 0.1796 0.0939
424.6326 0.2760 0.6509 -0.7657
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 12
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
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
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 13
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
j) Utworzyć macierz tłumienia C
Macierz tłumienia C wyznaczamy jako macierz tłumienia proporcjonalnego (macierz
Rayleigha): C =� a0M +� a1K.
Współczynniki a0 i a1 macierzy tłumienia wyznaczamy z równania:
-�1
1
�� 1
��ł�
w�i ł�
��
ę�w� w�1 ś�
w�i x�i
a0 �� �� a0 x�1
�� ł� �� ł� �� ��
1
1
��
ś�
=�
�� ż�.
ę�a ś���x� ż�, ę�a ś� =� 2ę� 1
1
2 ��
��
j
�� 1 �� �� 1 �� ��x�2 ��
�� ��
ę�w� w�j ś�
ę�w� w�2 ś�
j
�� 2 ��
����
jeżeli znane są dwie wartości liczby tłumienia x�i i x� , odpowiadające częstościom drgań
j
w�i oraz w� .
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);
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 14
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
l) Przyjąć dt Ł� dtkr, 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.
2000
1500
1000
500
0
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
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;
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 15
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
o) Zdefiniować wektory przemieszczeń i prędkości początkowych u0 oraz v0.
u0=zeros(1,Neq); % wektor przemieszczenia poczatkowego
v0=zeros(1,Neq); % wektor predkosci poczatkowej
p) Dokonać całkowania równań ruchu Mu +� Cu +� Ku =� p(t) metodą różnic centralnych:
[u,v,a]=mrc(M,C,K,P,t,u0,v0);
q) Wykreślić odpowiedz
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]')
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 16
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
2000
1000
0
0 1 2 3 4 5 6 7 8
-3
x 10
5
0
-5
0 1 2 3 4 5 6 7 8
0.1
0
-0.1
0 1 2 3 4 5 6 7 8
5
0
-5
0 1 2 3 4 5 6 7 8
t [s]
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 17
5
p [N]
12
u
[m]
12
v
[m/s]
2
12
a
[m/s ]
Dynamika Budowli  laboratorium Ćwiczenie 5 Magdalena Rucka
Zadanie do samodzielnego rozwiązania:
Politechnika Gdańska, Wydział Inżynierii Lądowej i Środowiska, Katedra Wytrzymałości Materiałów 18