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