Date 14.5.99
Student 45834H, Mikko Hämäläinen

Final Project

Link Structures

Here is a way to solve link structures with Matlab. To obtain some more information about the theory behind this Matlab code contact Heikki Apiola, the lecturer of the course.
File: ristikko.m



% Ratkaistaan ristikkorakenteen siirtymät elementtimenetelmällä
% e viittaa aina elementtikohtaisiin matriiseihin

clear;

% Parametrien maarittelyja
d=0.01;         % Sauvan halkaisija
A=pi*(d/2)^2;   % Poikkipinnan ala
E=200e9;        % Kimmokerroin
EA=E*A;         % Aksiaalijäykkyys

% Solmupisteet
Ns=9;
X=[0 0 1 1 2 2 3 3 4    % x-koordinaatit
   0 1 0 1 0 1 0 1 0];  % y-koordinaatit

% Reunaehdot ( 1=kiinnitetty, 0=vapaa )
D=[1 1 0 0 0 0 0 0 0    % x-siirtymä
   1 1 0 0 0 0 0 0 0];  % y-siirtymä

% Ulkoiset pistekuormat
L=[0 0 0 0 0 0 0 0 0            % x-voima
   0 0 0 0 0 0 0 0 -1e5];       % y-voima

% Elementit
Ne=15;
E=[1 1 2 2 3 4 4 5 6 6 7 8 3 5 7        % Solmusta Xk
   2 3 4 3 5 6 5 7 8 7 9 9 4 6 8];      % solmuun Xl

% Piirretään rakenne kuvaan 1
figure(1); clf; hold on;
for e=1:Ne;
        k=E(1,e); l=E(2,e);     % Elementin e solmujen numerot k ja l
        Xk=X(:,k); Xl=X(:,l);   % Solmupisteet Xk ja Xl
        plot([Xk(1) Xl(1)],[Xk(2) Xl(2)]);
end;
axis image; axis off;

disp('Paina namiskuukkelia');
pause

% Muodostetaan muunnosmatriisi T
Nv=0;
for i=1:Ns;
        if D(1,i) == 0; Nv=Nv+1; T(2*i-1,Nv)=1; end;
        if D(2,i) == 0; Nv=Nv+1; T(2*i  ,Nv)=1; end;
end;

% Muodostetaan vapaan rakenteen jäykkyysmatriisi elementeittäin
K=zeros(2*Ns,2*Ns);
for e=1:Ne;
        k=E(1,e); l=E(2,e);     % Elementin solmupisteiden numerot
        Xk=X(:,k); Xl=X(:,l);   % Elementin solmupisteiden koordinaatit
        he=norm(Xk-Xl);         % Elementin pituus
        de=(Xk-Xl)/he;          % Elementin yksikkösuuntavektori
        De=[de' [0 0]           % Suuntamatriisi
            [0 0] de'];
        ke=(EA/he)*[1 -1;-1 1]; % Aksiaalijäykkyysmatriisi
        Ke=De'*ke*De;           % Elementin jäykkyysmatriisi
        Ae=zeros(4,2*Ns);       % Asennusmatriisi
        Ae(1,2*k-1)=1; Ae(2,2*k)=1; Ae(3,2*l-1)=1; Ae(4,2*l)=1;
        K=K+Ae'*Ke*Ae;          % Summataan globaaliin jäykkyysmatriisiin
end;

% Muodostetaan vapaan rakenteen voimavektori (vain solmupistevoimia)
F=zeros(2*Ns,1);
for i=1:Ns; F(2*i-1)=L(1,i); F(2*i)=+L(2,i); end;

% Ratkaistaan vapausasteet (red=redusoitu matriisi)
Kred=T'*K*T; Fred=T'*F; Ured=inv(Kred)*Fred; U=T*Ured;

% Tulostetaan siirtymät
for e=1:Ne;
        k=E(1,e); l=E(2,e);     % Elementin solmupisteiden numerot
        Xk=X(:,k); Xl=X(:,l);   % Elementin solmupisteiden koordinaatit
        Ae=zeros(4,2*Ns);       % Asennusmatriisi
        Ae(1,2*k-1)=1; Ae(2,2*k)=1; Ae(3,2*l-1)=1; Ae(4,2*l)=1;
        Ue=Ae*U;
        plot([Xk(1)+Ue(1) Xl(1)+Ue(3)],[Xk(2)+Ue(2) Xl(2)+Ue(4)], '--');
end;

axis image; axis off;





I'm very sorry that the comments are in Finnish only... :(



Check also the PDE ToolBox work, which is much like the problem described, but NOT yet the same problem... ...the plate has a different kind of mesh and it doesn't behave like the link structure. That's partly because of the volume force instead of a nodal force.

The file is: pdeuy.m