%__________________________________________________________________________
%| 4NODE-20DOF CURVED PARAMETRIC SHELL FEM [A.] (C)(R) |
%|######################################## |
%| Element family:Low proclivity curved shell |
%| Shape function:Workable(geometric boundary condition)|
%| Stifness :Topology + Accumulate method |
%| Theoretic app.Total potantial enegry second theorem |
%| Solve Equation:Cholesky-[L][D][u] |
%| |
%|_______________________________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. |
%|_________________________________________________________________________|
% This finite element model`s describtion 4 Node and per node`s local
%displacements [ux,vy,wz,Qx,Qy] respectively. Element shape function emanated
%to two different region. First region; plane stress finite element motion ,
%(membran effect). Second region: Kirchoff bending plate element motion. Per
%region depend to rigidy body motion workable the element`s geometric boundary
%conditions. This model if components 20 >=fo/h >=6 and (h/Rmin)<=1/20 than
%moment and force supply kinematic equilibrium equation else it`s not
%conservation. Nevertheless this model if kx=ky=kxy=0 than transformed bending
%plate moniton.(Example cantilever beam analysis). Element`s local axis per node
%total 5 d.o.f but global node 6. d.o.f.. Element global axis curvature kgx,kgy,
%kgxy tranform local axis curvature kx,ky,kxy. This area kgxy=kxy=0 neglected
%because element is geometry symmetric. If kgxy or kxy different to zero than
%element shape function applied torsional curvature coefficient.
%========================================================This model`s REFERENCES
% [1] "Sonlu Elemanlar Metodu Temelleri ve Yap Mekaniinde Uygulamalar"
% (Fundamental Finite Element Method and Structural Application),Prof.Azer.A.
% KASIMZADE, Birsen yayn evi(press) ,Samsun 19 Mays University, 2004,
% ISBN:975-511-379-7.
%
% [2] Zienkiewich O.C and Taylor R.L,"The finite element Method Formulation and
% Linear Problems",Mc-Graw-Hill LONDON,1989, 4th edition
%
% [3] Zienkiewich O.C and Taylor R.L,"The finite element Method Formulation Solid
% Mechanic",volume II,ISBN: 0-7506-5055-9,Barcelona,Spain,MPG books Ltd, 1999 Aug.
%
%#############THEORICALY APPLICATION###########
%
%____________1.Describtion Element`s shape function
%
%Membrane effect u(x,y)=c1 + c2*x + c3*y + c4*xy
% v(x,y)=c5 + c6*x + c7*y + c8*xy
%
% 1
% Pascal Triangle ________ x y _____Lineer
% x^2\x*y/ y^2
% x^3 x^2*y x*y^2 y^3
% neglect[x^4] x^3*y [x^2y^2] x*y^3 [y^4]
%
%Bending effect (Smilarly kirchoff bending plate element motion)
% w(x,y)=c9 + c10*x + c11*y + c12*x^2 + c13*xy + c14*y^2 + c15*x^3+
% c16*x^2*y + c17*x*y^2 + c18*y^3 + c19*x^3*y + c20*x*y^3
%Ru(x,y)=c1+c2*x+c3*y
%Rv(x,y)=c5+c6*x+c7*y
%Rw(x,y)=c9+c10*x+c11*y
%
%____________2.Rigidy body motion depend element geometricaly boundary conditions
%
% Membran motion Curvature
%
% ex =du/dx - w*kx kx =df^2/dx^2
% ex =dv/dy - w*ky ky =df^2/dy^2
% exy=du/dy + dv/dx - w*kxy kxy=df^2/dxdy
%
% Rigdy body motion ex=0,ey=0,exy=0
%
%du/dx-Rw(x,y)*kx=0 , du/dx=[c9+c10*x+c11*y]*kx , uo=[c9*x+c10*x^2/2+c11*xy]*kx
%dv/dy-Rw(x,y)*ky=0 , dv/dy=[c9+c10*x+c11*y]*ky , vo=[c9*y+c10*xy+1/2*c11*y^2]*ky
%du/dy+dv/dx-Rw(x,y)*2kxy=0, du/dy+d/dx*vo-Rw(x,y)*2ky=0
%
% du/dy+d/dx*[c9*x+c10*x^2/2+c11*xy]*ky-[c9+c10*x+c11*y]*2ky=0
% uo=[c9*y+c10*xy+c11*y^2/2]*2kxy-[1/2*y^2*c10]*ky
%
%du/dy+dv/dx-Rw(x,y)*2kxy=0, dv/dx+d/dy*uo-Rw(x,y)*2ky=0
%
% vo=[c9*x+1/2*c10x^2+c11*xy]*2kxy-[1/2*c11*x^2]*kx
%
%uo'=[x*kx+y*2kxy]c9 + [0.5*x^2*kx+xy*2kxy-0.5*y^2*ky]c10 + [xy*kx+y^2*kxy]c11
%vo'=[y*ky+x*2kxy]c9 + [xy*ky+x^2*kxy]c10 + [0.5*y^2*ky+2*xy*kxy-0.5*x^2*kx]c11
%
%new deplacement functions
%u(x,y):=matrix([1 , x , y , x*y , 0 , 0 , 0 , 0 , x*kx+y*2*kxy , 0.5*x^2*kx+x*y*2*kxy-0.5*y^2*ky , x*y*kx+y^2*kxy , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0]);
%v(x,y):=matrix([0 , 0 , 0 , 0 , 1 , x , y , x*y , y*ky+x*2*kxy , x*y*ky+x^2*kxy , 0.5*y^2*ky+2*x*y*kxy-0.5*x^2*kx , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0]);
%w(x,y):=matrix([0 , 0 , 0 , 0 , 0 , 0 , 0 ,0, 1 , x , y , x^2 , x*y , y^2 , x^3 , x^2*y , x*y^2 , y^3 , x^3*y , x*y^3]);
%
%##################This procedure runing maxima#############
%____________________________________________[PROCEDURE(1)]
%u(x,y):=matrix([1 , x , y , x*y , 0 , 0 , 0 , 0 , x*kx+y*2*kxy , 0.5*x^2*kx+x*y*2*kxy-0.5*y^2*ky , x*y*kx+y^2*kxy , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0]);
%v(x,y):=matrix([0 , 0 , 0 , 0 , 1 , x , y , x*y , y*ky+x*2*kxy , x*y*ky+x^2*kxy , 0.5*y^2*ky+2*x*y*kxy-0.5*x^2*kx , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0]);
%w(x,y):=matrix([0 , 0 , 0 , 0 , 0 , 0 , 0 ,0, 1 , x , y , x^2 , x*y , y^2 , x^3 , x^2*y , x*y^2 , y^3 , x^3*y , x*y^3]);
%fortran(diff(u(x,y),x));
%fortran(diff(u(x,y),y));
%fortran(diff(v(x,y),x));
%fortran(diff(v(x,y),y));
%fortran(diff(w(x,y),x));
%fortran(diff(w(x,y),y));
%fortran(diff(w(x,y),x,2));
%fortran(diff(w(x,y),y,2));
%fortran(diff(w(x,y),x,1,y,1));
%ux(x,y):=MATRIX([0,1,0,y,0,0,0,0,kx,2*kxy*y+kx*x,kx*y,0,0,0,0,0,0,0,0,0]);
%uy(x,y):=MATRIX([0,0,1,x,0,0,0,0,2*kxy,2*kxy*x-1.0*ky*y,2*kxy*y+kx*x,0,0,0,0,0,0,0,0,0]);
%vx(x,y):=MATRIX([0,0,0,0,0,1,0,y,2*kxy,ky*y+2*kxy*x,2*kxy*y-1.0*kx*x,0,0,0,0,0,0,0,0,0]);
%vy(x,y):=MATRIX([0,0,0,0,0,0,1,x,ky,ky*x,ky*y+2*kxy*x,0,0,0,0,0,0,0,0,0]);
%wx(x,y):=MATRIX([0,0,0,0,0,0,0,0,0,1,0,2*x,y,0,3*x^2,2*x*y,y^2,0,3*x^2*y,y^3]);
%wy(x,y):=MATRIX([0,0,0,0,0,0,0,0,0,0,1,0,x,2*y,0,x^2,2*x*y,3*y^2,x^3,3*x*y^2]);
%wxx(x,y):=MATRIX([0,0,0,0,0,0,0,0,0,0,0,2,0,0,6*x,2*y,0,0,6*x*y,0]);
%wyy(x,y):=MATRIX([0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,2*x,6*y,0,6*x*y]);
%wxy(x,y):=MATRIX([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,2*x,2*y,0,3*x^2,3*y^2]);
%Bound(x,y):=matrix(
%[u(-a/2,-b/2)],
%[v(-a/2,-b/2)],
%[w(-a/2,-b/2)],
%[wy(-a/2,-b/2)],
%[-wx(-a/2,-b/2)],
%[u(a/2,-b/2)],
%[v(a/2,-b/2)],
%[w(a/2,-b/2)],
%[wy(a/2,-b/2)],
%[-wx(a/2,-b/2)],
%[u(a/2,b/2)],
%[v(a/2,b/2)],
%[w(a/2,b/2)],
%[wy(a/2,b/2)],
%[-wx(a/2,b/2)],
%[u(-a/2,b/2)],
%[v(-a/2,b/2)],
%[w(-a/2,b/2)],
%[wy(-a/2,b/2)],
%[-wx(-a/2,b/2)]);
%fortmx(B,Bound(x,y));
%____________________________________________[PROCEDURE(2)]
%A(x,y):=MATRIX([1,-a/2.0,-b/2.0,a*b/4.0,0,0,0,0,-b*kxy-a*kx/2.0,-0.125*b**2*ky+a*b*kxy/2.0+0.125*a**2*kx,b**2*kxy/4.0+a*b*kx/4.0,0,0,0,0,0,0,0,0,0],
%[0,0,0,0,1,-a/2.0,-b/2.0,a*b/4.0,-b*ky/2.0-a*kxy,a*b*ky/4.0+a**2*kxy/4.0,0.125*b**2*ky+a*b*kxy/2.0-0.125*a**2*kx,0,0,0,0,0,0,0,0,0],
%[0,0,0,0,0,0,0,0,1,-a/2.0,-b/2.0,a**2/4.0,a*b/4.0,b**2/4.0,-a**3/8.0,-a**2*b/8.0,-a*b**2/8.0,-b**3/8.0,a**3*b/16.0,a*b**3/16.0],
%[0,0,0,0,0,0,0,0,0,0,1,0,-a/2.0,-b,0,a**2/4.0,a*b/2.0,3.0*b**2/4.0,-a**3/8.0,(-3.0)*a*b**2/8.0],
%[0,0,0,0,0,0,0,0,0,-1,0,a,b/2.0,0,(-3.0)*a**2/4.0,-a*b/2.0,-b**2/4.0,0,3.0*a**2*b/8.0,b**3/8.0],
%[1,a/2.0,-b/2.0,-a*b/4.0,0,0,0,0,a*kx/2.0-b*kxy,-0.125*b**2*ky-a*b*kxy/2.0+0.125*a**2*kx,b**2*kxy/4.0-a*b*kx/4.0,0,0,0,0,0,0,0,0,0],
%[0,0,0,0,1,a/2.0,-b/2.0,-a*b/4.0,a*kxy-b*ky/2.0,a**2*kxy/4.0-a*b*ky/4.0,0.125*b**2*ky-a*b*kxy/2.0-0.125*a**2*kx,0,0,0,0,0,0,0,0,0],
%[0,0,0,0,0,0,0,0,1,a/2.0,-b/2.0,a**2/4.0,-a*b/4.0,b**2/4.0,a**3/8.0,-a**2*b/8.0,a*b**2/8.0,-b**3/8.0,-a**3*b/16.0,-a*b**3/16.0],
%[0,0,0,0,0,0,0,0,0,0,1,0,a/2.0,-b,0,a**2/4.0,-a*b/2.0,3.0*b**2/4.0,a**3/8.0,3.0*a*b**2/8.0],
%[0,0,0,0,0,0,0,0,0,-1,0,-a,b/2.0,0,(-3.0)*a**2/4.0,a*b/2.0,-b**2/4.0,0,3.0*a**2*b/8.0,b**3/8.0],
%[1,a/2.0,b/2.0,a*b/4.0,0,0,0,0,b*kxy+a*kx/2.0,-0.125*b**2*ky+a*b*kxy/2.0+0.125*a**2*kx,b**2*kxy/4.0+a*b*kx/4.0,0,0,0,0,0,0,0,0,0],
%[0,0,0,0,1,a/2.0,b/2.0,a*b/4.0,b*ky/2.0+a*kxy,a*b*ky/4.0+a**2*kxy/4.0,0.125*b**2*ky+a*b*kxy/2.0-0.125*a**2*kx,0,0,0,0,0,0,0,0,0],
%[0,0,0,0,0,0,0,0,1,a/2.0,b/2.0,a**2/4.0,a*b/4.0,b**2/4.0,a**3/8.0,a**2*b/8.0,a*b**2/8.0,b**3/8.0,a**3*b/16.0,a*b**3/16.0],
%[0,0,0,0,0,0,0,0,0,0,1,0,a/2.0,b,0,a**2/4.0,a*b/2.0,3.0*b**2/4.0,a**3/8.0,3.0*a*b**2/8.0],
%[0,0,0,0,0,0,0,0,0,-1,0,-a,-b/2.0,0,(-3.0)*a**2/4.0,-a*b/2.0,-b**2/4.0,0,(-3.0)*a**2*b/8.0,-b**3/8.0],
%[1,-a/2.0,b/2.0,-a*b/4.0,0,0,0,0,b*kxy-a*kx/2.0,-0.125*b**2*ky-a*b*kxy/2.0+0.125*a**2*kx,b**2*kxy/4.0-a*b*kx/4.0,0,0,0,0,0,0,0,0,0],
%[0,0,0,0,1,-a/2.0,b/2.0,-a*b/4.0,b*ky/2.0-a*kxy,a**2*kxy/4.0-a*b*ky/4.0,0.125*b**2*ky-a*b*kxy/2.0-0.125*a**2*kx,0,0,0,0,0,0,0,0,0],
%[0,0,0,0,0,0,0,0,1,-a/2.0,b/2.0,a**2/4.0,-a*b/4.0,b**2/4.0,-a**3/8.0,a**2*b/8.0,-a*b**2/8.0,b**3/8.0,-a**3*b/16.0,-a*b**3/16.0],
%[0,0,0,0,0,0,0,0,0,0,1,0,-a/2.0,b,0,a**2/4.0,-a*b/2.0,3.0*b**2/4.0,-a**3/8.0,(-3.0)*a*b**2/8.0],
%[0,0,0,0,0,0,0,0,0,-1,0,a,-b/2.0,0,(-3.0)*a**2/4.0,a*b/2.0,-b**2/4.0,0,(-3.0)*a**2*b/8.0,-b**3/8.0]);
%%____________________________________________[PROCEDURE(3)]
%Der(x,y):=matrix(
%[ux(x,y)-w(x,y)*kx],
%[vy(x,y)-w(x,y)*ky],
%[uy(x,y)+vx(x,y)-w(x,y)*2*kxy],
%[-wxx(x,y)],
%[-wyy(x,y)],
%[-2*wxy(x,y)]);
%fortmx(Drv,Der(x,y));
%Drv(x,y):=MATRIX([0,1,0,y,0,0,0,0,0,2*kxy*y,0,-kx*x**2,-kx*x*y,-kx*y**2,-kx*x**3,-kx*x**2*y,-kx*x*y**2,-kx*y**3,-kx*x**3*y,-kx*x*y**3],
%[0,0,0,0,0,0,1,x,0,0,2*kxy*x,-ky*x**2,-ky*x*y,-ky*y**2,-ky*x**3,-ky*x**2*y,-ky*x*y**2,-ky*y**3,-ky*x**3*y,-ky*x*y**3],
%[0,0,1,x,0,1,0,y,2*kxy,2*kxy*x,2*kxy*y,-2*kxy*x**2,-2*kxy*x*y,-2*kxy*y**2,-2*kxy*x**3,-2*kxy*x**2*y,-2*kxy*x*y**2,-2*kxy*y**3,-2*kxy*x**3*y,-2*kxy*x*y**3],
%[0,0,0,0,0,0,0,0,0,0,0,-2,0,0,-6*x,-2*y,0,0,-6*x*y,0],
%[0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,-2*x,-6*y,0,-6*x*y],
%[0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,-4*x,-4*y,0,-6*x**2,-6*y**2]);
%C(x,y):=matrix([KSI1,0,0,0,0,0],
% [0,KSI1,0,0,0,0],
% [0,0,KSI1,0,0,0],
% [0,0,0,KSI2,0,0],
% [0,0,0,0,KSI2,0],
% [0,0,0,0,0,KSI2]).matrix([1 ,v ,0 ,0 ,0 ,0 ],
% [v ,1 ,0 ,0 ,0 ,0 ],
% [0 ,0 ,(1-v)/2 ,0 ,0 ,0 ],
% [0 ,0 ,0 ,1 ,v ,0 ],
% [0 ,0 ,0 ,v ,1 ,0 ],
% [0 ,0 ,0 ,0 ,0,(1-v)/2]);
%Kd(x,y):=defint(defint(Transpose(Drv(x,y)).C(x,y).Drv(x,y),x,-a/2,a/2),y,-b/2,b/2);
%fortmx(Kd,Kd(x,y));
%__________________________________________________________________________
%| 4NODE-20DOF CURVED PARAMETRIC SHELL FEM [A.] (C)(R) |
%|######################################## |
%| Element family:Low proclivity curved shell |
%| Shape function:Workable(geometric boundary condition)|
%| Stifness :Topology + Accumulate method |
%| Theoretic app.Total potantial enegry second theorem |
%| Solve Equation:Cholesky-[L][D][u] |
%| |
%|_______________________________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. |
%|_________________________________________________________________________|
clc
clear
%##################INPUT DATA###############
%####################### MATERIALS PROPERTIES
E=2E6; %shell material elasticity constant;
v=0.3; %shell material poission ratio;
%#######################
%####################### AUTOMATIC MESH DATA
meshx = 2; %"x" axis direction mesh number
meshy = 2; %"y" axis direction mesh number
totala = 2.00; %"x" axis direction total length
totalb = 2.00; %"y" axis direction total width
th = 0.10; %thickness shell element;
%Notice: Shell`s surface function f(x,y):=-1/10y^2-1/10x^2+5
%kx=d^2/dx^2*f(x,y)=-2/10
%ky=d^2/dy^2*f(x,y)=-2/10
%kxy=d^2/dxdy*f(x,y)=0
kgx = -2/10; %x axis direction`s curvature;
kgy = -2/10; %y axis direction`s curvature;
kgxy = 0; %xy axis direction`s curvature;
%########################
gridspacex = totala/meshx; % along to "x" axis direction mesh space
gridspacey = totalb/meshy; % along to "y" axis direction mesh space
%##################### MESH SURFACE FUNCTIONS
%fshell=inline('-1/10*(x)^2-1/10*(y)^2+5');
%#####################
%Mesh error control
if meshx <1 || meshy <1;
display('Open dimension system mesh and this mesh value > 1')
% [meshx meshy]
error('mesh values')
else if meshx >10 || meshy>10;
display('system analysis mesh is big')
end
end
%_____________________________________GLOBAL COORDINATES
dim=1;
for sut=0:meshy
for sat=0:meshx %Element middle node reference
x=sat*gridspacex-meshx*gridspacex/2;
y=sut*gridspacey-meshy*gridspacey/2;
Cor(dim,1)=x; %"x" global coordinate value
Cor(dim,2)=y; %"y" global coordinate value
Cor(dim,3)=-1/10*(x)^2-1/10*(y)^2+5;%"z" global coordinates value
%Cor(dim,3)=fshell(x,y);%"z" global coordinates value
dim=dim+1;
end
end%___________________________________________|
Number=size(Cor);clear x y;
Node=Number(1); %Elements Number
%____________________________________POSITION MATRIX
no=1; %Pos(Shell Element No,:)=[Shell Element Node Number]
for i=1:meshy+1;
for j=1:meshx;
if i < meshy+1 && j<meshx+1;
Pos(no,1)=(meshx+1)*(i-1)+j;
Pos(no,2)=(meshx+1)*(i-1)+j+1;
Pos(no,3)=(meshx+1)*(i)+j+1;
Pos(no,4)=(meshx+1)*(i)+j;
end
no=no+1;
end
end%________________________________________|
Number=size(Pos);
No=Number(1); %Plane Elements Number
%Shell`s Node Topology
for i=1:Node;
Re(i,:)=[1 1 1 1 1 1];
end
%Re(Node number,:)=[ux vy wz theta beta]
%############################SYSTEM SUPPORTS
Re(1,:)=[0 0 0 0 0 0];
Re(3,:)=[0 0 0 0 0 0];
Re(7,:)=[0 0 0 0 0 0];
Re(9,:)=[0 0 0 0 0 0];
%#############################
Topology=Re;
%_________________________________ACCUMULATE
Number=size(Re);
Nom=Number(2); %Element node d.o.f Number
clear Number
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;
%_________________________________MODAL DISPLACEMENTS
for i=1:No
R(i,:)=[Re(Pos(i,1),:) Re(Pos(i,2),:) Re(Pos(i,3),:) Re(Pos(i,4),:)];
end
%for i=1:No;
% for j=1:Nom;
% R(i,j+Nom*0) = Re(Pos(i,1),j);
% R(i,j+Nom*1) = Re(Pos(i,2),j);
% R(i,j+Nom*2) = Re(Pos(i,3),j);
% R(i,j+Nom*3) = Re(Pos(i,4),j);
% end
%end%_________________________________________|
%######################SYSTEM LOADS
%P(Re(Node number, Freedoom number))=[u v w theta beta]
P(Item)=0;
P(Re(2,3))=-5;
P(Re(4,3))=-5;
P(Re(6,3))=-5;
P(Re(8,3))=-5;
%#######################
%__________________________________________________SHELL_SYSTEM_ANALYSE
%_____________________________LOCAL COORDINATES
%Element global node cartesian coordiantes moving local nodes
for Element=1:No;
px(Element,:)=[Cor(Pos(Element,1),1)...
Cor(Pos(Element,2),1)...
Cor(Pos(Element,3),1)...
Cor(Pos(Element,4),1)];
py(Element,:)=[Cor(Pos(Element,1),2)...
Cor(Pos(Element,2),2)...
Cor(Pos(Element,3),2)...
Cor(Pos(Element,4),2)];
pz(Element,:)=[Cor(Pos(Element,1),3)...
Cor(Pos(Element,2),3)...
Cor(Pos(Element,3),3)...
Cor(Pos(Element,4),3)];
end%____________________________|
%________________________________ELEMENT DIRECTION COSINESES
for Element=1:No;
pxo5=(px(Element,1)+px(Element,2))/2;
pyo5=(py(Element,1)+py(Element,2))/2;
pzo5=(pz(Element,1)+pz(Element,2))/2;
pxo6=(px(Element,2)+px(Element,3))/2;
pyo6=(py(Element,2)+py(Element,3))/2;
pzo6=(pz(Element,2)+pz(Element,3))/2;
pxo7=(px(Element,3)+px(Element,4))/2;
pyo7=(py(Element,3)+py(Element,4))/2;
pzo7=(pz(Element,3)+pz(Element,4))/2;
pxo8=(px(Element,1)+px(Element,4))/2;
pyo8=(py(Element,1)+py(Element,4))/2;
pzo8=(pz(Element,1)+pz(Element,4))/2;
Lu1=sqrt((pxo6-pxo8)^2+(pyo6-pyo8)^2+(pzo6-pzo8)^2);
directvectorcos1(Element,:)=[(pxo6-pxo8)/Lu1 ...
(pyo6-pyo8)/Lu1 ...
(pzo6-pzo8)/Lu1];
Lu2=sqrt((pxo7-pxo5)^2+(pyo7-pyo5)^2+(pzo7-pzo5)^2);
directvectorcos2(Element,:)=[(pxo7-pxo5)/Lu2 ...
(pyo7-pyo5)/Lu2 ...
(pzo7-pzo5)/Lu2];
%__________Shell middle surface thickness direction vector [V3](e3)
%__________Vectorel Multiply "x" direction and "y" position vector{e3}=e1*e2
%__[l1(i)+m1(j)+n1(k)]*[l2(i)+m2(j)+n2(k)]=
%__[l1*l2(i*i)+l1*m2(i*j)+l1*n2(i*k)]+
%__[m1*l2(j*i)+m1*m2(j*j)+m1*n2(j*k)]+
%__[n1*l2(k*i)+n1*m2(k*j)+n1*n2(k*k)]
%__[l1*m2(k)+l1*n2(-j)+m1*l2(-k)+m1*n2(i)+n1*l2(j)+n1*m2(-i)
%__i[m1*n2-n1*m2]+j[n1*l2-l1*n2]+k[l1*m2-m1*l2]
l1=directvectorcos1(Element,1); l2=directvectorcos2(Element,1);
m1=directvectorcos1(Element,2); m2=directvectorcos2(Element,2);
n1=directvectorcos1(Element,3); n2=directvectorcos2(Element,3);
Lu3=sqrt((m1*n2 - m2*n1)^2+ ...
(n1*l2 - n2*l1)^2+ ...
(l1*m2 - m1*l2)^2);
directvectorcos3(Element,:)=[(m1*n2-n1*m2)/Lu3 ...
(n1*l2-l1*n2)/Lu3 ...
(l1*m2-m1*l2)/Lu3];
end%________________________|
clear Lu1 Lu2 Lu3 pzo5 pzo6 pzo7 pzo8 pyo5 pyo6 pyo7 pyo8
clear pxo5 pxo6 pxo7 pxo8
%___________________________SPECIAL AXIS DIRECTION COSINESES
%This [lx,ly,lz] direction cosines by shell special axis
for Element=1:No;
pxn=[px(Element,1) px(Element,2) px(Element,3) px(Element,4)];
pyn=[py(Element,1) py(Element,2) py(Element,3) py(Element,4)];
pzn=[pz(Element,1) pz(Element,2) pz(Element,3) pz(Element,4)];
for j=1:4
lx=directvectorcos1(Element,1); ly=directvectorcos2(Element,1);
mx=directvectorcos1(Element,2); my=directvectorcos2(Element,2);
%Per unit shell`s elements node rotation [Jacobian Trans.]
%Notice this rotation angles is values get radii
%Shell surface function
% f(x,y)g = -1/10*(x)^2-1/10*(y)^2+5
% [df/dx]g = -2/10*(x)^1*1
% [df/dy]g = -2/10*(y)^1*1
%Moving shell global nodes cartesian node`s slopes , element special axis direction
%
% shellspecialaxis(xg) = lx*xs + ly*ys + lz*zs [l,m,n] special-->global axis dir. cos.
% shellspecialaxis(yg) = mx*xs + my*ys + mz*zs
% shellspecialaxis(zg) = nx*xs + ny*ys + nz*zs
%
% [df/dxs] [dxg/dxs dyg/dxs] [df/dxg]
% [ ] = [ ]* [ ]
% [df/dys] [dxg/dys dyg/dys] [df/dyg]
% __Jacobian transform__
%
% [df/dxs] [ lx mx ] [-2/10*(x)^1*1]
% [ ] = [ ]* [ ]
% [df/dys] [ ly my ] [-2/10*(y)^1*1]
%################## SHELL PARTIAL DERIVATIVES
tanalfax=(lx*-2*pxn(j)+mx*-2*pyn(j))/10;
tanalfay=(ly*-2*pxn(j)+my*-2*pyn(j))/10;
%##################
alfax=atan(tanalfax);
alfay=atan(tanalfay);
% Ctr2(:,:,Element)=[cos(alfax) 0 0;
% 0 cos(alfay) 0;
% sin(alfax) sin(alfay) 1];
directcurvecosa(:,j,Element)=[cos(alfax)
0
sin(alfax)];
directcurvecosb(:,j,Element)=[ 0
cos(alfay)
sin(alfay)];
end
end%___________________________|
clear tanalfax tanalfay sayman pxn pyn pzn px py pz
%_______________________ELEMENT TRANSFORMATION MATRIX
zero =[0 0 0 0 0 0;
0 0 0 0 0 0;
0 0 0 0 0 0;
0 0 0 0 0 0;
0 0 0 0 0 0;
0 0 0 0 0 0];
zero1 =[0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0 ];
%Shell element transform special axis to system global axis
for Element=1:No;
%Element axises for global system direction cosineses
l1=directvectorcos1(Element,1);
m1=directvectorcos1(Element,2);
n1=directvectorcos1(Element,3);
l2=directvectorcos2(Element,1);
m2=directvectorcos2(Element,2);
n2=directvectorcos2(Element,3);
l3=directvectorcos3(Element,1);
m3=directvectorcos3(Element,2);
n3=directvectorcos3(Element,3);
%Element main local axis to global axis convert matrix
Tc =[ l1 l2 l3 0 0 0 ;
m1 m2 m3 0 0 0 ;
n1 n2 n3 0 0 0 ;
0 0 0 l1 l2 l3 ;
0 0 0 m1 m2 m3 ;
0 0 0 n1 n2 n3];
%Element per node has curved transform direction cosineses.
la1=directcurvecosa(1,1,Element); lb1=directcurvecosb(1,1,Element);
ma1=directcurvecosa(2,1,Element); mb1=directcurvecosb(2,1,Element);
na1=directcurvecosa(3,1,Element); nb1=directcurvecosb(3,1,Element);
la2=directcurvecosa(1,2,Element); lb2=directcurvecosb(1,2,Element);
ma2=directcurvecosa(2,2,Element); mb2=directcurvecosb(2,2,Element);
na2=directcurvecosa(3,2,Element); nb2=directcurvecosb(3,2,Element);
la3=directcurvecosa(1,3,Element); lb3=directcurvecosb(1,3,Element);
ma3=directcurvecosa(2,3,Element); mb3=directcurvecosb(2,3,Element);
na3=directcurvecosa(3,3,Element); nb3=directcurvecosb(3,3,Element);
la4=directcurvecosa(1,4,Element); lb4=directcurvecosb(1,4,Element);
ma4=directcurvecosa(2,4,Element); mb4=directcurvecosb(2,4,Element);
na4=directcurvecosa(3,4,Element); nb4=directcurvecosb(3,4,Element);
% [Tson]:Local to Global Tranformation Matrix
%_________Per Unit Node transformation matrix
Tci=[ la1 lb1 0 0 0 ; %(1) Transformation depend node slopes
ma1 mb1 0 0 0 ;
na1 nb1 1 0 0 ;
0 0 0 la1 lb1 ;
0 0 0 ma1 mb1 ;
0 0 0 na1 nb1 ];
Tcj=[ la2 lb2 0 0 0 ;%(2) Transformation depend node slopes
ma2 mb2 0 0 0 ;
na2 nb2 1 0 0 ;
0 0 0 la2 lb2 ;
0 0 0 ma2 mb2 ;
0 0 0 na2 nb2 ];
Tck=[ la3 lb3 0 0 0 ;%(3) Transformation depend node slopes
ma3 mb3 0 0 0 ;
na3 nb3 1 0 0 ;
0 0 0 la3 lb3 ;
0 0 0 ma3 mb3 ;
0 0 0 na3 nb3 ];
Tcl=[ la4 lb4 0 0 0 ;%(4) Transformation depend node slopes
ma4 mb4 0 0 0 ;
na4 nb4 1 0 0 ;
0 0 0 la4 lb4 ;
0 0 0 ma4 mb4 ;
0 0 0 na4 nb4 ];
%Element local axis convert system global axis
Ts =[Tc , zero , zero , zero ;
zero , Tc , zero , zero ;
zero , zero , Tc , zero ;
zero , zero , zero , Tc ];
%Curved node convertion matrix
Tson2 =[ Tci , zero1 , zero1 , zero1 ;
zero1 , Tcj , zero1 , zero1 ;
zero1 , zero1 , Tck , zero1 ;
zero1 , zero1 , zero1 , Tcl ];
Tson(:,:,Element)=Ts*Tson2;
end%___________________________|
clear Tci Tcj Tck Tcl Ts Tson2 Tc l1 l2 l3 m1 m2 m3
clear la1 la2 la3 la4 lb1 lb2 lb3 lb4 zero zero1
clear na1 na2 na3 na4 nb1 nb2 nb3 nb4 n1 n2 n3
clear mb1 mb2 mb3 mb4 ma1 ma2 ma3 ma4
%_______________________________BOUNDARY CONDITIONS COEFFICIENT MATRIX
KSI1= E*th/(1-v^2);%Elasticity constants
KSI2= E*th^3/(12*(1-v^2));
%Open dimension as shell element local axis stiffness matrix
Kd(20,20)=0;
for s=1:No;
lx=directvectorcos1(Element,1);
mx=directvectorcos1(Element,2);
nx=directvectorcos1(Element,3);
ly=directvectorcos2(Element,1);
my=directvectorcos2(Element,2);
ny=directvectorcos2(Element,3);
%Element global axis curvatures transform to local axis
kx =lx*kgx + mx*kgy + nx*kgxy;
ky =ly*kgx + my*kgy + ny*kgxy;
kxy=0;%lz*kgx + mz*kgy + nz*kgxy;(neglected)[system symmetric]
%[kx ky kxy]
a=gridspacex;
b=gridspacey;
%________________________________________________________[PROCEDURE(3) OUTPUT]
Kd(2,2,s) = a*b*KSI1;
Kd(2,7,s) = a*b*KSI1*v;
Kd(2,12,s) = -b*(a^3*KSI1*ky*v+a^3*KSI1*kx)/12.0;
Kd(2,14,s) = -a*(b^3*KSI1*ky*v+b^3*KSI1*kx)/12.0;
Kd(3,3,s) = a*b*KSI1*(1-v)/2.0;
Kd(3,6,s) = a*b*KSI1*(1-v)/2.0;
Kd(3,9,s) = a*b*KSI1*kxy*(1-v);
Kd(3,12,s) = -a^3*b*KSI1*kxy*(1-v)/12.0;
Kd(3,14,s) = -a*b^3*KSI1*kxy*(1-v)/12.0;
Kd(4,4,s) = (-a^3*b*KSI1*v-(-2*a*b^3-a^3*b)*KSI1)/24.0;
Kd(4,10,s) = (-a^3*b*KSI1*kxy*v-(-2*a*b^3-a^3*b)*KSI1*kxy)/12.0;
Kd(4,15,s) = (a^5*b*KSI1*kxy*v-a^5*b*KSI1*kxy)/80.0;
Kd(4,16,s) = -(a^3*b^3*KSI1*ky*v+a^3*b^3*KSI1*kx)/144.0;
Kd(4,17,s) = (a^3*b^3*KSI1*kxy*v-a^3*b^3*KSI1*kxy)/144.0;
Kd(4,18,s) = -(a*b^5*KSI1*ky*v+a*b^5*KSI1*kx)/80.0;
Kd(6,3,s) = a*b*KSI1*(1-v)/2.0;
Kd(6,6,s) = a*b*KSI1*(1-v)/2.0;
Kd(6,9,s) = a*b*KSI1*kxy*(1-v);
Kd(6,12,s) = -a^3*b*KSI1*kxy*(1-v)/12.0;
Kd(6,14,s) = -a*b^3*KSI1*kxy*(1-v)/12.0;
Kd(7,2,s) = a*b*KSI1*v;
Kd(7,7,s) = a*b*KSI1;
Kd(7,12,s) = -b*(a^3*KSI1*kx*v+a^3*KSI1*ky)/12.0;
Kd(7,14,s) = -a*(b^3*KSI1*kx*v+b^3*KSI1*ky)/12.0;
Kd(8,8,s) = -(a*b^3*KSI1*v+(-a*b^3-2*a^3*b)*KSI1)/24.0;
Kd(8,11,s) = -(a*b^3*KSI1*kxy*v+(-a*b^3-2*a^3*b)*KSI1*kxy)/12.0;
Kd(8,15,s) = -(a^5*b*KSI1*kx*v+a^5*b*KSI1*ky)/80.0;
Kd(8,16,s) = (a^3*b^3*KSI1*kxy*v-a^3*b^3*KSI1*kxy)/144.0;
Kd(8,17,s) = -(a^3*b^3*KSI1*kx*v+a^3*b^3*KSI1*ky)/144.0;
Kd(8,18,s) = (a*b^5*KSI1*kxy*v-a*b^5*KSI1*kxy)/80.0;
Kd(9,3,s) = a*b*KSI1*kxy*(1-v);
Kd(9,6,s) = a*b*KSI1*kxy*(1-v);
Kd(9,9,s) = 2*a*b*KSI1*kxy^2*(1-v);
Kd(9,12,s) = -a^3*b*KSI1*kxy^2*(1-v)/6.0;
Kd(9,14,s) = -a*b^3*KSI1*kxy^2*(1-v)/6.0;
Kd(10,4,s) = (-a^3*b*KSI1*kxy*v-(-2*a*b^3-a^3*b)*KSI1*kxy)/12.0;
Kd(10,10,s)= (-a^3*b*KSI1*kxy^2*v-(-2*a*b^3-a^3*b)*KSI1*kxy^2)/6.0;
Kd(10,15,s)= (a^5*b*KSI1*kxy^2*v-a^5*b*KSI1*kxy^2)/40.0;
Kd(10,16,s)= -(a^3*b^3*KSI1*kxy*ky*v+a^3*b^3*KSI1*kx*kxy)/72.0;
Kd(10,17,s)= (a^3*b^3*KSI1*kxy^2*v-a^3*b^3*KSI1*kxy^2)/72.0;
Kd(10,18,s)= -(a*b^5*KSI1*kxy*ky*v+a*b^5*KSI1*kx*kxy)/40.0;
Kd(11,8,s) = -(a*b^3*KSI1*kxy*v+(-a*b^3-2*a^3*b)*KSI1*kxy)/12.0;
Kd(11,11,s)= -(a*b^3*KSI1*kxy^2*v+(-a*b^3-2*a^3*b)*KSI1*kxy^2)/6.0;
Kd(11,15,s)= -(a^5*b*KSI1*kx*kxy*v+a^5*b*KSI1*kxy*ky)/40.0;
Kd(11,16,s)= (a^3*b^3*KSI1*kxy^2*v-a^3*b^3*KSI1*kxy^2)/72.0;
Kd(11,17,s)= -(a^3*b^3*KSI1*kx*kxy*v+a^3*b^3*KSI1*kxy*ky)/72.0;
Kd(11,18,s)= (a*b^5*KSI1*kxy^2*v-a*b^5*KSI1*kxy^2)/40.0;
Kd(12,2,s) = -b*(a^3*KSI1*ky*v+a^3*KSI1*kx)/12.0;
Kd(12,3,s) = -a^3*b*KSI1*kxy*(1-v)/12.0;
Kd(12,6,s) = -a^3*b*KSI1*kxy*(1-v)/12.0;
Kd(12,7,s) = -b*(a^3*KSI1*kx*v+a^3*KSI1*ky)/12.0;
Kd(12,9,s) = -a^3*b*KSI1*kxy^2*(1-v)/6.0;
Kd(12,12,s)= b*((2*a^5*KSI1*kx*ky-2*a^5*KSI1*kxy^2)*v+a^5*KSI1*ky^2+2*a^5*KSI1*kxy^2+a^5*KSI1*kx^2+320*a*KSI2)/80.0;
Kd(12,14,s)= ((2*a^3*b^3*KSI1*kx*ky-2*a^3*b^3*KSI1*kxy^2+576*a*b*KSI2)*v+a^3*b^3*KSI1*ky^2+2*a^3*b^3*KSI1*kxy^2+a^3*b^3*KSI1*kx^2)/144.0;
Kd(13,13,s)= ((2*a^3*b^3*KSI1*kx*ky-2*a^3*b^3*KSI1*kxy^2-288*a*b*KSI2)*v+a^3*b^3*KSI1*ky^2+2*a^3*b^3*KSI1*kxy^2+a^3*b^3*KSI1*kx^2+288*a*b*KSI2)/144.0;
Kd(13,19,s)= ((2*a^5*b^3*KSI1*kx*ky-2*a^5*b^3*KSI1*kxy^2-480*a^3*b*KSI2)*v+a^5*b^3*KSI1*ky^2+2*a^5*b^3*KSI1*kxy^2+a^5*b^3*KSI1*kx^2+480*a^3*b*KSI2)/960.0;
Kd(13,20,s)= ((2*a^3*b^5*KSI1*kx*ky-2*a^3*b^5*KSI1*kxy^2-480*a*b^3*KSI2)*v+a^3*b^5*KSI1*ky^2+2*a^3*b^5*KSI1*kxy^2+a^3*b^5*KSI1*kx^2+480*a*b^3*KSI2)/960.0;
Kd(14,2,s) = -a*(b^3*KSI1*ky*v+b^3*KSI1*kx)/12.0;
Kd(14,3,s) = -a*b^3*KSI1*kxy*(1-v)/12.0;
Kd(14,6,s) = -a*b^3*KSI1*kxy*(1-v)/12.0;
Kd(14,7,s) = -a*(b^3*KSI1*kx*v+b^3*KSI1*ky)/12.0;
Kd(14,9,s) = -a*b^3*KSI1*kxy^2*(1-v)/6.0;
Kd(14,12,s)= ((2*a^3*b^3*KSI1*kx*ky-2*a^3*b^3*KSI1*kxy^2+576*a*b*KSI2)*v+a^3*b^3*KSI1*ky^2+2*a^3*b^3*KSI1*kxy^2+a^3*b^3*KSI1*kx^2)/144.0;
Kd(14,14,s)= a*((2*b^5*KSI1*kx*ky-2*b^5*KSI1*kxy^2)*v+b^5*KSI1*ky^2+2*b^5*KSI1*kxy^2+b^5*KSI1*kx^2+320*b*KSI2)/80.0;
Kd(15,4,s) = (a^5*b*KSI1*kxy*v-a^5*b*KSI1*kxy)/80.0;
Kd(15,8,s) = -(a^5*b*KSI1*kx*v+a^5*b*KSI1*ky)/80.0;
Kd(15,10,s)= (a^5*b*KSI1*kxy^2*v-a^5*b*KSI1*kxy^2)/40.0;
Kd(15,11,s)= -(a^5*b*KSI1*kx*kxy*v+a^5*b*KSI1*kxy*ky)/40.0;
Kd(15,15,s)= b*((2*a^7*KSI1*kx*ky-2*a^7*KSI1*kxy^2)*v+a^7*KSI1*ky^2+2*a^7*KSI1*kxy^2+a^7*KSI1*kx^2+1344*a^3*KSI2)/448.0;
Kd(15,17,s)= ((2*a^5*b^3*KSI1*kx*ky-2*a^5*b^3*KSI1*kxy^2+960*a^3*b*KSI2)*v+a^5*b^3*KSI1*ky^2+2*a^5*b^3*KSI1*kxy^2+a^5*b^3*KSI1*kx^2)/960.0;
Kd(16,4,s) = -(a^3*b^3*KSI1*ky*v+a^3*b^3*KSI1*kx)/144.0;
Kd(16,8,s) = (a^3*b^3*KSI1*kxy*v-a^3*b^3*KSI1*kxy)/144.0;
Kd(16,10,s)= -(a^3*b^3*KSI1*kxy*ky*v+a^3*b^3*KSI1*kx*kxy)/72.0;
Kd(16,11,s)= (a^3*b^3*KSI1*kxy^2*v-a^3*b^3*KSI1*kxy^2)/72.0;
Kd(16,16,s)= ((2*a^5*b^3*KSI1*kx*ky-2*a^5*b^3*KSI1*kxy^2-640*a^3*b*KSI2)*v+a^5*b^3*KSI1*ky^2+2*a^5*b^3*KSI1*kxy^2+a^5*b^3*KSI1*kx^2+(320*a*b^3+640*a^3*b)*KSI2)/960.0;
Kd(16,18,s)= ((2*a^3*b^5*KSI1*kx*ky-2*a^3*b^5*KSI1*kxy^2+960*a*b^3*KSI2)*v+a^3*b^5*KSI1*ky^2+2*a^3*b^5*KSI1*kxy^2+a^3*b^5*KSI1*kx^2)/960.0;
Kd(17,4,s) = (a^3*b^3*KSI1*kxy*v-a^3*b^3*KSI1*kxy)/144.0;
Kd(17,8,s) = -(a^3*b^3*KSI1*kx*v+a^3*b^3*KSI1*ky)/144.0;
Kd(17,10,s)= (a^3*b^3*KSI1*kxy^2*v-a^3*b^3*KSI1*kxy^2)/72.0;
Kd(17,11,s)= -(a^3*b^3*KSI1*kx*kxy*v+a^3*b^3*KSI1*kxy*ky)/72.0;
Kd(17,15,s)= ((2*a^5*b^3*KSI1*kx*ky-2*a^5*b^3*KSI1*kxy^2+960*a^3*b*KSI2)*v+a^5*b^3*KSI1*ky^2+2*a^5*b^3*KSI1*kxy^2+a^5*b^3*KSI1*kx^2)/960.0;
Kd(17,17,s)= ((2*a^3*b^5*KSI1*kx*ky-2*a^3*b^5*KSI1*kxy^2-640*a*b^3*KSI2)*v+a^3*b^5*KSI1*ky^2+2*a^3*b^5*KSI1*kxy^2+a^3*b^5*KSI1*kx^2+(640*a*b^3+320*a^3*b)*KSI2)/960.0;
Kd(18,4,s) = -(a*b^5*KSI1*ky*v+a*b^5*KSI1*kx)/80.0;
Kd(18,8,s) = (a*b^5*KSI1*kxy*v-a*b^5*KSI1*kxy)/80.0;
Kd(18,10,s)= -(a*b^5*KSI1*kxy*ky*v+a*b^5*KSI1*kx*kxy)/40.0;
Kd(18,11,s)= (a*b^5*KSI1*kxy^2*v-a*b^5*KSI1*kxy^2)/40.0;
Kd(18,16,s)= ((2*a^3*b^5*KSI1*kx*ky-2*a^3*b^5*KSI1*kxy^2+960*a*b^3*KSI2)*v+a^3*b^5*KSI1*ky^2+2*a^3*b^5*KSI1*kxy^2+a^3*b^5*KSI1*kx^2)/960.0;
Kd(18,18,s)= a*((2*b^7*KSI1*kx*ky-2*b^7*KSI1*kxy^2)*v+b^7*KSI1*ky^2+2*b^7*KSI1*kxy^2+b^7*KSI1*kx^2+1344*b^3*KSI2)/448.0;
Kd(19,13,s)= ((2*a^5*b^3*KSI1*kx*ky-2*a^5*b^3*KSI1*kxy^2-480*a^3*b*KSI2)*v+a^5*b^3*KSI1*ky^2+2*a^5*b^3*KSI1*kxy^2+a^5*b^3*KSI1*kx^2+480*a^3*b*KSI2)/960.0;
Kd(19,19,s)= ((10*a^7*b^3*KSI1*kx*ky-10*a^7*b^3*KSI1*kxy^2-6048*a^5*b*KSI2)*v+5*a^7*b^3*KSI1*ky^2+10*a^7*b^3*KSI1*kxy^2+5*a^7*b^3*KSI1*kx^2+(6720*a^3*b^3+6048*a^5*b)*KSI2)/26880.0;
Kd(19,20,s)= ((2*a^5*b^5*KSI1*kx*ky-2*a^5*b^5*KSI1*kxy^2+800*a^3*b^3*KSI2)*v+a^5*b^5*KSI1*ky^2+2*a^5*b^5*KSI1*kxy^2+a^5*b^5*KSI1*kx^2+800*a^3*b^3*KSI2)/6400.0;
Kd(20,13,s)= ((2*a^3*b^5*KSI1*kx*ky-2*a^3*b^5*KSI1*kxy^2-480*a*b^3*KSI2)*v+a^3*b^5*KSI1*ky^2+2*a^3*b^5*KSI1*kxy^2+a^3*b^5*KSI1*kx^2+480*a*b^3*KSI2)/960.0;
Kd(20,19,s)= ((2*a^5*b^5*KSI1*kx*ky-2*a^5*b^5*KSI1*kxy^2+800*a^3*b^3*KSI2)*v+a^5*b^5*KSI1*ky^2+2*a^5*b^5*KSI1*kxy^2+a^5*b^5*KSI1*kx^2+800*a^3*b^3*KSI2)/6400.0;
Kd(20,20,s)= ((10*a^3*b^7*KSI1*kx*ky-10*a^3*b^7*KSI1*kxy^2-6048*a*b^5*KSI2)*v+5*a^3*b^7*KSI1*ky^2+10*a^3*b^7*KSI1*kxy^2+5*a^3*b^7*KSI1*kx^2+(6048*a*b^5+6720*a^3*b^3)*KSI2)/26880.0;
%____________________________________________________[PROCEDURE(2)]__[OUTPUT]
%Element connection matrix
A(:,:,s)=[ 1 -a/2.0 -b/2.0 a*b/4.0 0 0 0 0 -b*kxy-a*kx/2.0 -0.125*b^2*ky+a*b*kxy/2.0+0.125*a^2*kx b^2*kxy/4.0+a*b*kx/4.0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 -a/2.0 -b/2.0 a*b/4.0 -b*ky/2.0-a*kxy a*b*ky/4.0+a^2*kxy/4.0 0.125*b^2*ky+a*b*kxy/2.0-0.125*a^2*kx 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 -a/2.0 -b/2.0 a^2/4.0 a*b/4.0 b^2/4.0 -a^3/8.0 -a^2*b/8.0 -a*b^2/8.0 -b^3/8.0 a^3*b/16.0 a*b^3/16.0;
0 0 0 0 0 0 0 0 0 0 1 0 -a/2.0 -b 0 a^2/4.0 a*b/2.0 3.0*b^2/4.0 -a^3/8.0 (-3.0)*a*b^2/8.0;
0 0 0 0 0 0 0 0 0 -1 0 a b/2.0 0 (-3.0)*a^2/4.0 -a*b/2.0 -b^2/4.0 0 3.0*a^2*b/8.0 b^3/8.0;
1 a/2.0 -b/2.0 -a*b/4.0 0 0 0 0 a*kx/2.0-b*kxy -0.125*b^2*ky-a*b*kxy/2.0+0.125*a^2*kx b^2*kxy/4.0-a*b*kx/4.0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 a/2.0 -b/2.0 -a*b/4.0 a*kxy-b*ky/2.0 a^2*kxy/4.0-a*b*ky/4.0 0.125*b^2*ky-a*b*kxy/2.0-0.125*a^2*kx 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 a/2.0 -b/2.0 a^2/4.0 -a*b/4.0 b^2/4.0 a^3/8.0 -a^2*b/8.0 a*b^2/8.0 -b^3/8.0 -a^3*b/16.0 -a*b^3/16.0;
0 0 0 0 0 0 0 0 0 0 1 0 a/2.0 -b 0 a^2/4.0 -a*b/2.0 3.0*b^2/4.0 a^3/8.0 3.0*a*b^2/8.0;
0 0 0 0 0 0 0 0 0 -1 0 -a b/2.0 0 (-3.0)*a^2/4.0 a*b/2.0 -b^2/4.0 0 3.0*a^2*b/8.0 b^3/8.0;
1 a/2.0 b/2.0 a*b/4.0 0 0 0 0 b*kxy+a*kx/2.0 -0.125*b^2*ky+a*b*kxy/2.0+0.125*a^2*kx b^2*kxy/4.0+a*b*kx/4.0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 a/2.0 b/2.0 a*b/4.0 b*ky/2.0+a*kxy a*b*ky/4.0+a^2*kxy/4.0 0.125*b^2*ky+a*b*kxy/2.0-0.125*a^2*kx 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 a/2.0 b/2.0 a^2/4.0 a*b/4.0 b^2/4.0 a^3/8.0 a^2*b/8.0 a*b^2/8.0 b^3/8.0 a^3*b/16.0 a*b^3/16.0;
0 0 0 0 0 0 0 0 0 0 1 0 a/2.0 b 0 a^2/4.0 a*b/2.0 3.0*b^2/4.0 a^3/8.0 3.0*a*b^2/8.0;
0 0 0 0 0 0 0 0 0 -1 0 -a -b/2.0 0 (-3.0)*a^2/4.0 -a*b/2.0 -b^2/4.0 0 (-3.0)*a^2*b/8.0 -b^3/8.0;
1 -a/2.0 b/2.0 -a*b/4.0 0 0 0 0 b*kxy-a*kx/2.0 -0.125*b^2*ky-a*b*kxy/2.0+0.125*a^2*kx b^2*kxy/4.0-a*b*kx/4.0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 -a/2.0 b/2.0 -a*b/4.0 b*ky/2.0-a*kxy a^2*kxy/4.0-a*b*ky/4.0 0.125*b^2*ky-a*b*kxy/2.0-0.125*a^2*kx 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 -a/2.0 b/2.0 a^2/4.0 -a*b/4.0 b^2/4.0 -a^3/8.0 a^2*b/8.0 -a*b^2/8.0 b^3/8.0 -a^3*b/16.0 -a*b^3/16.0;
0 0 0 0 0 0 0 0 0 0 1 0 -a/2.0 b 0 a^2/4.0 -a*b/2.0 3.0*b^2/4.0 -a^3/8.0 (-3.0)*a*b^2/8.0;
0 0 0 0 0 0 0 0 0 -1 0 a -b/2.0 0 (-3.0)*a^2/4.0 a*b/2.0 -b^2/4.0 0 (-3.0)*a^2*b/8.0 -b^3/8.0];
end%________________________________|
clear a b kx ky kxy kgx kgy kgxy lx ly mx my nx ny
%_______________________SYSTEM RIGID BODY MOTION CONTROL
for Element=1:No;
equation=size(A);
if equation(1)~=rank(A(:,:,Element))
display('This system connection matrix is badly scaled')
error('Control system shape function and geometricaly boundary conditions')
else
B=A(:,:,Element)^-1;
K(:,:,Element) =B'*Kd(:,:,Element)*B;
end%__________________________________________________|
clear equation Ku
end
%Shell stiffness matrix transform local axis to global axis
for Element=1:No;
Kg(:,:,Element)=Tson(:,:,Element)*K(:,:,Element)*Tson(:,:,Element)';
end
clear Element
%___________________________________SUPERPOSITION
%System stiffness matrix superposition
Ksis(Item,Item)=0;
for n=1:No;
for sat=1:24;
for sut=1:24;
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%_______________________|
clear sat sut n
%Ksis:System global stiffness matrix
%______________________________SYSTEM DISPLACEMENT
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 Ku
%_________________________________MODAL DISPLACEMENT
for s = 1 : No;
for m = 1 :24;
u = R(s, m);
if u ~=0; Hu(s,m,:) = D(u) ;
else Hu(s,m,:) = 0 ;
end%
end%
end%__________________________________________________|
%Hu:per element node`s displacements
clear s m u
%Elements global and local node reactions
for v=1 :No;
Pg(:,v) = Kg(:,:,v)*Hu(v,:)';
Pl(:,v) = Tson(:,:,v)'*Pg(:,v);
end
clear Node Nom No v
%display('GLOBAL SYSTEM NODE DISPLACEMENTS (metre)');
%for s=1:size(D,1);
%fprintf(' freedom =% .f % .7f \n',s,D(s));
%end
%fprintf(' \n \n \n')
display('ELEMENT DISPLACEMENTS (metre)');
for s=1:size(Hu',1);
fprintf(' [% .7f] [% .7f] [% .7f] [% .7f] \n',Hu(:,s)');
if mod(s,6)==0; fprintf(' \n');end;
end
fprintf(' \n \n \n')
display('LOCAL AXIS ELEMENTS NODE REACTIONS (SI)Kn');
for s=1:size(Pl,1);
fprintf(' [% .7f] [% .7f] [% .7f] [% .7f] \n',Pl(s,:));
if mod(s,5)==0; fprintf(' \n');end;
end
fprintf(' \n \n \n')
display('GLOBAL AXIS ELEMENTS NODE REACTIONS (SI)Kn');
for s=1:size(Pg,1);
fprintf(' [% .7f] [% .7f] [% .7f] [% .7f] \n',Pg(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,6)==0; fprintf(' \n');end;
end
clear s ans i j no alfax alfay dim
%=========================ANALYSIS RESULTS(Global reactions)
% [1] [2] [3] [4]
%[ 8.2120] [ 6.4526] [-0.0000] [ 1.9166] fx
%[ 8.2120] [ 0.0000] [ 6.4526] [ 1.9166] fy
%[ 4.9024] [-2.5000] [-2.5000] [-0.0000] fz
%[ 0.0533] [-0.0000] [-0.3543] [ 0.0178] Mxx
%[-0.0533] [ 0.3543] [-0.0000] [-0.0178] Myy
%[ 0.0000]1 [ 0.1084]2 [-0.1084]4 [-0.0000]5 Mzz
%[-6.4526] [-8.2120] [-1.9166] [-0.0000]
%[-0.0000] [ 8.2120] [ 1.9166] [ 6.4526]
%[-2.5000] [ 4.9024] [ 0.0000] [-2.5000]
%[ 0.0000] [ 0.0533] [ 0.0178] [-0.3543]
%[-0.3543] [ 0.0533] [ 0.0178] [-0.0000]
%[-0.1084]2 [ 0.0000]3 [-0.0000]5 [ 0.1084]6
%[-1.9166] [-0.0000] [-6.4526] [-8.2120]
%[-1.9166] [-6.4526] [ 0.0000] [-8.2120]
%[-0.0000] [-2.5000] [-2.5000] [ 4.9024]
%[-0.0178] [ 0.3543] [-0.0000] [-0.0533]
%[ 0.0178] [ 0.0000] [-0.3543] [ 0.0533]
%[ 0.0000]5 [-0.1084]6 [ 0.1084]8 [-0.0000]9
%[ 0.0000] [ 1.9166] [ 8.2120] [ 6.4526]
%[-6.4526] [-1.9166] [-8.2120] [-0.0000]
%[-2.5000] [-0.0000] [ 4.9024] [-2.5000]
%[ 0.3543] [-0.0178] [-0.0533] [ 0.0000]
%[ 0.0000] [-0.0178] [-0.0533] [ 0.3543]
%[ 0.1084]4 [ 0.0000]5 [-0.0000]7 [-0.1084]8
%____________Node (1)
% [1]
%[ 8.2120]
%[ 8.2120]
%[ 4.9024] = SUPPORT
%[ 0.0533]
%[-0.0533]
%[ 0.0000]1
%____________Node (2) [K][D]+{f}=[P]
% [1] [2] Out Force
%[-6.4526] [ 6.4526] [0.0000]
%[-0.0000] [ 0.0000] [0.0000]
%[-2.5000] + [-2.5000] + [5.0000] = zero
%[ 0.0000] [-0.0000] [0.0000]
%[-0.3543] [ 0.3543] [0.0000]
%[-0.1084]2 [ 0.1084]2 [0.0000]
%____________Node (3)
% [3]
%[-8.2120]
%[ 8.2120]
%[ 4.9024] =SUPPORT
%[ 0.0533]
%[ 0.0533]
%[ 0.0000]3
%
%____________Node (4)
% [3] out load
%[-0.0000] [ 0.0000]
%[ 6.4526] [-6.4526]
%[-2.5000] + [-2.5000] +5.0000 = zero
%[-0.3543] [ 0.3543]
%[-0.0000] [ 0.0000]
%[-0.1084]4 [ 0.1084]
