DAFTAR PUSTAKA 1.
Gallagher, Richard H.(1975), Finite Element Analysis Fundamental, PrenticeHall, Englewood Clifs, N J.
2.
Gallagher, Richard H., Ziemian, Ronald D. and McGuire, William (2000), Matrix Structural Analysis, John Wiley & Sons, New York.
3.
Chen, W.F. and Lui, E.M. (1987), Structural Stability, Elsevier, New York.
4.
Hariandja, Binsar H. (1997), Analisis Struktur Berbentuk Rangka Dalam Formulasi Matriks, Aksara Hutasada, Bandung, Jawa Barat.
5.
Kassimali, Aslam (1999), Matrix Analysis of Structures, International Thomson Publishing, USA.
6.
Simitses, George J. (1976), An Introduction to the Elastic Stability of Structures, Prentice-Hall, Englewood Clifs, N J.
7.
Chajes, Alexander (1974), Principles of Structural Stability Theory, PrenticeHall, Englewood Clifs, N J.
8.
Kirby, P.A and Nethercot, D.A. (1979), Design for Structural Stability, Collins, London.
9.
Tauchert, Theodore R.(1974), Energy Principles in Structural Mechanics, McGraw-Hill.
10. Segerlind, Larry J. (1984), Applied Finite Element Analysis, John Wiley & Sons, USA 11. Zienkiewicz, O.C. (1989), The Finite Element Method, McGraw-Hill. 12. Wong, Y.W and Bang, Hyochoong (2000), The Finite Element Method Using MATLAB Second edition, CRC Press, USA. 13. Alkaff, M.Firdaus (2004), Matlab 6 untuk Teknik Sipil, Maxikom, Palembang. 14. Cook, Robert D. and Young, Warren C. (1985), Advanced Mechanics of Materials, Macmilan Publishing Company, New York, USA. 15. Palm, William J. (2001), Introduction to Matlab 6 for Engineers, Mc Graw Hill, New York, USA. 16. Cook, Robert D. (1981), Konsep dan Aplikasi Metode Elemen Hingga, P.T. Eresco, Bandung, Indonesia.
DP
17. Cook, Robert D., Malkus, David S., and Plesha, Michael E. (1989), Concepts and Applications of Finite Element Analysis, Third edition, John Wiley and Sons, USA.
DP
A-1 LAMPIRAN A. BLOK DIAGRAM PROGRAM PENENTUAN BEBAN KRITIS ELASTIK. Input Koordinat Member dan Properti elemen
restraint
Cases
Proses Perhitungan awal Jumlah elemen Nodes tiap elemen Jumlah d.o.f per node Total node pada system Jumlah total d.o.f pada system Tempat untuk matrik kekakuan lentur struktur Tempat untuk matrik geometri struktur Tempat untuk matrik index
[ ]
Matriks Transformasi T
Matriks kekakuan elemen dalam koordinat lokal
[ S ]g [ S ]l Matriks kekakuan elemen dalam koordinat global
[ k ]l [ k ]g [ ]l
Matriks kekakuan struktur K
K =0
Output : Pcritis
[ ]g
& K
B-1 LAMPIRAN B. BLOK DIAGRAM PROGRAM GAYA DALAM Input
Forces diambil dari nilai Pcritis
Proses Perhitungan awal Jumlah elemen Nodes tiap elemen Jumlah d.o.f per node Total node pada system Jumlah total d.o.f pada system Tempat untuk matrik kekakuan lentur struktur Tempat untuk matrik index
[ ]
Matriks Transformasi T
Matriks kekakuan elemen dalam koordinat lokal
[ S ]l Matriks kekakuan elemen dalam koordinat global
[ k ]l [ ]l
Matriks kekakuan struktur K Deformasi Global dikekang = delta
d.o.f.
yang
tidak
Deformasi Global semua d.o.f = DELTAR Invers Transformasi Deformasi local per elemen =deltlocal
Output : Gaya Dalam Penampilan gaya dalam semua elemen portal bidang
C1-1 LAMPIRAN C1. FLOWCHART SUBROUTINE DARI PROGRAM BEBAN KRITIS ELASTIS INDEXS(iel,mem,ndof)
str1=(mem(iel,2)-1)*ndof
str2=(mem(iel,3)-1)*ndof
index= [str1+1 str1+2 str2+1 str2+2 ]
kembali
feasmbl1SS(kk,k,index)
edof=length(index)
i=1:edof
ii=index(i)
j=1:edof
jj=index(j) kk(ii,jj)=kk(ii,jj)+k(i,j)
kembali
C1-2
feframeS(el,xi,leng,a1,rho,beta,ipt,G,ident)
a=el*a1/leng
c=el*xi/(leng^3)
[ S ]l [T ]
[ k ]l = [T ] [ S ]l [T ] T
gb=6/(5*leng) gc=1/10 gd=2*1*leng/15 ge=1*leng/30
[ S ]g [ k ]g = [T ] [ S ]g [T ] T
Kembali
C1-3 feasmbl1s(kk,k,index)
edof=length(index)
i=1:edof
ii=index(i)
j=1:edof
jj=index(j) kk(ii,jj)=kk(ii,jj)+k(i,j)
kembali
C1-4 Bcdofs(res,ndof)
n=size(res) nn=n(1)*(n(2)-2) bcdof=zeros(1,nn)
i=1:n(1)
resi=res(i,:) c=resi(1) ident=[ resi(3) resi(4) ];
ident==[ 1 1]
ya
bcdof(2*i-1)=ndof*c-1 bcdof(2*i)=ndof*c
ident==[ 1 0]
ya
bcdof(2*i-1)=ndof*c-1 bcdof(2*i)=0 bcdof(2*i-1)=ndof*c-1 bcdof(2*i)=0
Kembali
C1-5 RESTRAINTS(kk,kkG,bcdof)
nnn=length(bcdof)
i=1:nnn
elbc=bcdof(i)
elbc==0
YA
kk=kk kkG=kkG kk(elbc-(i-1),:)=[] kk(:,elbc-(i-1))=[] kkG(elbc-(i-1),:)=[] kkG(:,elbc-(i-1))=[]
Kembali
C2-1 LAMPIRAN C2. FLOWCHART SUBROUTINE DARI PROGRAM GAYA DALAM INDEXSGD(iel,mem,ndof)
str1=(mem(iel,2)-1)*ndof
str2=(mem(iel,3)-1)*ndof
index= [str1+1 str1+2 str1+3 str2+1 str2+2 str2+3]
kembali
feframe2SS(el,xi,leng,a1,beta)
a=el*a1/leng c=el*xi/(leng^3)
[ S ]l [T ]
[ k ]l = [T ] [ S ]l [T ] T
kembali
C2-2 feasmbl1G(kkG,kG,index)
edof=length(index)
i=1:edof
ii=index(i)
j=1:edof
jj=index(j) kkG(ii,jj)=kkG(ii,jj)+kG(i,j)
kembali
C2-3 Bcdof_gayadalam(res,ndof)
n=size(res) nn=n(1)*(n(2)-2) bcdof=zeros(1,nn)
i=1:n(1)
resi=res(i,:) c=resi(1) ident=[ resi(2) resi(3) resi(4) ]
ident==[ 1 11]
ya bcdof(3*i-2)=ndof*c-2; bcdof(3*i-1)=ndof*c-1; bcdof(3*i)=ndof*c;
ident==[ 1 1 0]
ya bcdof(3*i-2)=ndof*c-2 bcdof(3*i-1)=ndof*c-1
bcdof(3*i)=0 bcdof(3*i-2)=ndof*c-2 bcdof(3*i-1)=ndof*c-1 bcdof(3*i)=0
Kembali
C2-4 RESTRAINT_gayadalam(kk,kkG,bcdof)
nnn=length(bcdof)
i=1:nnn
elbc=bcdof(i)
elbc==0
YA
kk=kk ff=ff kk(elbc-(i-1),:)=[] kk(:,elbc-(i-1))=[] ff(elbc-(i-1))=[]
Kembali
C2-5
feframe_gayadalam(el,xi,leng,a1,beta)
a=el*a1/leng c=el*xi/(leng^3)
[ S ]l [T ]
[ k ]l = [T ] [ S ]l [T ] T
kembali
DELT(DELTAR,ndof,node1,node2)
i=1:ndof
delt1(i)=DELTAR(node1*3-(3-i))
j=1:ndof
delt2(j)=DELTAR(node2*3-(3-j))
delt = [delt1 delt2]
kembali
C-1 LAMPIRAN C. FLOWCHART PROGRAM BEBAN KRITIS ELASTIS MULAI
Baca: Koordinat (Coord) Properti elemen dan member (mem) Restraint (res) Cases (cs)
Jumlah elemen (nel)
Node per elemen (nnel)
Jumlah d.o.f. pertitik (ndof)
Total node pada struktur (nnode)
Total d.o.f dalam struktur (sdof)
Matriks nol kekakuan struktur lentur (kk)
Matriks nol kekakuan struktur geometri (kkG)
Matriks nol gaya (ff)
Matriks nol indexs
iel=1:nel
INDEXS(iel,mem,ndof)
A1
B1
C-2 A1
node1=mem(iel,2 node2=mem(iel,3) x1=x(node1); y1=y(node1) x2=x(node2); y2=y(node2) leng=sqrt((x2-x1)^2+(y2-y1)^2) el=mem(iel,4) G=mem(iel,5) a1=mem(iel,6) xi=mem(iel,7) rho=mem(iel,8) ipt=mem(iel,9)
(y2-y1)>0
ya
Beta = (y2-y1)<0
ya
Beta = Beta = a tan ( ( y2 − y1 ) ( x2 − x1 ) )
ident=cs(iel,2)
feframeS(el,xi,leng,a1,rho,beta,ipt,G,ident)
feasmbl1G(kkG,kG,index)
feasmbl1s(kk,k,index)
B1
C1
3 pi 2
pi 2
C-3
C1
Bcdofs(res,ndof)
RESTRAINTS(kk,kkG,bcdof)
Critis=eig(kk,kkG)
G
C-4
G
Baca : Forces (forces) dari nilai Pcritis
Jumlah elemen (nel)
Node per elemen (nnel)
Jumlah d.o.f per titik (ndof)
Total node pada struktur (nnode)
Total d.o.f. struktur (sdof)
Matrik nol kekakuan struktur lentur (kk)
Matrik nol index (index)
iel=1:nel
INDEXSGD (iel,mem,ndof)
Node 1 Node 2
x1 , y1 x2 , y2
A
B
C-5 A
leng, el, G, a1, xi, rho, ipt
( y2 − y1 ) > 0
ya
( y2 − y1 ) < 0
Beta = Beta = a tan ( ( y2 − y1 ) ( x2 − x1 ) )
feframe (el,xi,leng,a1,beta)
feasmbl1SS(kk,k,index
B
Bcdof gayadalam(res,ndof)
RESTRAINT gayadalam(kk,ff,bcdof)
delta
C
Beta =
ya
3 pi 2
pi 2
C-6 C
Matriks nol DELTAR
i=1:length(bcdof)
id = bcdof(i)
id=0
ya DELTAR=DELTAR
DELTAR(id)=id
zero=0
j=1:sdof
idd=DELTAR(j)
idd ==0
DELTAR(j)=0 zero=zero+1
edof=nnel*ndof SGD=zeros(edof,nel)
D
ya
DELTAR(j)=delta(j-zero)
C-7 D
iel=1:nel
iel
node1=mem(iel,2) node2=mem(iel,3) x1=x(node1); y1=y(node1) x2=x(node2); y2=y(node2) leng=sqrt((x2-x1)^2+(y2-y1)^2) el=mem(iel,4) G=mem(iel,5 a1=mem(iel,6) xi=mem(iel,7) rho=mem(iel,8) ipt=mem(iel,9);
(y2-y1)>0
ya
Beta = (y2-y1)<0
ya
Beta = Beta = a tan ( ( y2 − y1 ) ( x2 − x1 ) )
feframe gayadalam(el,xi,leng,a1,beta)
rinv=inv(r)
E
3 pi 2
pi 2
C-8 E
DELT(DELTAR,ndof,node1,node2)
deltlocal=rinv*delt'
GD=kl*deltlocal
SGD(:,iel)=GD
SELESAI
D-1 LAMPIRAN D. PROGRAM BEBAN KRITIS ELASTIS DAN GAYA DALAM DARI BEBAN KRITIS ELASTIS %============================================================% % PROGRAM FOR CALCULATING ELASTIC CRITICAL LOAD AND % % INTERNAL FORCE % % BY % % FRANSISCA MARIA FARIDA (25004045) % % BANDUNG INSTITUTE OF TECHNOLOGY % % 2007 % %============================================================% % PROGRAM FOR CALCULATING ELASTIC CRITICAL LOAD % %============================================================% clear all clc format short; %============================================================% % Inputting the magnitude of elastic critical load %============================================================% load cases.txt; cs=cases; % Unit of measurement : Kgf,m % ============================================================% % Inputing coordinates using notepad % ============================================================% load coordinates.txt; coord=coordinates; j=coord(:,1);%j=joint x=coord(:,2);%x=x coordinate of joint y=coord(:,3);%y=y coordinate of joint %============================================================% % Inputing members and its properties using notepad %============================================================% load members1.txt; mem=members1; %============================================================% % Inputing restraints using notepad %============================================================% % 0=release % 1=restrainted load restraints.txt; res=restraints; jr=res(:,1);%jr=joint xr=res(:,2);%xr=restraint at x direction yr=res(:,3);%yr=restraint at y direction teta_r=res(:,4);%teta_r=restraint at teta direction
D-2 %============================================================% % Initial Calculation %============================================================% jm=length(mem(:,1)); jc=length(coord(:,1)); nel=jm ; % nel=number of element nnel=2; % nnel=nodes per element ndof=2; % ndof=number of degree of freedom pernode nnode=jc; % nnode=total number of nodes in system sdof=nnode*ndof; % sdof=sum of degree of freedom kk=zeros(sdof,sdof);% kk=flexure structural stiffness matrix kkG=zeros(sdof,sdof); % mm=mass structural stiffness matrix ff=zeros(sdof,1); index=zeros(nnel*ndof,1);% index for matrix assemblying %============================================================% % Element stiffness matrix and Assembling Structural stiffness matrix %============================================================% for iel=1:nel index=INDEXS(iel,mem,ndof); % index for assembly node1=mem(iel,2); % node1 per element node2=mem(iel,3); % node2 per element x1=x(node1); y1=y(node1); % coordinate of node1 per element x2=x(node2); y2=y(node2); % coordinate of node2 per element leng=sqrt((x2-x1)^2+(y2-y1)^2); % length per element el=mem(iel,4); % el=elasticity modulus G=mem(iel,5); % G=shear modulus a1=mem(iel,6); % a1=area of section xi=mem(iel,7); % xi=inertia of section rho=mem(iel,8); % rho=mass density ipt=mem(iel,9); % mass matrix choice if(y2-y1)>0; % angle of each element beta=pi/2; elseif (y2-y1)<0 beta=3*pi/2; else beta=atan((y2-y1)/(x2-x1)); end ident=cs(iel,2); [k,kG]=feframeS(el,xi,leng,a1,rho,beta,ipt,G,ident) %k=global flexure stiffness matrix %and kG=global geometry stiffness matrix kkG=feasmbl1G(kkG,kG,index);%kkG=geometry structural stiffness matrix kk=feasmbl1s(kk,k,index)% kk=flexure structural stiffness matrix end [bcdof]=Bcdofs(res,ndof) ; % vector containing restraint dof [kk,kkG]=RESTRAINTS(kk,kkG,bcdof); % Restraint
D-3 Critis=eig(kk,kkG) %============================================================% % PROGRAM FOR CALCULATING INTERNAL FORCE % % WITH % % ELASTIC CRITICAL LOAD AS EXTERNAL FORCE % %============================================================% % Inputting force using notepad %============================================================% load forces.txt; ffGD=forces; ff=forces; %============================================================ % % calculation %============================================================ % jm=length(mem(:,1)); jc=length(coord(:,1)); nel=jm ; % nel=number of element nnel=2; % nnel=nodes per element ndof=3; % ndof=number of degree of freedom pernode nnode=jc; % nnode=total number of nodes in system sdof=nnode*ndof; % sdof=sum of degree of freedom kk=zeros(sdof,sdof);% kk=flexure structural stiffness matrix index=zeros(nnel*ndof,1);% index for matrix assemblying for iel=1:nel %iel=increment index=INDEXSGD(iel,mem,ndof); % index for assembly node1=mem(iel,2); % node1 per element node2=mem(iel,3); % node2 per element x1=x(node1); y1=y(node1); % coordinate of node1 per element x2=x(node2); y2=y(node2); % coordinate of node2 per element leng=sqrt((x2-x1)^2+(y2-y1)^2); % length per element el=mem(iel,4); % el=elasticity modulus G=mem(iel,5); % G=shear modulus a1=mem(iel,6); % a1=area of section xi=mem(iel,7); % xi=inertia of section rho=mem(iel,8); % rho=mass density ipt=mem(iel,9); % mass matrix choice if(y2-y1)>0; % angle of each element beta=pi/2; elseif (y2-y1)<0 beta=3*pi/2; else beta=atan((y2-y1)/(x2-x1)); end k=feframe2SS(el,xi,leng,a1,beta) %k=global flexure stiffness matrix
D-4 kk=feasmbl1SS(kk,k,index) end
% kk=flexure structural stiffness matrix
[bcdof]=Bcdof_gayadalam(res,ndof)%bcdof=Boundary condition degrees of freedom [kk,ff]=RESTRAINT_gayadalam(kk,ff,bcdof)%RESTRAINT=flexure and geometry structural stiffness matrix without restraint delta=inv(kk)*ff;%delta=displacement DELTAR =zeros(sdof,1) for i=1:length(bcdof) id = bcdof(i) if id==0 DELTAR=DELTAR else DELTAR(id)=id end end
%Initialization of displacement vector %length of bcdof vector (bcdof : vector %containing restrained dof)
zero=0 for j=1:sdof idd=DELTAR(j) if idd ==0 DELTAR(j)=delta(j-zero) else DELTAR(j)=0 zero=zero+1 end end edof=nnel*ndof; SGD=zeros(edof,nel); for iel=1:nel iel node1=mem(iel,2); node2=mem(iel,3); x1=x(node1); y1=y(node1); x2=x(node2); y2=y(node2); leng=sqrt((x2-x1)^2+(y2-y1)^2); el=mem(iel,4); G=mem(iel,5); a1=mem(iel,6); xi=mem(iel,7); rho=mem(iel,8); ipt=mem(iel,9); if(y2-y1)>0;
%iel=increment % node1 per element % node2 per element % coordinate of node1 per element % coordinate of node2 per element % length per element % el=elasticity modulus % G=shear modulus % a1=area of section % xi=inertia of section % rho=mass density
D-5 beta=pi/2; % beta=transformation angle elseif (y2-y1)<0 beta=3*pi/2; % beta=transformation angle else beta=atan((y2-y1)/(x2-x1)); % beta=transformation angle end [kl,r]=feframe_gayadalam(el,xi,leng,a1,beta) ;%k=local flexure stiffness matrix rinv=inv(r) delt=DELT(DELTAR,ndof,node1,node2);%delt=deformasi global per elemen deltlocal=rinv*delt'%deltlocal=deformasi lokal per elemen GD=kl*deltlocal; %GD=internal forces per element SGD(:,iel)=GD end
D1-1 LAMPIRAN D1. SUBROUTINE BEBAN KRITIS ELASTIK DAN GAYA DALAM SUBROUTINE INDEXS function[index]=INDEXS(iel,mem,ndof) str1=(mem(iel,2)-1)*ndof str2=(mem(iel,3)-1)*ndof index= [str1+1 str1+2 str2+1 str2+2 ] SUBROUTINE feframeS function [k,kG]=feframeS(el,xi,leng,a1,rho,beta,ipt,G,ident) a=el*a1/leng; c=el*xi/(leng^3); kl = [ 12*c -6*leng*c -12*c -6*leng*c ;... -6*leng*c 4*leng^2*c 6*leng*c 2*leng^2*c;... -12*c 6*leng*c 12*c 6*leng*c ;... -6*leng*c 2*leng^2*c 6*leng*c 4*leng^2*c]; r=[ cos(beta) 0 0 0;... 0 1 0 0;... 0 0 cos(beta) 0;... 0 0 0 1]; k=r'*kl*r; % global stiffness matrix %ga=(30*1)/(30*leng); gb=6/(5*leng); gc=1/10; gd=2*1*leng/15; ge=1*leng/30; %if beta==0 % kG=zeros(4,4); %else klG = ident*[ gb -gc -gb -gc ;... -gc gd gc -ge ;... -gb gc gb gc ;... -gc -ge gc gd]; %klG=local geometry stiffness matrix %r=[ cos(beta) 0 0 0;... % 0 1 0 0;... % 0 0 cos(beta) 0;... % 0 0 0 1]; %r=transformation matrix kG=r'*klG*r; % kG=global geometry stiffness matrix
D1-2 SUBROUTINE feasmbl1G function [kkG]=feasmbl1G(kkG,kG,index) edof=length(index); for i=1:edof; ii=index(i); for j=1:edof; jj=index(j); kkG(ii,jj)=kkG(ii,jj)+kG(i,j) end end SUBROUTINE feasmbl1s function [kk]=feasmbl1s(kk,k,index) edof=length(index); for i=1:edof; ii=index(i); for j=1:edof; jj=index(j); kk(ii,jj)=kk(ii,jj)+k(i,j) end end SUBROUTINE Bcdofs function [bcdof]=Bcdofs(res,ndof) n=size(res); nn=n(1)*(n(2)-2); bcdof=zeros(1,nn); for i=1:n(1) resi=res(i,:) c=resi(1) ident=[ resi(3) resi(4) ]; if ident==[ 1 1] bcdof(2*i-1)=ndof*c-1; bcdof(2*i)=ndof*c; elseif ident==[ 1 0] bcdof(2*i-1)=ndof*c-1; bcdof(2*i)=0; else bcdof(2*i-1)=ndof*c-1; bcdof(2*i)=0; end end
D1-3 SUBROUTINE RESTRAINTS function [kk,kkG]=RESTRAINTS(kk,kkG,bcdof) nnn=length(bcdof); for i=1:nnn elbc=bcdof(i) if elbc==0 kk=kk kkG=kkG else kk(elbc-(i-1),:)=[] kk(:,elbc-(i-1))=[] kkG(elbc-(i-1),:)=[] kkG(:,elbc-(i-1))=[] end end SUBROUTINE INDEXSGD function[index]=INDEXSGD(iel,mem,ndof) str1=(mem(iel,2)-1)*ndof str2=(mem(iel,3)-1)*ndof index= [str1+1 str1+2 str1+3 str2+1 str2+2 str2+3] SUBROUTINE feframe2SS function [k]=feframe2SS(el,xi,leng,a1,beta) a=el*a1/leng;%a=(EA)/L c=el*xi/(leng^3);%c=(EI)/(L^3) kl = [ a 0 0 -a 0 0 ;... 0 12*c 6*leng*c 0 -12*c 6*leng*c ;... 0 6*leng*c 4*leng^2*c 0 -6*leng*c 2*leng^2*c;... -a 0 0 a 0 0 ;... 0 -12*c -6*leng*c 0 12*c -6*leng*c ;... 0 6*leng*c 2*leng^2*c 0 -6*leng*c 4*leng^2*c]; %kl=local flexure stiffness matrix r=[ cos(beta) sin(beta) 0 0 0 0;... -sin(beta) cos(beta) 0 0 0 0;... 0 0 1 0 0 0;... 0 0 0 cos(beta) sin(beta) 0;... 0 0 0 -sin(beta) cos(beta) 0;... 0 0 0 0 0 1]; %r=transformation matrix k=r'*kl*r; % k=global flexure stiffness matrix SUBROUTINE feasmbl1SS function [kk]=feasmbl1SS(kk,k,index) edof=length(index);
D1-4 for i=1:edof; ii=index(i); for j=1:edof; jj=index(j); kk(ii,jj)=kk(ii,jj)+k(i,j); end end SUBROUTINE Bcdof_gayadalam function [bcdof]=Bcdof_gayadalam(res,ndof) n=size(res); nn=n(1)*(n(2)-2); bcdof=zeros(1,nn); for i=1:n(1) resi=res(i,:) c=resi(1) ident=[ resi(2) resi(3) resi(4) ]; if ident==[ 1 1 1] bcdof(3*i-2)=ndof*c-2; bcdof(3*i-1)=ndof*c-1; bcdof(3*i)=ndof*c; elseif ident==[ 1 1 0] bcdof(3*i-2)=ndof*c-2; bcdof(3*i-1)=ndof*c-1; bcdof(3*i)=0; else bcdof(3*i-2)=ndof*c-2; bcdof(3*i-1)=ndof*c-1; bcdof(3*i)=0; end end SUBROUTINE RESTRAINT_gayadalam function [kk,ff]=RESTRAINT_gayadalam(kk,ff,bcdof) nnn=length(bcdof); for i=1:nnn elbc=bcdof(i) if elbc==0 kk=kk ff=ff else kk(elbc-(i-1),:)=[] kk(:,elbc-(i-1))=[] ff(elbc-(i-1))=[] end
D1-5 end SUBROUTINE feframe_gayadalam function [kl,r]=feframe_gayadalam(el,xi,leng,a1,beta) %[k,m] a=el*a1/leng; c=el*xi/(leng^3); kl = [ a 0 0 -a 0 0 ;... 0 12*c 6*leng*c 0 -12*c 6*leng*c ;... 0 6*leng*c 4*leng^2*c 0 -6*leng*c 2*leng^2*c;... -a 0 0 a 0 0 ;... 0 -12*c -6*leng*c 0 12*c -6*leng*c ;... 0 6*leng*c 2*leng^2*c 0 -6*leng*c 4*leng^2*c]; r=[ cos(-beta) sin(-beta) 0 0 0 0 ;... -sin(-beta) cos(-beta) 0 0 0 0 ;... 0 0 1 0 0 0 ;... 0 0 0 cos(-beta) sin(-beta) 0 ;... 0 0 0 -sin(-beta) cos(-beta) 0 ;... 0 0 0 0 0 1]; SUBROUTINE DELT function delt=DELT(DELTAR,ndof,node1,node2) for i=1:ndof delt1(i)=DELTAR(node1*3-(3-i)); end % for j=1:ndof delt2(j)=DELTAR(node2*3-(3-j)); end delt = [delt1 delt2];
