try%!!!!!!!After clear this line
clc%!!!!!!!!!
disp('PLEASE GLADE MATHLAB LINE 334-343')%!!!!!!!!!!
catch%!!!!!!!!!!!!!!
% _________________________________________________________________________
%|____10NODE_20DOF PARAMETRIC TRIANGULAR PLANE STRESS FEM_______[A.]_(C)(R)|
%| |
%| Shape function:Homogen |
%| 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_report() |
%| provides bug reporting information. |
%| |
%|_____PARAMETRIC FEM ANALYSIS SUBPROGRMAM_______Maxima5.9.0.9beta2_[GPL]__|
%Information:
% Pascal Triangle
%
% 1
% x y
% x^2 x*y y^2
% x^3 x^2*y x*y^2 y^3
%
% This finite element model shape function is very complex. Because system
% coefficient matrix dimension is size 6*6 and invers matrix not
% satisfactorly element stiffness matrix for maxima code. High degree
% triangular finite element model`s generally application Enegry II. therom.
% Nevertless shape function`s description homogen axis depend the element
% cartesian local coordinates. This method is more simplest than other
% methods.
%
% Element Homogen Shape Functions
%
% --------------(3)---------------e3=1 (homogen axis)
% / \
% ---------- (8)---(7)------------e3=2/3
% / \
% --------(9)---(10)--(6)---------e3=1/3
% / \
% -------(1)--(4)--(5)--(2)-------e3=0
%
% N3=c*L[1-2]*L[9-6]*L[8-7];
% N3=c*(e3)*(e3-1/3)*(e3-2/3);
% if N3=1 for e3=1 than
% 1=c*1/3*2/3=> c=9/2 and N3=9/2*(e3)*(e3-1/3)*(e3-2/3) other shape function is cyclic permutation
%
% N4=c*L[1-3]*L[5-8]*L[2-3];
% N4=c*(e2)*(e1-1/2)*(e1);
% if N4=1 for e2=2/3,e1=1/3 than
% 1=c*1/3*1/3*2/2 => c=27/2 and N4=27/2*e1*e2*(e1-1/3) other shape function is cyclic permutation
% N10=c*L[1-2]*L[2-3]*L[1-3]
% N10=c*e1*e2*e3
% if N3=1 for e3=1/3,e2=1/3,e3=1/3 than
%
% 1=c*1/3*1/3*1/3* => c=27 and N10=27*e1*e2*e3
%
% _________________________________________________________________________
%|____10NODE_20DOF PARAMETRIC TRIANGULAR PLANE STRESS FEM______[A.]_(C)(R)|
%| |
%| Shape function:Homogen |
%| Stifness :Topology+Accumulate method |
%|_______________________________Solve Equation:Cholesky-[L][D][u]_________|
%| |
%|_____PARAMETRIC FEM ANALYSIS SUBPROGRMAM_______Maxima5.9.0.9beta2_[GPL]__|
clc
clear
%===============================================INPUT DATA
%=============MATERIAL PROPERTIES
E=2E6; %Element elasticity constant;
v=0.3; %Element material poission ratio;
th=0.10; %Element thickness
%=============
%==============ELEMENT POSITION MATRIX
%Pos(Element No,:)=[Element Node Number]
Pos(1,:)= [1 34 31 12 23 33 32 21 11 22];
Pos(2,:)= [1 4 34 2 3 14 24 23 12 13];
Pos(3,:)= [4 37 34 15 26 36 35 24 14 25];
Pos(4,:)= [4 7 37 5 6 17 27 26 15 16];
Pos(5,:)= [7 40 37 18 29 39 38 27 17 28];
Pos(6,:)= [7 10 40 8 9 20 30 29 18 19];
Pos(7,:)= [31 34 51 32 33 43 52 51 41 42];
Pos(8,:)= [61 34 64 52 43 44 54 63 62 53];
Pos(9,:)= [34 37 64 35 36 46 55 54 44 45];
Pos(10,:)=[64 37 67 55 46 47 57 66 65 56];
Pos(11,:)=[37 40 67 38 39 49 58 57 47 48];
Pos(12,:)=[67 40 70 58 48 50 60 69 68 59];
%==============
Number=size(Pos);
No=Number(1);
%_________________Automatic coordinate function
%Cor=Elemaent Node Castesian Coordinate for Global System Axis
% Cor=[ xi yi ]
nomin=0;
for sut=1:7;
for sat=1:10;
nomin=nomin+1;
Cor(nomin,1)=0.25*(sat-1);
Cor(nomin,2)=0.20*(sut-1);
end
end
clear sut sat nomin;
Number=size(Cor);
Node=Number(1); %System Node Number
for i=1:Node;
Re(i,:)=[1 1];
end
%===========SUPPORT REACTION
%Re(Node number,:)=[ux vy]
Re(1,:) =[0 0];
Re(10,:) =[1 0];
Re(61,:) =[0 1];
%============
Number=size(Re);
Nom=Number(2); %Plane Element node d.o.f Number
%______________________Topology and Accumulate method
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; %Plane system sum deplacement value
%_________________________________Modal deplacement parameter
for i=1:No;
for j=1:Nom;
R(i,j+0*Nom)= Re(Pos(i,1),j);
R(i,j+1*Nom)= Re(Pos(i,2),j);
R(i,j+2*Nom)= Re(Pos(i,3),j);
R(i,j+3*Nom)= Re(Pos(i,4),j);
R(i,j+4*Nom)= Re(Pos(i,5),j);
R(i,j+5*Nom)= Re(Pos(i,6),j);
R(i,j+6*Nom)= Re(Pos(i,7),j);
R(i,j+7*Nom)= Re(Pos(i,8),j);
R(i,j+8*Nom)= Re(Pos(i,9),j);
R(i,j+9*Nom)= Re(Pos(i,10),j);
end
end
%================SYSTEM LOAD
%Re(Node number, Freedoom number)[u v]
P(Item)=0;
P(Re(40,1)) =-5;
P(Re(64,2))=-10;
P(Re(67,2))=-10;
%================
%___________________________________________________END_INPUT_DATA
%===================================================SYSTEM_ANALYSE
KSI=E*th/(1-v^2);%Elasticity constant
%Element local axis stiffness matrix
for s=1:No;%Element number
X1=Cor(Pos(s,1),1); Y1=Cor(Pos(s,1),2);
X2=Cor(Pos(s,2),1); Y2=Cor(Pos(s,2),2);
X3=Cor(Pos(s,3),1); Y3=Cor(Pos(s,3),2);
Area=0.5*(-X1*(Y3-Y2)+X2*Y3-X3*Y2+(X3-X2)*Y1);
Ai=( X2*Y3-X3*Y2);
Aj=(-X1*Y3+X3*Y1);
Ak=( X1*Y2-X2*Y1);
Bi=-(Y3-Y2);
Bj=-(Y1-Y3);
Bk=-(Y2-Y1);
Ci=+(X3-X2);
Cj=+(X1-X3);
Ck=+(X2-X1);
%=============================================[PROCEDURE(2)OUTPUT]
%K(1,1,s) = (-17.0)*(Ci^2*v-Ci^2-2*Bi^2)/(160.0*Area);
%K(1,2,s) = 17.0*Bi*Ci*(v+1)/(160.0*Area);
%K(1,3,s) = (-7.0)*(Ci*Cj*v-Ci*Cj-2*Bi*Bj)/(320.0*Area);
%K(1,4,s) = 7.0*(2*Bi*Cj*v-Bj*Ci*v+Bj*Ci)/(320.0*Area);
%K(1,5,s) = (-7.0)*(Ci*Ck*v-Ci*Ck-2*Bi*Bk)/(320.0*Area);
%K(1,6,s) = 7.0*(2*Bi*Ck*v-Bk*Ci*v+Bk*Ci)/(320.0*Area);
%__________________________________________________________________________
% |THIS REGION RUN MAXIMA(5.0.9)& AND ELEMENT STIFFNESS |
%__________|_____________________________________________________|_________
% | Element total Node = 12 |
% | Element total degree of freedom =12*[ux,vy]=24 |
% | Stiffness matrix =K(24,24) |
% |=====================================================|
% |PLEASE RUNING MAXIMA FILE:[maximafulltrnparplane3.m] |
% |=====================================================|
%__________|_____________________________________________________|_________
%K(20,14,s) = (-81.0)*(Bj*Bk*v+2*Bi*Bk*v+Bj^2*v+Bi*Bj*v-2*Cj*Ck-4*Ci*Ck-2*Cj^2-2*Ci*Cj-Bj*Bk-2*Bi*Bk-Bj^2-Bi*Bj)/(320.0*Area);
%K(20,15,s) = (-81.0)*(2*Bj*Ck*v+Bi*Ck*v-4*Bk*Cj*v-2*Bi*Cj*v-2*Bk*Ci*v+Bj*Ci*v-Bi*Ci*v-2*Bj*Ck-Bi*Ck-Bj*Ci-Bi*Ci)/(320.0*Area);
%K(20,16,s) = (-81.0)*(2*Bj*Bk*v+Bi*Bk*v+Bi*Bj*v+Bi^2*v-4*Cj*Ck-2*Ci*Ck-2*Ci*Cj-2*Ci^2-2*Bj*Bk-Bi*Bk-Bi*Bj-Bi^2)/(320.0*Area);
%K(20,17,s) = 81.0*(Bk*Ck*v-Bj*Ck*v+2*Bi*Ck*v+2*Bk*Cj*v+4*Bi*Cj*v-Bk*Ci*v-2*Bj*Ci*v+Bk*Ck+Bj*Ck+Bk*Ci+2*Bj*Ci)/(320.0*Area);
%K(20,18,s) = (-81.0)*(Bk^2*v+Bj*Bk*v+Bi*Bk*v+2*Bi*Bj*v-2*Ck^2-2*Cj*Ck-2*Ci*Ck-4*Ci*Cj-Bk^2-Bj*Bk-Bi*Bk-2*Bi*Bj)/(320.0*Area);
%K(20,19,s) = 81.0*(2*Bk*Ck+Bj*Ck+Bi*Ck+Bk*Cj+2*Bj*Cj+Bi*Cj+Bk*Ci+Bj*Ci+2*Bi*Ci)*(v+1)/(160.0*Area);
%K(20,20,s) = (-81.0)*(Bk^2*v+Bj*Bk*v+Bi*Bk*v+Bj^2*v+Bi*Bj*v+Bi^2*v-2*Ck^2-2*Cj*Ck-2*Ci*Ck-2*Cj^2-2*Ci*Cj-2*Ci^2-Bk^2-Bj*Bk-Bi*Bk-Bj^2-Bi*Bj-Bi^2)/(80.0*Area);
K(:,:,s)=KSI*K(:,:,s);
end
%open dimension system stiffness matrix
Ksis(Item,Item)=0;
for n=1:No;
for sat=1:20;
for sut=1:20;
if (R(n,sat)~=0)
if (R(n,sut)~=0);
Ksis(R(n,sut),R(n,sat))=Ksis(R(n,sut),R(n,sat)) + K(sat,sut,n);
end
end
end
end
end
%System global stiffness matris 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 topolog matrix
for v = 1 : No;
for m = 1 :20;
u = R(v, m);
if u ~=0
Hu(v, m) = D(u) ;
else
Hu(v,m)=0;
end
end
end
%Global node displacement is moving element local node
for s=1:No;
for i=1:10*Nom;
if R(s,i)~=0 ;
Dep(s,i)=D(R(s,i));
else
Dep(s,i)=0;
end
end
end
display('Global displacement [u(x) v(y)]T');
Dep'
for v=1 :No;
Pg(:,v) = K(:,:,v)*Hu(v,:)';
end
display('Global system node reactions');
Pg
contourf(Ksis)
title('Global system stiffness band matrix');
%==========================================Manual Patch Test
% [1] [2] [3] [4] [5] [6]
%[ 3.2056] [ 3.6147] [ 0.1102] [-0.9130] [ 0.1634] [-1.2703]
%[ 9.2480]1 [ 2.3245]1 [ 0.1683]4 [ 0.0281]4 [ 0.0467]7 [-0.2409]7
%[-0.5523] [ 0.8029] [-0.1623] [ 1.1068] [-1.0886] [ 0.0000]
%[-0.6890]34[-0.1964]4 [ 0.0178]7 [ 0.1942]7 [-0.1183]40 [ 8.4275]10
%[-0.4381] [-0.9990] [-0.2137] [ 0.2313] [ 0.8412] [-0.5903]
%[-1.2221]31[-0.4747]34 [-0.1577]37 [-0.2652]37 [-0.6446]37 [-3.4487]40
%[ 2.7847] [ 0.0000] [ 3.1421] [ 0.0000] [ 2.7068] [-0.0000]
%[ 0.8902]12[-0.0000]2 [-0.4679]5 [-0.0000]5 [-1.7349]18 [ 0.0000]8
%[ 1.6300] [ 0.0000] [-0.2929] [-0.0000] [-4.5674] [-0.0000]
%[-1.7195]23[ 0.0000]3 [ 0.7095]6 [ 0.0000]6 [ 6.4727]29 [ 0.0000]9
%[-3.6309] [ 1.5914] [-0.1652] [ 3.3741] [ 1.6148] [-0.0000]
%[-2.3295]33[-0.3016]14 [-1.2392]17 [-0.3760]17 [-3.5490]39 [ 0.0000]20
%[-2.9989] [-0.5953] [-1.4221] [-0.9501] [ 2.7539] [ 0.0000]
%[-4.1781]32[-2.1811]24 [-1.5136]27 [ 0.6604]27 [-0.1882]38 [-0.0000]30
%[-0.0000] [-1.6300] [ 0.5953] [ 0.2929] [ 0.9501] [ 4.5674]
%[-0.0000]21[ 1.7195]23 [ 2.1811]26 [-0.7095]26 [-0.6604]27 [-6.4727]29
%[ 0.0000] [-2.7847] [-1.5914] [-3.1421] [-3.3741] [-2.7068]
%[ 0.0000]11[-0.8902]12 [ 0.3016]15 [ 0.4679]15 [ 0.3760]17 [ 1.7349]18
%[ 0.0000] [ 0.0000] [ 0.0000] [ 0.0000] [-0.0000] [ 0.0000]
%[ 0.0000]22[-0.0000]13 [ 0.0000]16 [-0.0000]16 [ 0.0000]28 [-0.0000]19
% [7] [8] [9] [10] [11] [12]
%[ 0.4381] [-1.8203] [ 1.2397] [ 2.4183] [ 0.1083] [ 2.0129]
%[ 1.2221]31[ 0.0000]61 [ 1.4287]34 [-0.8946]64 [ 0.3294]37 [-1.6317]67
%[-0.0808] [ 0.6062] [-0.6622] [-0.3564] [-1.0495] [-2.2715]
%[-0.5318]34[ 0.4245]34 [ 0.0844]37 [ 0.4781]37 [ 0.1990]40 [ 3.3681]40
%[-0.7305] [-3.7072] [ 1.2890] [-3.3020] [ 1.2891] [-0.0000]
%[ 0.0145]51[-5.4230]64 [-3.6824]64 [-4.7366]67 [-3.6318]67 [ 0.0000]70
%[ 2.9989] [ 1.2165] [ 1.4221] [ 2.5946] [-2.7539] [-1.8204]
%[ 4.1781]32[ 2.3893]52 [ 1.5136]35 [ 0.8069]55 [ 0.1882]38 [ 3.0499]58
%[ 3.6309] [ 5.7706] [ 0.1652] [ 2.9249] [-1.6148] [ 2.0790]
%[ 2.3295]33[ 4.8087]43 [ 1.2392]36 [ 1.9760]46 [ 3.5490]39 [-4.7864]48
%[-5.7706] [-0.5265] [-2.9249] [-1.9634] [ 0.0000] [-0.0000]
%[-4.8087]43[-1.7933]44 [-1.9760]46 [ 1.1817]47 [-0.0000]49 [-0.0000]50
%[-1.2165] [-1.5393] [-2.5946] [-2.3160] [ 1.8204] [-0.0000]
%[-2.3893]52[-0.4060]54 [-0.8069]55 [ 1.1886]57 [-3.0499]58 [ 0.0000]60
%[ 0.7305] [-0.0000] [ 1.5393] [ 0.0000] [ 2.3160] [ 0.0000]
%[-0.0145]51[-0.0000]63 [ 0.4060]54 [ 0.0000]66 [-1.1886]57 [ 0.0000]69
%[ 0.0000] [ 0.0000] [ 0.5265] [ 0.0000] [ 1.9634] [-0.0000]
%[ 0.0000]41[ 0.0000]62 [ 1.7933]44 [ 0.0000]65 [-1.1817]47 [ 0.0000]68
%[-0.0000] [-0.0000] [ 0.0000] [-0.0000] [-2.0790] [ 0.0000]
%[-0.0000]42[ 0.0000]53 [ 0.0000]45 [-0.0000]56 [ 4.7864]48 [-0.0000]59
%=================================Node Patch Testing Code
for sutun=1:No;
for sat=1:10;
Pgx(sat,sutun)=Pg(2*sat-1,sutun);
Pgy(sat,sutun)=Pg(2*sat,sutun);
end
end
Pdx(Node)=0;
Pdy(Node)=0;
for Elemanno=1:No;
for sutun=1:10;
Pdx(Pos(Elemanno,sutun))=Pdx(Pos(Elemanno,sutun))+Pgx(sutun,Elemanno);
Pdy(Pos(Elemanno,sutun))=Pdy(Pos(Elemanno,sutun))+Pgy(sutun,Elemanno);
end
end
%Information:This Pdx,Pdy different the zero terms Support Reactions
('Manual Patch Test (node kin. eq. met. [Pgx Pgy])')
[Pdx' Pdy']
end%!!!!!!!!!!!!!!!!!!!!!!
