Politechnika Krakowska

Wydział Inżynierii Lądowej

Instytut Technologii Informatycznej

METODY OBLICZENIOWE

PROJEKT NR II

SPRAWOZDANIE

Wykonał:

Radosław Burcek

gr. 10

Napisać w języku MATLAB program wyznaczający ugięcia belki metodą elementów skończonych. Narysować wykres sił przekrojowych.

Belka:

Tok postępowania:

1. Preprocessing
2. Processing

3. Postprocessing

Preprocessing polega na dyskretyzacji, czyli podziale belki na skończoną ilość elementów określonejdługości. Określenie liczby stopni swobodyw kazdym węźle, układzie i elemencie. Następnie generujemy macierz sztywności i wektor obciążeń węzłowych dla elementów. Macierz sztywnościgenerowana jest z wykorzystaniem sześciennych funkcji kształtu Hermite’a (funkcje wpływu). Ostatnim etapem preprocessingu jestokreśleniekinematycznych warunków brzegowych.

Processing to generowanieglobalnej macierzsztywności i wektora obciążenia przez agregacje wartości z każdego elementu. Następnym krokiem jest podstawienie warunków brzegowych i rozwiązanie układu równań KG*u=f, gdzie KG – globalna macierz sztywności, u – wektor przemieszczenia, f – wektor obciążenia (po jego rozwiązaniu otrzymujemy wartość reakcji i przemieszczenia.

Postprocessing, to generowanie ugięcia i macierzy sił przywęzłowych.

Model matematyczny

a) równanie różniczkowe IV rzędu

EI*(d⁴w/dx⁴) + q(x)=0

b) warunki brzegowe kinematyczne

w(0)=0

w’(0)=0

c) warunki brzegowe statyczne

M(L)=0

Q(L)=0

Przejście z modelu matematycznego w model dyskretny

Algorytm MES

1) SFORMUŁOWANIE PROBLEMU - PREPROCESSING

Wartości parametrów:

E= 200 GPa

I= 171 cm⁴

EI= (E*I)/100;

L= 4 m

q= 5 kN/m

Dyskretyzacja:

n=4

Nel = 4*n; ilość elementów

Nodes= Nel+1 ilość węzłów

hel=L/Nel; długość pojedynczego elementu

coord= linspace(0,L,Nodes) globalna macierz współrzędnych węzłów

node_dof= 2 liczba stopni swobody w węźle (ugięcie, kąt obrotu)

Ndofs= Nodes*node_dof liczba stopni swobody całego elementu

edof=[1 1 2 3 4] macierz stopni swobody w elemencie 1

Edof=[(1:Nel)',(1:2:Ndofs-3)',globalna macierz stopni swobody

(2:2:Ndofs-2)',(3:2:Ndofs-1)',

(4:2:Ndofs)']

Generowanie macierzy sztywności i wektora obciążenia węzłowego dla elementów

 Kel=double(subs(Kelsym,{h,EIs},{hel,EI}))   macierz sztywności poszczególnych elementów
 qel=double(subs(qelsym,h,hel)) obciążenie dla poszczególnych elementów

Q(1:n)=0; obciążenie dla zadanej belki

Q(n+1:2*n)=0;

Q(2*n+1:3*n)=0;

Q(3*n+1:4*n)=1;

kinematic_bc = kinematyczne warunki brzegowe dla zadanej belki

[1 0;2 0;Ndofs-1 0;]

2) PROCESSING

KG=zeros(Ndofs,Ndofs) globalna macierz sztywności

f=zeros(Ndofs,1) wektor sił przywęzłowych

for iel=1:Nel proces agregacji - pętla
     Qel = Q(iel)* qel;             
     edof=Edof(iel,:);                
     [KG,f]=assem(edof,KG,Kel,f,Qel);

 end

[u,r]=solveq(KG,f,kinematic_bc) wektor wartości stopni swobody

3) POSTPROCESSING

figure(1); rysowanie wykresu ugięcia belki dla punktów

plot(coord, u(1:2:Ndofs) ,'r*'),

axis ij, hold on, grid on

for iel=1:Nel pętla rysująca wykres aproksymacji - linia łącząca punkty

n1=iel;

x1=coord(n1);

n2=n1+1;

x2=coord(n2);

ndof_el=Edof(iel,2:5);

uel=u(ndof_el);

w_hel=subs(N,h,(x2-x1))*uel;

w_hel_m=matlabFunction(w_hel);

xx=linspace(x1,x2);

w=w_hel_m(xx-x1);

plot(xx,w,'b')

end

for iel=1:Nel obliczenie macierzy sił przywęzłowych f_el
     ndof_el=Edof(iel,2:5)
     uel=u(ndof_el) ndof_el - globalna liczba stopni swobody elementów
     fel(1:4,iel)=Kel*uel-Q(iel)*qel
 end

figure(2) obliczanie oraz rysowanie wykresów momentów zginających i siły tnących

subplot(2,1,1), title('bending moment [kNm]')

hold on

plot([0,L],[0,0],'k','LineWidth',2)

axis ij

subplot(2,1,2), title('shear force [kN]')

hold on

plot([0,L],[0,0],'k','LineWidth',2)

for iel=1:Nel

n1=iel;

x1=coord(n1);

n2=n1+1;

x2=coord(n2);

ndof_el=Edof(iel,2:5);

uel=u(ndof_el);

M_el=subs(B,h,(x2-x1))*uel;

M=matlabFunction(M_el);

Q_el=diff(M_el);

Q=@(x)double(Q_el);

xx=linspace(x1,x2);

subplot(2,1,1),plot(xx,M(xx-x1))

subplot(2,1,2),plot(xx,Q(xx-x1))

end

Wyniki:

Sprawdzenie wyników:

wnioski: W moim sprawozdaniu uwzględniłem dwa wykresy dla różniej ilości elementów. Pierwszy wykres obrazuje podział przy 4 elementach, natomiast drugi przy 100. Widać, że przy zwiększeniu ilości elementów nasz wykres lepiej odwzorowuje rzeczywiste wartości.