Strona 1 z 6
Opracował: mgr inż. Jerzy Wołoszyn (
jwoloszy@agh.edu.pl
) mgr inż. Michał Ryś (
mrys@agh.edu.pl
)
Instrukcja do Laboratorium 2
W zadaniu 1 (rys. 1) należy wyznaczyć:
częstości drgań własnych oraz przemieszczenie pręta,
narysować wykres błędu względnego w zależności od liczby elementów (dla
elementu prętowego liniowego, parabolicznego i sześciennego)
Rys. 1
Dane do zadania :
p=1000 [N]
– siła ,
E=2*10e11 [Pa]
– moduł Younga,
L=2 [m] -
długość pręta,
ro=7850 [kg/m3]
– gęstość materiału,
A=0.003 [m2]
– przekrój poprzeczny elementu.
Rozwiązanie :
Zadanie zostanie rozwiązane przy wykorzystaniu elementu skończonego typu pręt:
a) o liniowej funkcji kształtu,
b) o parabolicznej fukncji kształtu,
c) o sześciennej fukcji kształtu.
Kod Programu (Matlab)
Opis i interpretacja Graficzna
p = -1000;
E = 2e11;
L = 2;
ro = 7850;
A = 0.003;
1. Wprowadzenie danych materia
łowych oraz
wymiarów.
el = 1
2. Wybór ilości elementów skończonych.
a)
Rozwiązanie z zastosowaniem elementu prętowego o liniowej funkcji kształtu.
n=el+1;
3
. Wyznaczenie liczby węzłów.
k = [1 -1;
-1 1 ];
m = [ 2 1;
1 2 ];
kk = ( (E*A) / ( L/(n-1) ) )*k;
4. Budowa macierzy bezwładności (mas) oraz
macierzy sztywności dla elementu prętowego
liniowego
w
jego
lokalnym
układzie
współrzędnych.
2
1
1
1
1
1
2
1
U
U
L
A
E
=
kk
U
U
Strona 2 z 6
Opracował: mgr inż. Jerzy Wołoszyn (
jwoloszy@agh.edu.pl
) mgr inż. Michał Ryś (
mrys@agh.edu.pl
)
mm=( ( ro*A* ( L/(n-1) ) ) / 6 )*m;
2
1
6
2
1
1
2
2
1
U
U
L
A
=
mm
U
U
%%%%%%%%% agregacja
K = zeros(n)
for i = 2:n
pk = i-1:i
K(pk,pk) = K(pk,pk) + kk;
end
M = zeros(n);
for i = 2:n
pm = i-1:i;
M(pm,pm) = M(pm,pm) + mm;
end
5. Budowa globalnej macierzy sztywności oraz
bezwładności (agregacja macierzy)
nm
nm
nm
nm
n
n
n
n
k
k
k
k
k
k
k
k
k
k
k
k
K
0
0
0
0
0
0
2
2
2
2
12
12
12
12
gdzie: n=3,
4,…; m=4,5,….
k
nm
– współczynnik z lokalnej macierzy
sztywności elementu o węzłach n,m
Budowa globalnej macierzy bezwładności jest
analogiczna. Przy budowie globalnej macierzy
bezwładności oraz sztywności sumują się tylko
elementy o wspólnych węzłach jak wyżej.
b)
Rozwiązanie z zastosowaniem elementu prętowego o parabolicznej funkcji kształtu.
n=2*el+1;
3. Wyznaczenie liczby węzłów.
k = [ 7 -8 1 ;
-8 16 -8 ;
1 -8 7 ]
m = [ 4 2 -1; 2 16 2 ; -1 2 4 ]
kk = ( ( E*A ) / ( (3*L) / el ) )*k;
mm = ( ( ro*A* ( L/el ) ) / 30 )*m;
4. Budowa ma
cierzy bezwładności (mas) oraz
macierzy sztywności dla elementu prętowego w
jego lokalnym układzie współrzędnych.
7
8
1
8
16
8
1
8
7
3 L
A
E
=
kk
4
2
1
2
16
2
1
2
4
30
L
A
=
mm
K = zeros(n);
for i = 3:2:n
pk = i-2:i;
K(pk,pk) = K(pk,pk) + kk
end
M = zeros(n);
for i = 3:2:n
pm = i-2:i;
M(pm,pm) = M(pm,pm) + mm;
end
5. Budowa globalnej macierzy sztywności
(agregacja macierzy)
Strona 3 z 6
Opracował: mgr inż. Jerzy Wołoszyn (
jwoloszy@agh.edu.pl
) mgr inż. Michał Ryś (
mrys@agh.edu.pl
)
c)
Rozwiązanie z zastosowaniem elementu prętowego o sześciennej funkcji kształtu.
n=3*el+1;
3. Wyznaczenie liczby węzłów.
k =[ 37/10 -189/40 27/20 -13/40 ;
-189/40 54/5 -297/40 27/20 ;
27/20 -297/40 54/5 -189/40 ;
-13/40 27/20 -189/40 37/10 ];
m=[[ 8/105 33/560 -3/140 7/619 ;
33/560 27/70 -27/560 -3/140 ;
-3/140 -27/560 27/70 33/560 ;
7/619 -3/140 33/560 8/105 ];
kk = ( E*A / (L/el) )*k;
mm = ( ro*A * (L/el) )*m;
4. Budowa macierzy bezwładności (mas) oraz
macierzy sztywności dla elementu prętowego w
jego lokalnym układzie współrzędnych.
10
37
40
189
20
27
40
13
40
189
5
54
40
297
20
27
20
27
40
297
5
54
40
189
40
13
20
27
40
189
10
37
L
A
E
=
kk
105
8
560
33
140
3
619
7
560
33
70
27
560
27
140
3
140
3
560
27
70
27
560
33
619
7
140
3
560
33
105
8
L
A
=
mm
K = zeros(n);
for i = 4:3:n
pk = i-3:i;
K(pk,pk) = K(pk,pk) + kk;
end
M = zeros(n);
for i = 4:3:n
pm = i-3:i;
M(pm,pm) = M(pm,pm) + mm;
end
5. Budowa globalnej macierzy sztywności
(agregacja macierzy)
P=zeros(1,n)';
P(n,1)=p
6. Definicja siły
K(:,1) = 0;
K(1,:) = 0;
K(1,1) = 1
M(:,1) = 0;
M(1,:) = 0;
M(1,1) = 1
7.
Odebranie
stopni
swobody
(brak
przemieszczeń na kierunku U w punkcie 1)
umes = inv(K)*P
omega2 = eig(K,M);
fmes = sqrt(omega2(2))/(2*pi)
8. Rozwiązanie numeryczne
u =( L / (A*E) )*p
omega = ( pi/(2*L) ) * sqrt(E/ro);
f = omega/(2*pi)
9. Rozwiązanie analityczne
Strona 4 z 6
Opracował: mgr inż. Jerzy Wołoszyn (
jwoloszy@agh.edu.pl
) mgr inż. Michał Ryś (
mrys@agh.edu.pl
)
W zadaniu 2
(rys. 2) należy wyznaczyć częstości drgań własnych oraz przemieszczenia.
Rys.2
Dane do zadania:
p = 1000 [N]
– siła ,
E = 2e11 [Pa]
– moduł Younga,
L = 2 [m] -
długość pręta,
ro = 7850 [kg/m3]
– gęstość materiału,
A = 0.003 [m2]
– przekrój poprzeczny elementu.
Rozwiązanie:
Zadanie zostanie rozwiązane przy wykorzystaniu elementu skończonego typu pręt o liniowej
funkcji kształtu.
Programy pomocnicze do obliczenia sin
– ów i cos – ów :
1.
policz_cosinusa.m
function cosinus = policz_cosinusa(x1,x2,y1,y2)
cosinus = (x2-x1)/sqrt((x2 - x1)^2 + (y2 - y1)^2);
2.
policz_sinusa.m
function sinus = policz_sinusa(x1,x2,y1,y2)
sinus = (y2-y1)/sqrt((x2 - x1)^2 + (y2 - y1)^2 );
Kod Programu (Matlab)
Interpretacja Graficzna
%============ dane:
p = 1000;
E = 2e11;
L = 2;
ro = 7850;
A = 0.003;
1. Wprowadzenie danych materiałowych oraz
wymiarów.
%============
współrzędne:
% 1 2 3
x = [0, 0.5*L, L];
y = [0, (sqrt(3)/2)*L, 0];
2.
Budowa
geometrii
kratownicy
(wprowadzenie współrzędnych)
Strona 5 z 6
Opracował: mgr inż. Jerzy Wołoszyn (
jwoloszy@agh.edu.pl
) mgr inż. Michał Ryś (
mrys@agh.edu.pl
)
k=[1 0 -1 0 ;
0 0 0 0 ;
-1 0 1 0 ;
0 0 0 0] ;
m=[2 0 1 0 ;
0 2 0 1 ;
1 0 2 0 ;
0 1 0 2] ;
%============
k w ukł.lokalnym:
k12 = ( ( E*A )/L )*k;
k23 = ( ( E*A )/L )*k;
m12 = ( ( ro*A*L )/6 )*m;
m23 = ( ( ro*A*L )/6 )*m;
3. Zdefiniowanie macierzy sztywności i
macierzy
bezwładności
dla
elementu
prętowego
w jego lokalnym
układzie
współrzędnych
0
0
0
0
0
1
0
1
0
0
0
0
0
1
0
1
L
A
E
=
k
2
0
1
0
0
2
0
1
1
0
2
0
0
1
0
2
6
L
A
=
m
c12 = policz_cosinusa(x(1),x(2),y(1),y(2))
s12 = policz_sinusa(x(1),x(2),y(1),y(2))
c23 = policz_cosinusa(x(2),x(3),y(2),y(3))
s23 = policz_sinusa(x(2),x(3),y(2),y(3))
DC12=[ c12 s12 0 0 ;
-s12 c12 0 0 ;
0 0 c12 s12 ;
0 0 -s12 c12 ];
DC23=[ c23 s23 0 0 ;
-s23 c23 0 0 ;
0 0 c23 s23 ;
0 0 -s23 c23 ];
%============
k w ukł.globalnym:
ko12 = DC12' * k12 * DC12;
ko23 = DC23' * k23 * DC23;
4. Transformacja z lokalnego układu
współrzędnych
do
układu
globalnego,
macierz
bezładności
nie
podlega
transformacji
c
12
= cos β
1
s
12
= sin β
1
c
23
= cos β
2
s
23
= sin β
2
c
s
s
c
c
s
s
c
DC
0
0
0
0
0
0
0
0
DC
k
DC
=
k
T
o
%============ agregacja
Kg = zeros(6);
gdzie = [1 2 3 4];
Kg(gdzie, gdzie) = Kg(gdzie, gdzie) + ko12;
gdzie =[3 4 5 6];
Kg(gdzie, gdzie) = Kg(gdzie, gdzie) + ko23;
5. Budowa globalnej macierzy sztywności
i
bezwładności (agregacja macierzy)
3
3
2
2
1
1
0
0
0
0
0
0
0
0
0
23
0
23
0
23
0
23
0
23
0
23
0
23
0
23
0
23
0
23
0
23
0
12
0
23
0
12
0
12
0
12
0
23
0
23
0
23
0
12
0
23
0
12
0
12
0
12
0
12
0
12
0
12
0
12
3
3
2
0
12
2
0
12
1
0
12
1
0
12
V
U
V
U
V
U
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
=
k
V
U
V
U
V
U
o
g
Strona 6 z 6
Opracował: mgr inż. Jerzy Wołoszyn (
jwoloszy@agh.edu.pl
) mgr inż. Michał Ryś (
mrys@agh.edu.pl
)
Mg = zeros(6);
gdzie =[1 2 3 4];
Mg(gdzie, gdzie) = Mg(gdzie, gdzie) + m12;
gdzie =[3 4 5 6];
Mg(gdzie, gdzie) = Mg(gdzie, gdzie) + m23;
3
3
2
2
1
1
0
0
0
0
0
0
0
0
23
23
23
23
23
23
23
23
23
23
23
12
23
12
12
12
23
23
23
12
23
12
12
12
12
12
12
12
3
3
2
12
2
12
1
12
1
12
V
U
V
U
V
U
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
=
m
V
U
V
U
V
U
g
%============ odebranie stopni swobody
Kg(:,[1 2 5 6])=0;
Kg([1 2 5 6],:)=0;
Kg(1,1)=1;
Kg(2,2)=1;
Kg(5,5)=1;
Kg(6,6)=1;
Kg
Mg(:,[1 2 5 6])=0;
Mg([1 2 5 6],:)=0;
Mg(1,1)=1;
Mg(2,2)=1;
Mg(5,5)=1;
Mg(6,6)=1;
Mg
6. Odebranie stopni swobody
3
3
2
2
1
1
1
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
1
0
23
0
12
0
23
0
12
0
23
0
12
0
23
0
12
3
3
2
2
1
1
V
U
V
U
V
U
k
k
k
k
k
k
k
k
=
k
V
U
V
U
V
U
o
g
- analogicznie dla macierzy
bezwładność
u = inv(Kg) * [0;0;0;p;0;0]
omega2 = eig(Kg,Mg);
fmes = sqrt( omega2(2) ) / (2*pi)
7. Wyzn
aczenie przemieszczeń węzłowych
Zadanie 3
(rys.3) do samodzielnego rozwiązania.
Wyznacz przemieszczenia oraz częstości drgań własnych przyjmując element skończony jako
element prętowy liniowy.
Dane do zadania:
p=3 [kN]
– siła ,
E=2e11 [Pa]
– moduł Younga,
L=2 [m] -
długość pręta,
ro=7850 [kg/m3]
– gęstość materiału,
A=0.003 [m2]
– przekrój poprzeczny elementu.