I.Obliczenia wykonano z wykorzystaniem programów MATHCAD14 i RM-WIN.

1. Dyskretyzacja układu.

Przyjęty przekrój: HEB200; STAL-S275JO

L = 3,6m

P= 90kN

q = 45kN/m

1. Uwolnione przemieszczenia na węzłach: u , v , φ,

2. Podział elementów.

   1. Element 1.

1. Schemat:

1. Wykres momentów M[kn/m]:

1. Wykres sił tnących T [kN]:

1. Element 2:

1. Schemat:

1. Wykres momentów:

1. Wykres sił tnących:

3.3 Element 3:

a) schemat:

b) wykresy momentów M[kN/m]:

c) wykresy sił tnących T [kN]:

3.4 Element 4 :

a) Schemat:

b) Wykres momentów:

c) Wykres sił tnących:

4. Przyjmuję następujące dane:

A= 78,1 mm²

I_y= 5696 mm⁴*10⁴

I_x=2003 mm⁴*10⁴

E = 2,1*10^-5 MPa,

5. Określenie współrzędnych:

WĘZEŁ	WSPÓŁRZĘDNE
	X
1	0
2	10,8
3	3,6
4	10,8
5	10,8

6. Dane geometryczne i cosinusy kierunkowe elementów.

ELEMENT	GEOMETRIA ELEMENTÓW	DŁUGOŚĆ	COSINUSY KIERUNKOWE
	PRZEKRÓJ POPRZECZNY	MOMENT BEZWŁADNOŚCI
1	0,078m	5,69*10^-3m^4	10,8
2	0,078m	5,69*10^-3m^4	7,2
3	0,078m	5,69*10^-3m^4	7,2
4	0,078m	5,69*10^-3m^4	3,6

7. Wektory równoważnych obciążeń węzłowych poszczególnych elementów w lokalnym układzie współrzędnych.

N 0 0 0 0

T 66,667 0 162,0 0

M -144,0 0 -194,4 0

0 0 0 0

23,222 0 162,0 0

72,0 0 194,4 0

1 2 3 4

8. Wektory równoważnych obciążeń węzłowych poszczególnych elementów w lokalnym układzie współrzędnych.

R_p^e = Q R_p^e

R_p^e = Q^-1 R_p^e = Q^T R_p^e

cosα sinα 0 0 0 0

-sinα cosα 0 0 0 0

0 0 1 0 0 0

Q= 0 0 0 cosα sinα 0

0 0 0 -sinα cosα 0

0 0 0 0 0 1

66,667 0 162,0

0 1 0 2 0 3

R_P ¹= -144 R_p²= R_p⁴= 0 R_P³= -194,4

23,333 0 162,0

0 2 0 3 0 4

72 0 194,4

66,667

0 1

-144,0

23,333

0 2

72,0

R= $\sum_{e}^{}R$_P^e = 162,0

0 3

-194,4

162,0

0 4

194,4

0

0 5

0

Macierz sztywności elementu.

$\frac{AL^{2}}{J}$ 0 0 $\ \ \ \ - \frac{AL^{2}}{J}$ 0 0

0 12 6L 0 -12 6L

K^e= $\frac{\text{EJ}}{L^{3}}$ 0 6L 4L² 0 -6L 2L²

- $\frac{AL^{2}}{J}$ 0 0 $\frac{AL^{2}}{J}$ 0 0

0 -12 -6L 0 12 -6L

0 6L 2L² 0 -6L 4L

K^e= Q^TK^eQ

Globalna macierz sztywności.

K= $\sum_{e = 1}^{n}K$_e

	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
1	1.54	0	0	1.54	0	0	0	0	0	0	0	0	0	0	0
2	0	12	63.6	0	-12	63.6	0	0	0	0	0	0	0	0	0
3	0	63.6	449.44	0	-63.6	213.27	0	0	0	0	0	0	0	0	0
4	-1.54	0	0	0.71	0	0	0.71	0	0	0	0	0	0	0	0
5	0	-12	-63.6	0	12	63.6	0	-12	43.2	0	0	0	0	0	0
6	0	63.6	213.27	0	63.6	43.2	0	43.2	103.68	0	0	0	0	0	0
7	0	0	0	0.71	0	0	0.71	0	0	0.71	0	0	0	0	0
8	0	0	0	0	-12	-43.2	0	12	43.2	0	-12	43.2	0	0	0
9	0	0	0	0	43.2	103.68	0	43.2	26.4	0	-43.2	103.68	0	0	0
10	0	0	0	0	0	0	0.71	0	0	0.17	0	0	-0.17	0	0
11	0	0	0	0	0	0	0	-12	43.2	0	12	21.6	0	-12	21.6
12	0	0	0	0	0	0	0	43.2	103.68	0	21.6	51.84	0	-21.6	25.92
13	0	0	0	0	0	0	0	0	0	0.17	0	0	0.17	0	0
14	0	0	0	0	0	0	0	0	0	0	-12	-21.6	0	12	-21.6
15	0	0	0	0	0	0	0	0	0	0	21.6	25.92	0	-21.6	14.4

Wyniki sumowania macierzy :

0,71 0 0 0,71 0 0

0 12 63,6 0 -12 43,2

K_1,2= 0 63,6 43,2 0 43,2 103,68 *10^-6

0,71 0 0 0,71 0 0

0 -12 - 43,2 0 12 43,2

0 43,2 103,68 0 43,2 26,4

0,71 0 0 0,71 0 0 u₄ 1,316*10^-4m

0 12 43,2 0 -12 43,2 u₅ -5,288*10^-3m

K*u= R* 0 43,2 207,36 0 43,2 103,68 *10^-6 * φ₆ = -5,225*10^-3 rad

0,71 0 0 0,71 0 0 u₇ 1,3616*10^-4m

0 -12 - 43,2 0 12 43,2 u₈ -0,012m

0 43,2 103,68 0 43,2 207,36 φ₉ -0,033 rad

Siły w globalnym układzie współrzędnych:

F^e= R_P^e+K^e-u

69,851 20,149 -150,967 34,535

15,125 15,125 -19,409 -171,116

F₁= -156,552 F₂ = 50,159 F₃= -143,426 F₄= -84,680

-20,149 15,125 173,033 34,535

-15,125 -20,149 -19,409 -171,116

-50,159 58,742 -222,859 39,644

II. Dyskretyzacja układu z wykorzystaniem programu Matlab.

Kod programu MATLAB:

format short;

E=2.1e8;

I=0.005696

A=0.00781;

l=3.6;

P=90;

q=45;

le=5; %liczba prętów

N=18; %liczba elementów do obliczenia

lw=le+1; %liczba węzłów

ien=[1 2;2 3;5 4;6 4;4 3]'

Ee=[E,E,E,E,E];

Ie=[I,I,2*I,2*I,2*I];

Ae=[A,A,A,A,A];

xxyy=[0 3.6;3.6 10.8;10.8 10.8;3.6 3.6;10.8 3.6;10.8 0]; % wektor poszczególnych współrzędnych

% obliczenie długości elementów

e=1; L1=3.6

e=2; L2=7.2

e=3; L3=7.2

e=4; L4=7.2

e=5; L5=3.6

lE=[L1,L2,L3,L4,L5]

e=1; k1r=sztywelramylok(e,Ee,Ie,Ae,lE) ;

e=2; k2r=sztywelramylok(e,Ee,Ie,Ae,lE) ;

e=3; k3r=sztywelramylok(e,Ee,Ie,Ae,lE) ;

e=4; k4r=sztywelramylok(e,Ee,Ie,Ae,lE) ;

e=5; k5r=sztywelramylok(e,Ee,Ie,Ae,lE) ;

%obliczenie macierzy sztywności elementów w układzie globalnym

e=1; alfa1=0;

theta1=macierztransform(alfa1)

k1rg=theta1'*k1r*theta1

e=2; alfa2=0;

theta2=macierztransform(alfa2)

k2rg=theta2'*k2r*theta2

e=3; alfa3=90;

theta3=macierztransform(alfa3)

k3rg=theta3'*k3r*theta3

e=4; alfa4=0;

theta4=macierztransform(alfa4)

k4rg=theta4'*k4r*theta4

e=5; alfa5=90;

theta5=macierztransform(alfa5)

k5rg=theta5'*k5r*theta5

% agregacja macierzy sztywności elem. krg1-krg3 do macierzy globalnej K

K=zeros(18,18);

%lm=pomocagregacji(2,3,lw,ien,le)

lm=[1 4 13 10 7;2 5 14 11 8;3 6 15 12 9;4 7 10 7 16;5 8 11 8 17;6 9 12 9 18]

e=1; K=agregacjamacglob(K,e,k1rg,3,2,le,lw,lm) ;

e=2; K=agregacjamacglob(K,e,k2rg,3,2,le,lw,lm) ;

e=3; K=agregacjamacglob(K,e,k2rg,3,2,le,lw,lm) ;

e=4; K=agregacjamacglob(K,e,k3rg,3,2,le,lw,lm) ;

e=5; K=agregacjamacglob(K,e,k5rg,3,2,le,lw,lm) ;

% obliczenie wektora sił węzłowych układu dla siły równomiernie rozłożonej

fr4=[0,-q*lE(4)/2,-q*lE(4)*lE(4)/12,0,-q*lE(4)/2,q*lE(4)*lE(4)/12]'

fr4g=theta4'*fr4

V=zeros(N,1) ;

e=4; V=agregacjawekobc(V,e,3,2,lm,fr4g)

%obliczanie wektora sił wezłowych od siły skupionej

VP=[0;0;0;0;-90;0;0;0;0;0;0;0;0;0;0;0;0;0]

V=V+VP

%uwzględnienie warunków brzegowych u1=w1=FI1=u3=w3=FI3=u4=w4=FI4=0

K(1,:)=[]

K(:,1)=[]

K(1,:)=[]

K(:,1)=[]

K(11,:)=[]

K(:,11)=[]

K(11,:)=[]

K(:,11)=[]

K(11,:)=[]

K(:,11)=[]

K(12,:)=[]

K(:,12)=[]

V(1)=[];

V(1)=[];

V(11)=[];

V(11)=[];

V(11)=[]

V(12)=[];

%Rozwiązanie układu równań KU=F

U=K\V

%obliczenie sił wewnętrznych

%element 1

u1g=[0;0;U(1);U(2);U(3);U(4)]

F1g=k1rg*u1g

%element 2

u2g=[U(2);U(3);U(4);U(5);U(6);U(7)]

F2g=k2rg*u2g

%element 3

u3g=[U(5);U(6);U(7);U(8);U(9);U(10)]

F3g=k3rg*u3g

%element 4

u4g=[0;0;0;U(8);U(9);U(10);]

F4g=-fr4g+k4rg*u4g

%element 5

u5g=[U(10);U(11);U(12);0;0;0]

F5g=k5rg*u5g

ex1=[xxyy(ien(1,1),1) xxyy(ien(2,1),1)]

ex2=[xxyy(ien(1,2),1) xxyy(ien(2,2),1)]

ex3=[xxyy(ien(1,3),1) xxyy(ien(2,3),1)]

ex4=[xxyy(ien(1,4),1) xxyy(ien(2,4),1)]

ex5=[xxyy(ien(1,5),1) xxyy(ien(2,5),1)]

ey1=[xxyy(ien(1,1),2) xxyy(ien(2,1),2)]

ey2=[xxyy(ien(1,2),2) xxyy(ien(2,2),2)]

ey3=[xxyy(ien(1,3),2) xxyy(ien(2,3),2)]

ey4=[xxyy(ien(1,4),2) xxyy(ien(2,4),2)]

ey5=[xxyy(ien(1,5),2) xxyy(ien(2,5),2)]

eq1=[0 0]

eq2=[0 0]

eq3=[0 0]

eq4=[0 0]

eq5=[0 0]

format long

ep1=[Ee(1) Ae(1) Ie(1)];

ep2=[Ee(2) Ae(2) Ie(2)];

ep3=[Ee(3) Ae(3) Ie(3)];

ep4=[Ee(4) Ae(4) Ie(5)];

ep5=[Ee(5) Ae(5) Ie(5)];

es1=beam2s(ex1,ey1,ep1,u1g ’,eq1,10);

es2=beam2s(ex2,ey2,ep2,u2g’,eq2,10);

es3=beam2s(ex3,ey3,ep3,u3g’,eg3,10);

es4=beam2s(ex4,ey4,ep4,u4g’,eg4,10);

es5=beam2s(ex5,ey5,ep5,u5g’,eg5,10);

plotpar=[2 1];

sfac=scalfact2(ex5,ey5,es5(:,3),0.2)

clf;

eldia2(ex1,ey1,es1(:,3),plotpar,sfac);tilte(‘ELEMENT 1’)

eldia2(ex2,ey2,es2(:,3),plotpar,sfac);tilte(‘ELEMENT 2’)

eldia2(ex3,ey3,es3(:,3),plotpar,sfac);tilte(‘ELEMENT 1 2 3’)

eldia2(ex4,ey4,es4(:,3),plotpar,sfac);tilte(‘ELEMENT 1 2 3 4’)

eldia2(ex5,ey5,es5(:,3),plotpar,sfac);tilte(‘MOMENTY ZGINAJACE-ELEMENT 1 2 3 4 5’);

Niezbędne dodatki do programu:

*agregacjamacglob.m

function K=agregacjamacglob(K,e,ke,ssw,lwe,le,lw,lm)

for m=1:lwe*ssw

im=lm(m,e);

s=im;

if s~=0

for n=1:lwe*ssw

jm=lm(n,e);

r=jm;

if r~=0

K(im,jm)=K(im,jm)+ke(m,n);

end

end

end

end

end

*agregacjawekobc.m

function V=agregacjawekobc(V,e,ssw,lwe,lm,fre)

for n=1:ssw*lwe

j=lm(n,e);

ss=j;

if ss~=0

V(j,1)=V(j,1)+fre(n);

end

end

*dlugoscelem.m

function y=dlugoscelem(xxyy,e,ien)

a=ien(2,e);

b=ien(1,e);

xx=xxyy(a,1)-xxyy(b,1);

yy=xxyy(a,2)-xxyy(b,2);

y=sqrt(xx*xx+yy*yy);

*macierztransform.m

function thetae=macierztransform(alfa)

thetae=zeros(6,6);

x=alfa*pi/180;

C=cos(x);

S=sin(x);

thetae=[C,S,0,0,0,0;-S,C,0,0,0,0;0,0,1,0,0,0;0,0,0,C,S,0;0,0,0,-S,C,0;0,0,0,0,0,1];

*sztywelramylok.m

function ker=sztywelramylok(e,Ee,Ie,Ae,lE)

ker=zeros(6,6);

w=Ae(1,e)*lE(1,e)*lE(1,e)/Ie(1,e);

ker=Ee(1,e)*Ie(1,e)/(lE(1,e)*lE(1,e)*lE(1,e))*[w,0,0,-w,0,0;0,12,6*lE(1,e),0,-12,6*lE(1,e);0,6*lE(1,e),4*lE(1,e)*lE(1,e),0,-6*lE(1,e),2*lE(1,e)*lE(1,e);-w,0,0,w,0,0;0,-12,-6*lE(1,e),0,12,-6*lE(1,e);0,6*lE(1,e),2*lE(1,e)*lE(1,e),0,-6*lE(1,e),4*lE(1,e)*lE(1,e)];

Otrzymane wyniki w programie RM-WIN:

Wykresy:

M[kN/m]

T [kN]

N [kN]

Podsumowanie i wnioski:

Zastosowanie obliczeń układu Metodą Elementów Skończonych umożliwia rozwiązanie przybliżone do wyników otrzymanych w programie obliczeniowym RM-WINi MATHCAD14. Maksymalny moment elementu obciążonego równomiernie uzyskany metodą MES wynosi 227,1946 podczas gdy program podaje wartość 222,859 co stanowi 2% różnicę. Podobnie otrzymujemy wyniki dla sił tnących i kątów obrotu na poziomie < 6%. Oznacza to poprawne rozwiązanie układu metodą MES.