MK-8
STATYKA USTROJÓW PRĘTOWYCH W UJĘCIU MES
Dla kratownicy przedstawionej na rysunku wyznaczyć przemieszczenia węzłów i siły w prętach
wywołane ciężarem własnym konstrukcji.
> with(LinearAlgebra):
> m:=4; n:=14; e:=11; w:=7;
# m - wymiar wektora funkcji kształtu, n - liczba
współrzędnych globalnych całej struktury, e - liczba elementów , w - liczba węzłów kratownicy
> interface(rtablesize=n):
wektor funkcji kształtu dla przemieszczeń osiowych (pu) i poprzecznych (pv)
> pu:=Vector([1-x/h,0,x/h,0]); pv:=Vector([0,1-x/h,0,x/h]);
# h - długość
elementu
Wyznaczenie macierzy sztywności K, macierzy obrotu R i wektora sił f typowego elementu prętowego
> puprim:=map(diff,pu,x);
> K:=map(int,E*A*puprim.Transpose(puprim),x=0..h);
> R:=Matrix([[cos(alpha),sin(alpha),0,0],[-sin(alpha),cos(alpha),0,0],
[0,0,cos(alpha),sin(alpha)],[0,0,-sin(alpha),cos(alpha)]]);
> f:=-q0*map(int,sin(alpha)*pu+cos(alpha)*pv,x=0..h);
Wprowadzenie danych opisujących poszczególne elementy – uzupełnić na podstawie rysunku
> dane:=table([1=[0,l],...
# tablica danych o elementach:. numer_elementu = [kąt, długość]
> W:=Matrix([[1,2,3,4],...
# macierz wskaźników współrzędnych węzłowych
Wyznaczenie macierzy sztywności, macierzy obrotu i wektorów sił poszczególnych elementów
> for i to e do
K||i:=eval(K,h=dane[i][2]);
R||i:=eval(R,alpha=dane[i][1]);
f||i:=eval(f,[alpha=dane[i][1],h=dane[i][2]]);
end do:
∫
+
−
=
⌡
⌠
∂
∂
∂
∂
=
h
v
u
h
T
u
u
dx
q
dx
x
x
EA
0
0
0
)
cos
sin
(
p
p
f
p
p
K
α
α
1
s
q
2
s
q
3
s
q
4
s
q
5
s
q
6
s
q
7
s
q
8
s
q
9
s
q
10
s
q
11
s
q
12
s
q
13
s
q
14
s
q
I
II
III
IV
V
VI
VII
1
2
3
4
5
6
7
8
9
10
11
Obrót macierzy sztywności i wektorów sił poszczególnych elementów do układu globalnego
> for i to e do
Kg||i:=Transpose(R||i).K||i.R||i;
fg||i:=Transpose(R||i).f||i;
end do:
Złożenie macierzy sztywności i wektora sił całej struktury
> for i to e do B||i:=Matrix(4,n) end do:
> for i to e do
for j to m do B||i[j,W[i,j]]:=1 end do;
# - macierze logiczne
end do;
> Ks:=eval(add(Transpose(B||i).Kg||i.B||i,i=1..e)):
# macierz sztywności
całej struktury
> fs:=eval(add(Transpose(B||i).fg||i,i=1..e));
# wektor sił całej struktury
Redukcja macierzy sztywności, wektora sił i wektora przemieszczeń
> Kr:=SubMatrix(Ks,[3..7,9..n],[3..7,9..n]);
# zredukowana macierz sztywności
> fr:=SubVector(fs,[3..7,9..n]);
# zredukowany wektor sił
> qs:=Vector(n,symbol=q):
# wektor przemieszczeń węzłowych całej struktury
> qr:=SubVector(qs,[3..7,9..n]);
# zredukowany wektor przemieszczeń węzłowych
Rozwiązanie układu równań: Kr.wr = fr
> roz:=solve(convert(Kr.qr-fr,set),convert(qr,set));
> assign(roz);
> q[1],q[2],q[8]:=0,0,0;
# więzy nałożone na układ
> qs:=Vector([seq(q[i],i=1..n)]):
# wektor wyliczonych przemieszczeń węzłowych
Wyliczenie sił w poszczególnych prętach (rozwiązanie symboliczne)
> for i to e do
f||i:=simplify(K||i.R||i.B||i.qs):
end do:
Dane liczbowe
> A:=0.001;g:=9.81;rho:=8000;q0:=rho*g*A;E:=2e11;l:=1;
Wartości sił w poszczególnych prętach
> for i to e do evalf(f||i[3]); end do;
Wizualizacja wyników
> qp:=evalf(10000*qs):
# przeskalowany wektor przemieszczeń węzłowych
Geometria kratownicy (współrzędne poszczególnych węzłów)
> P1:=[0,0];P2:=[l,0];P3:=[2*l,0];P4:=[3*l,0];P5:=[l/2,l*sqrt(3)/2];
P6:=[3/2,l*sqrt(3)/2];P7:=[5/2,l*sqrt(3)/2];
f
R
f
KR
R
K
T
g
T
g
=
=
∑
∑
=
=
=
=
e
i
g
T
i
s
e
i
i
g
T
i
s
i
i
1
1
f
B
f
B
K
B
K
s
i
i
i
i
q
B
R
K
f
=
Rysunek kratownicy przed obciążeniem
> q1:=plot([[P1,P2],[P2,P3],[P3,P4],[P1,P5],[P2,P5],[P2,P6],[P3,P6],
[P3,P7],[P5,P6],[P6,P7],[P4,P7]],color=red,axes=none,scaling=constrain
ed):
Dodanie przemieszczeń poszczególnych węzłów do ich współrzędnych
> j:=0:
for i to w do
P||i[1]:=P||i[1]+qp[i+j]:
P||i[2]:=P||i[2]+qp[i+1+j]:
j:=j+1:
end do:
Rysunek kratownicy po obciążeniu
> q2:=plot([[P1,P2],[P2,P3],[P3,P4],[P1,P5],[P2,P5],[P2,P6],
[P3,P6],[P3,P7],[P5,P6],[P6,P7],[P4,P7]],color=blue,axes=none,
scaling=constrained):
> plots[display]({q1,q2});
> unassign('q');
Zadanie
Wyliczyć przemieszczenia węzłów i siły w prętach kratownicy zamocowanej na obu podporach
nieprzesuwnych. Kratownica obciążona jest ciężarem własnym (wartość q0 jak w poprzednim zadaniu)
oraz siłą skupioną o wartości P = 1000 N, przyłożoną do węzła II.
P
1
s
q
2
s
q
3
s
q
4
s
q
5
s
q
6
s
q
7
s
q
8
s
q
9
s
q
10
s
q
11
s
q
12
s
q
13
s
q
14
s
q
I
II
III
IV
V
VI
VII
1
2
3
4
5
6
7
8
9
10
11