% krata - rozwiazanie MES

% przy uzyciu pakietu CALFEM 

%

% opracowal: Adam Wosatko

clear;



% dane

% sztywnosc

E=1e+5; A=1e-2; ep=[E A];

% wymiary

a=2; b=2; c=2;

% obciazenie

P=10; dx=0.05;



% dyskretyzacja

LE=7

LW=5

LSSW=2;

LSSU=LW*LSSW

K=zeros(LSSU);

F=zeros(LSSU,1);



% wspolrzedne

Coord=[0   0;

       a   0;

       a+b 0;

       0   c;

       a   c];

% topologia (przynaleznosc stopni swobody do ES)

% numery st.sw. dla wezlow

Dof=[ 1 2;

      3 4;

      5 6;

      7 8;

      9 10];

% numery st.sw. dla poszczegolnych ES  

Edof=[1  1 2   7 8;

      2  1 2   3 4;

      3  7 8   3 4;

      4  7 8   9 10;

      5  3 4   9 10; 

      6  3 4   5 6;

      7  9 10  5 6];

[Ex,Ey]=coordxtr(Edof, Coord, Dof, 2)

pause;



% rysunek

% siatka ES

Elnum=Edof(:,1);

eldraw2(Ex,Ey,[1,3,1],Elnum);



% si³a przywezlowa

F(10)=P;



% macierze ES, agregacja

for i=1:LE

    Ke=bar2e(Ex(i,:),Ey(i,:),ep);

    K=assem(Edof(i,:),K,Ke);

end



% warunki brzegowe + wymuszenie kinematyczne

bc=[1 0; 2 0; 6 -dx; 7 0];



% rozwiazanie

[a,R]=solveq(K,F,bc)

pause;



% powrot do ES

% przemieszczenia w poszczegolnych ES

Ed=extract(Edof,a)



% rysunek

% przemieszczenia

[sfac]=scalfact2(Ex,Ey,Ed,0.2);

figure(1)

eldisp2(Ex,Ey,Ed,[1,2,1],sfac);

pltscalb2(sfac,[5e-2 3 2]);

title('deformacja')



% si³y w pretach

for i=1:LE

    N(i)=bar2s(Ex(i,:),Ey(i,:),ep,Ed(i,:));

end

N_ES = [N; N]



% rysunek

% sily w pretach

hold off;

figure(2)

plotpar=[2 1];

sfac=scalfact2(Ex(7,:),Ey(7,:),N(7),0.2);

for i=1:LE

eldia2(Ex(i,:),Ey(i,:),N_ES(:,i),plotpar,sfac);

end

pltscalb2(sfac,[10 3 2]);

title('sily w pretach')