%_______________________________________________________________
%2Node-6DOF SPACE TRUSS SYSTEM FEM [A.] (C)(R) |
%############################# |
% Shape function : Lineer (Axial normal) |
% Stifness : Topology +Accumulate method |
% Solve Equation : Cholesky-[L][D][u] |
%_______________________________________________________________|
% Maxima 5.9.0.9beta2 http://maxima.sourceforge.net |
% Using Lisp Kyoto Common Lisp GCL 2.6.3 (aka GCL) |
% Distributed under the GNU Public License. See the file COPYING|
% Dedicated to the memory of William Schelter. |
% This is a development version of Maxima. The function bug_rep |
% provides bug reporting information. |
%_______________________________________________________________|
clc
clear
%## INPUT DATA ##
%% ## MATERIAL PROPERTIES ##
%Materials matrix
% Cross-sectional areas(m**2) Elasticity (KN/m**2)
%_______________________________________Example (1)
%Properties=[0.050 2E7 %1 Element No
% 0.050 2E7 %2
% 0.050 2E7 %3
% 0.050 2E7 ];%4
%___________________________|
%_______________________________________Example (2)
%Properties=[0.020 2E7 %1 Element No
% 0.020 2E7 %2
% 0.020 2E7 %3
% 0.020 2E7 %4
% 0.020 2E7 ];%5
%___________________________|
%_______________________________________Example (3)
%Properties=[0.060 2E7 %1 Element No
% 0.060 2E7 %2
% 0.060 2E7 %3
% 0.060 2E7 %4
% 0.060 2E7 %5
% 0.060 2E7 %6
% 0.060 2E7 %7
% 0.060 2E7 ];%8
%___________________________|
%_______________________________________Example (4)
%Properties=[0.050 2E7 %1 Element No
% 0.050 2E7 %2
% 0.050 2E7 %3
% 0.050 2E7 %4
% 0.050 2E7 ] %5
%___________________________|
%_______________________________________Example (5)
Properties=[0.050 2E7 %1 Element No
0.050 2E7 %2
0.050 2E7 %3
0.050 2E7 %4
0.050 2E7 %5
0.050 2E7 %6
0.050 2E7 ]%7
%for i=1:7
% Properties(i,:)=[0.050 2E7];
%end
%___________________________|
No=size(Properties,1); %Elements Number
%All element's Areas A(i) and Elasticity E(i)
for i=1:No
A(i,:)=Properties(i,1);
E(i,:)=Properties(i,2);
end
%% ## GLOBAL NODE COORDINATES ##
%_______________________________________Example (1)
% Cor=[ X Y Z ]
% Cor=[ 0.00 3.00 4.00
% 0.00 3.00 -4.00
% 0.00 6.00 0.00
% 0.00 0.00 0.00
% 4.00 3.00 0.00 ];
%_______________________________|
%_______________________________________Example (2)
%Cor=[ X Y Z ]
% Cor=[ 0.00 4.00 4.00
% 8.00 4.00 0.00
% 4.00 4.00 -4.00
% 4.00 0.00 4.00 ];
%_______________________________|
%_______________________________________Example (3)
% Cor=[ X Y Z ]
% Cor=[ 0.00 12.00 0.00
% 3.00 3.00 0.00
% 6.00 12.00 0.00
% 3.00 0.00 4.00
% 3.00 0.00 -4.00 ];
%_______________________________|
%_______________________________________Example (4)
% Cor=[ X Y Z ]
% Cor=[ 0.00 6.00 4.00
% 16.00 6.00 4.00
% 16.00 6.00 0.00
% 0.00 6.00 0.00
% 8.00 3.00 0.00
% 8.00 0.00 4.00 ];
%_______________________________|
%_______________________________________Example (5)
% Cor=[ X Y Z ]
Cor=[ 3.00 0.00 4.00
5.00 0.00 4.00
4.00 0.00 -4.00
0.00 4.00 0.00
4.00 8.00 0.00
8.00 4.00 0.00 ];
%_______________________________|
Node=size(Cor,1);
%% ## ELEMENT POSITION MATRIX ##
%_______________________________________Example (1)
%Pos=[3 5
% 2 5
% 1 5
% 4 5];
%________________________________|
%_______________________________________Example (2)
%Pos=[1 3
% 1 2
% 4 3
% 3 2
% 4 2];
%________________________________|
%_______________________________________Example (3)
%Pos=[4 1
% 1 5
% 4 2
% 2 3
% 4 3
% 5 2
% 5 3
% 2 1];
%________________________________|
%_______________________________________Example (4)
%Pos=[4 5
% 1 6
% 5 3
% 6 2
% 6 5];
%________________________________|
%_______________________________________Example (5)
Pos=[4 5
1 4
1 5
6 5
2 5
2 6
3 5];
%________________________________|
%Element`s Node Topology
for i=1:Node;
Re(i,:)=[1 1 1];
end
%% ## SYSTEM SUPPORTS ##
%Re(Node number,:)=[ux vy wz](0=Connected,1=free)
%_______________________________________Example (1)
% Re(1,:)= [0 0 0];
% Re(2,:)= [0 0 0];
% Re(3,:)= [0 0 0];
% Re(4,:)= [0 0 0];
%_______________________|
%_______________________________________Example (2)
% Re(1,:)= [0 0 0];
% Re(2,:)= [0 0 1];
% Re(4,:)= [0 0 0];
%_______________________|
%_______________________________________Example (3)
% Re(1,:)= [1 0 0];
% Re(3,:)= [1 0 0];
% Re(4,:)= [0 0 0];
% Re(5,:)= [0 0 0];
%_______________________|
%_______________________________________Example (4)
% Re(1,:)= [0 0 0];
% Re(2,:)= [0 0 0];
% Re(3,:)= [0 0 0];
% Re(4,:)= [0 0 0];
% Re(6,:)= [1 0 0];
%_______________________|
%_______________________________________Example (5)
Re(1,:)= [0 0 0];
Re(2,:)= [0 0 0];
Re(3,:)= [0 0 0];
Re(4,:)= [1 0 0];
Re(6,:)= [1 0 0];
%_______________________|
Topology=Re;
Nom=size(Re,2); % Element node d.o.f Number
%Topology---->Accumulate
sayman=0;
for i=1:Node;
for j=1:Nom;
if Re(i,j)==1 ;
sayman = sayman +1;
Re(i,j) = sayman;
end
end
end
Item=sayman; %System total displacement value
for i=1:No
R(i,:)=[Re(Pos(i,1),:) Re(Pos(i,2),:)];
end
P(Item)=0;
%% ## SYSTEM GLOBAL LOAD ##
%P(Re(Node number, freedom number))=LOAD
%_______________________________________Example (1)
%P(Re(5,1))= 10.00;
%P(Re(5,2))= -8.00;
%__________________________|
%_______________________________________Example (2)
%P(Re(3,2))= -10.00;
%__________________________|
%_______________________________________Example (3)
%P(Re(1,1))= -10.00;
%P(Re(2,2))= -30.00;
%P(Re(3,1))= 10.00;
%__________________________|
%_______________________________________Example (4)
% P(Re(5,1))= 10.00;
% P(Re(5,2))= -10.00;
% P(Re(6,1))= -10.00;
%__________________________|
%_______________________________________Example (5)
P(Re(4,1))= 10.00 ;
P(Re(5,2))=-30.00 ;
P(Re(6,1))=-10.00 ;
%__________________________|
% ###########
%SYSTEM_ANALYSIS
%
%% Global system node is moving per element nodes
for s=1:No;
px1(s,:)=Cor(Pos(s,1),1);
px2(s,:)=Cor(Pos(s,2),1);
py1(s,:)=Cor(Pos(s,1),2);
py2(s,:)=Cor(Pos(s,2),2);
pz1(s,:)=Cor(Pos(s,1),3);
pz2(s,:)=Cor(Pos(s,2),3);
end
%% Element Stiffness Applications
for s=1:No
%Element length
Ln = sqrt( (px2(s)-px1(s))^2+(py2(s)-py1(s))^2+(pz2(s)-pz1(s))^2 );
%Element direction cosinesses of global coordinat system
Dcos= [ (px2(s)-px1(s)) %ly
(py2(s)-py1(s)) %my
(pz2(s)-pz1(s))]/Ln; %ny
Dataport(s,:)=[Ln Dcos(1) Dcos(2) Dcos(3)];
%Analytical tranformation matrix
%[T]= [ ly(i) my(i) ny(i) 0 0 0
% -my(i) ly(i) 0 0 0 0
% -ny(i) 0 ly(i) 0 0 0
% 0 0 0 ly(i) my(i) ny(i)
% 0 0 0 -my(i) ly(i) 0
% 0 0 0 ny(i) 0 ly(i) ]
%Global-->Local axis transformation matrix
T(:,:,s)= [ Dcos(1) Dcos(2) Dcos(3) 0 0 0
-Dcos(2) Dcos(1) 0 0 0 0
-Dcos(3) 0 Dcos(1) 0 0 0
0 0 0 Dcos(1) Dcos(2) Dcos(3)
0 0 0 -Dcos(2) Dcos(1) 0
0 0 0 -Dcos(3) 0 Dcos(1)];
%Per element local axis stiffness matrix
K(:,:,s)=[ A(s)*E(s)/Ln 0.00 0.00 -A(s)*E(s)/Ln 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
-A(s)*E(s)/Ln 0.00 0.00 A(s)*E(s)/Ln 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 ];
%Element global axis stiffness matrix
%[K]global=[C].[K]local.[C]' or
%[K]global=[T]'[K]local.[T]
Kg(:,:,s)=T(:,:,s)'*K(:,:,s)*T(:,:,s);
end
%% System Stiffness Matrix (Superposition)
%Element stiffness matrix superposition of system global axis
Ksis(Item,Item)=0;
for n=1:No;
for sat=1:6;
for sut=1:6;
if (R(n,sat)~=0)
if (R(n,sut)~=0);
Ksis(R(n,sut),R(n,sat))=Ksis(R(n,sut),R(n,sat)) + Kg(sat,sut,n);
end
end
end
end
end
%% System global stiffness matrix singularity check
equation=size(Ksis);
if equation(1)~=rank(Ksis)
display('This system stiffness matrix is badly scaled')
R
error('Control system support boundary conditions')
else
Ku = inv(Ksis);
D = Ku *P';
end
clear equation
%Element per unit node topology matrix
for v = 1 : No;
for m = 1 :2*Nom;%6 degree of freedom
u = R(v, m);
if u ~=0
Hu(v, m,:) = D(u) ;
else
Hu(v,m)=0;
end
end
end
%% Element global and local node's reactions
for s=1 :No;
%Element global node reactions
Pg(:,s) = Kg(:,:,s)*Hu(s,:)';
%Element local node reactions
Pl(:,s) = T(:,:,s)*Pg(:,s);
end
Pf(3*Node)=0;
Af(3*Node)=0;
for s=1:No; %Total element number
for Nm=1:2; %One node d.o.f
Poss=Pos(s,Nm);
Pf((Poss-1)*3+1:Poss*3)=Pg((Nm-1)*3+1:Nm*3,s);
Af((Poss-1)*3+1:Poss*3)=Af((Poss-1)*3+1:Poss*3)+Pf((Poss-1)*3+1:Poss*3);
end
end
clear Nm
% for j=1:size(Re,1)
% if Topology(j,:)~=[1 1 1]
% Af(3*(j-1)+1:3*(j-i)+3)'
% end
% end
%% Analysis Results
%_______________________________SYSTEM DESCRIPTION
display ('====SYSTEM ANALYSIS PROPERTIES');
fprintf(' \n');
fprintf(' System total degree of freedom= [%.f] \n',Item);
fprintf(' System total element value = [%.f] \n',No);
fprintf(' One node degree of freedom = [%.f] \n',Nom);
%whos
fprintf('\n \n ');
%_______________________________PROPERTIES
display ('====PER ELEMENTS PROPERTIES');
fprintf(' \n');
fprintf(' Area(m**2) E. Module (KN/m**2) \n');
for s=1:size(Properties,1);
fprintf(' (%.f) [% .4f] [% .f] \n',s,Properties(s,:));
end
fprintf('\n \n ');
%_______________________________COORDINATES
display ('====ALL NODE GLOBAL CARTESIAN COORDINATE (metre)');
fprintf(' \n');
fprintf(' Xglobal Yglobal Zglobal \n');
for s=1:size(Cor,1);
fprintf(' (%.f) [% .4f] [% .4f] [% .4f] \n',s,Cor(s,:));
end
fprintf('\n \n ');
%_______________________________ELEMENT Length and Direction Coninesses
display ('====ELEMENT SPACE POSITON (metre)');
fprintf(' \n');
fprintf(' Length(m.) ly my ny (Direction Cosinesses) \n');
for s=1:size(Dataport,1);
fprintf(' (%.f) [% .2f] [% .4f] [% .4f] [% .4f] \n',s,Dataport(s,:));
end
fprintf('\n \n ');
%_______________________________POSITION MATRIX
display ('====ALL ELEMENT (i) and (j) NODE POSITION ');
fprintf(' \n');
fprintf(' Pos(i) Pos(j) \n');
for s=1:size(Pos,1);
fprintf(' (%.f) (% .f)o--->---o(% .f) \n',s,Pos(s,:));
end
fprintf('\n \n ');
%_______________________________TOPOLOGY
display ('====SYSTEM ALL NODE TOPOLOGY MATRIX ');
fprintf(' \n');
fprintf(' u(x) v(y) w(z) (Connected=0 Free=1) \n');
for s=1:size(Topology,1);
fprintf(' (%.f) [% .f] [% .f] [% .f] \n',s,Topology(s,:));
end
fprintf('\n \n ');
%_______________________________REOLOGY
display ('====SYSTEM ALL NODE REOLOGY MATRIX ');
fprintf(' \n');
fprintf(' u(x) v(y) w(z) (Connected=0 else Free) \n');
for s=1:size(Re,1);
fprintf(' (%.f) [% .f] [% .f] [% .f] \n',s,Re(s,:));
end
fprintf('\n \n ');
%_______________________________LOCAL [K]local
display ('====PER ELEMENT LOCAL STIFFNESS MATRIX ');
fprintf(' \n');
K
fprintf('\n \n ');
%_______________________________TRANSFORM [T]
display ('====PER ELEMENT GLOBAL to LOCAL TRANSFORM MATRIX ');
fprintf(' \n');
T
fprintf('\n \n ');
%_______________________________GLOBAL [K]global
display ('====PER ELEMENT GLOBAL STIFFNESS MATRIX ');
fprintf(' \n');
Kg
fprintf('\n \n ');
%_______________________________SYSTEM STIFFNESS [K]system
display ('====SYSTEM STIFFNESS MATRIX ');
fprintf(' \n');
Ksis
fprintf('\n \n ');
%_______________________________SYSTEM FORCES [P]
display ('====SYSTEM STIFFNESS MATRIX ');
fprintf(' \n');
Ksis
fprintf('\n \n ');
%_______________________________GLOBAL NODE DISPLACEMENT
display ('====PER ELEMENTS NODE GLOBAL DISPLACEMENT===[Du Dv Dw]global (metre)');
fprintf('\n');
fprintf(' u(i) v(i) w(i) u(j) v(j) w(j) \n');
for s=1:size(Hu,1);
fprintf(' (%.f) [% .7f] [% .7f] [% .7f] [% .7f] [% .7f] [% .7f] \n',s,Hu(s,:));
%fprintf(' [% .7f]%.f [% .7f]%.f [% .7f]%.f [% .7f]%.f \n',Pg(s,1),R(1,s)',Pg(s,2),R(2,s)',Pg(s,3),R(3,s)',Pg(s,4),R(4,s)');
%if mod(s,3)==0; fprintf(' \n');end;
end
fprintf('\n \n ');
%_______________________________GLOBAL NODE DISPLACEMENT
display ('====PER ELEMENTS NODE GLOBAL DISPLACEMENT===[Du Dv Dw]global (metre)');
fprintf('\n');
fprintf(' u(i) v(i) w(i) u(j) v(j) w(j) \n');
for s=1:size(Hu,1);
fprintf(' (%.f) [% .7f] [% .7f] [% .7f] [% .7f] [% .7f] [% .7f] \n',s,Hu(s,:));
%fprintf(' [% .7f]%.f [% .7f]%.f [% .7f]%.f [% .7f]%.f \n',Pg(s,1),R(1,s)',Pg(s,2),R(2,s)',Pg(s,3),R(3,s)',Pg(s,4),R(4,s)');
%if mod(s,3)==0; fprintf(' \n');end;
end
fprintf('\n \n ');
%_______________________________GLOBAL NODE REACTIONS
display ('====PER ELEMENTS GLOBAL NODE REACTION===[Fx Fy Fz]global (KN)');
fprintf('\n');
fprintf(' [fx(i) fy(i) fz(i)]o-------o[fx(j) fy(j) fz(j)] \n');
for s=1:size(Pg,2);
fprintf(' (%.f) [% .5f] [% .5f] [% .5f] [% .5f] [% .5f] [% .5f] \n',s,Pg(:,s)');
end
fprintf('\n \n ');
%_______________________________LOCAL NODE REACTIONS
display ('====PER ELEMENTS LOCAL NODE REACTION===[Fx Fy Fz]global (KN)');
fprintf('\n');
fprintf(' [ fnormal(i) ]o--->----o[ fnormal(j) ] \n');
for s=1:size(Pl,2);
if Pl(1,s) <= 0
fprintf(' (%.f) [% .5f] [% .5f] (Tansion effect) stretch motion \n',s,Pl(1,s)',Pl(4,s)');
else
fprintf(' (%.f) [% .5f] [% .5f] (Pressure effect) buckling motion\n',s,Pl(1,s)',Pl(4,s)');
end
end
fprintf('\n \n ');
%_______________________________NODE and SUPPORT REACTIONS (superposition)
display ('====SYSTEM ALL NODE and SUPPORT REACTIONS===[Fx Fy Fz]global (KN)');
fprintf('\n');
fprintf(' fx(g) fy(g) fz(g) \n');
for s=1:size(Af)';
fprintf(' [% .5f] [% .5f] [% .5f] \n',Af');
end
fprintf('\n \n ');
