belka.m

% rozwiazanie MES
% przy uzyciu pakietu CALFEM 
% (belka - ale ES ramowe!)
%
% opracowal: Piotr Pluciński
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% dane
E=250e4;
I=1e-4;
A=1e-2;
ep = [E A I]
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% dyskretyzacja
LE=4; LW=LE+1; LSSW=3; 
LSSU=LW*LSSW;
K=zeros(LSSU);
F=zeros(LSSU,1);
 
%wspolrzedne
Coord=[0 0; 1 0; 2 0; 6 0; 8 0];
 
% topologia 
% (przynaleznosc stopni swobody do ES)
% numery st.sw. dla wezlow
Dof=reshape((1:LSSU),LSSW,length(Coord))';
 
% macierz topologii
for i=1:LE
   Etop(i,:)=[i, i+1];
end
 
% definicja macierzy stopni swobody dla elementów
for i=1:LE
   Edof(i,:)=[i,Dof(Etop(i,1),:),Dof(Etop(i,2),:)];
end
 
[Ex,Ey]=coordxtr(Edof, Coord, Dof, 2)
 
% rysunek siatki ES
elnum=Edof(:,1);
eldraw2(Ex,Ey,[1,3,1],elnum);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% obciazenie:
% siły węzłowe
F(5)=-20;
F(12)=5;
 
% równomierne obcišżenie
eq=zeros(LE,2);
eq(3,:)=[0,-10];
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% macierze i wektory ES
% agregacja
 
for i=1:LE
    [Ke, Fe]=beam2e(Ex(i,:),Ey(i,:),ep,eq(i,:));
    [K,F]=assem(Edof(i,:),K,Ke,F,Fe);
end %for 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% warunki brzegowe
bc=[1 0; 2 0; 3 0; 8 0; 11 0; 13 0; 15 0];
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%rozwiazanie
[Q,R]=solveq(K,F,bc)
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% postprocessing
% powrot do ES
% przemieszczenia w poszczegolnych ES
Qe=extract_ed(Edof,Q);
% rysunek zdeformowanej siatki ES
figure(1)
[sfac]=scalfact2(Ex,Ey,Qe,0.2);
eldisp2(Ex,Ey,Qe,[2,1,1],sfac);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% obliczenie sił przywezłowych
for i=1:LE
  fe(:,:,i)=beam2s(Ex(i,:),Ey(i,:),ep,Qe(i,:),eq(i,:),20);
end
fe
plotpar=[2 1];
 
% Siły poprzeczne
figure(2)
[fmax,nmax]=max(max(fe(:,2,:)));
scal=scalfact2(Ex(nmax,:),Ey(nmax,:),fe(:,2,nmax),0.2);
for i=1:LE
eldia2(Ex(i,:),Ey(i,:),fe(:,2,i),plotpar,scal);
end
axis([-1 9 -1 1]);
title('sily poprzeczne')
 
% Momenty
figure(3)
[fmax,nmax]=max(max(fe(:,3,:)));
scal=scalfact2(Ex(nmax,:),Ey(nmax,:),fe(:,3,nmax),0.2);
for i=1:LE
eldia2(Ex(i,:),Ey(i,:),fe(:,3,i),plotpar,scal);
end
axis([-1 9 -1 1]);
title('momenty');