16.6 Computer Program for Finite Element Free Vibration Analysis of Plates
597
Table 16.9 (continued)
write(8,*) '
the overall stiffness matrix'
call crunch (nt,nbh,ntot,gstif,kind)
c****************************************************************************
c
c
fitting of boundary conditions
c
c****************************************************************************
nd=0
nx21=nfree*(nx+nx)+1
call bound (nbon1,1,nx21,nfree,gstif,kind,nd)
nd=1
call bound (nbon2,nx21,nt,nxf,gstif,kind,nd)
nin=nxf-nfree+1
nst = nt-nx21-nfree
call bound (nbon2,nin,nst,nxf,gstif,kind,nd)
nd=0
nst = nst+2
call bound (nbon3,nst,nt,nfree,gstif,kind,nd)
nd=1
call bound (nbon4,1,nst,nxf,gstif,kind,nd)
nx21 = nx21+ nfree
nst=nt-nxf+1
call bound (nbon4,nx21,nst,nxf,gstif,kind,nd)
call r8usiv(nt,nrqd,gstif,kind,amaso,mind,bd,nfree,W)
stop
end
c************************************************************************
c
c
main program ends
c
c*************************************************************************
c*************************************************************************
c
SHAPE SUBROUTINE IS GIVEN BELOW
c
c**************************************************************************
subroutine shape (zi,eta,zx,zy,sfd,sfdr,sfdz,sfdzr,sfde,sfder,
* bb,zac,xg,yg,nn,nenp,nfree,elmas,rho)
dimension sfd(8),sfde(8),sfdz(8),sfder(8),sfdzr(8),sfdr(8),
* bb(5,24),zx(8),zy(8),xg(8),yg(8),elmas(24)
do 1 i=1,nn
1
elmas(i)=0.0
sum=0.0
do 517 k=1,7,2
ziz= zi*zx(k)
etaz= eta*zy(k)
sfd(k)= (1.0+ziz)*(1.0+etaz)*(ziz+etaz-1.0)/4.0
sfdz(k)= (1.0+etaz)*zx(k)*(2.0*ziz+etaz)/4.0
517
sfde(k)= (1.0+ziz)*zy(k)*(2.0*etaz+ziz)/4.0
do 518 ll=2,6,4
etaz=eta*zy(ll)
sfd(ll)= (1.0-zi*zi)*(1.0+etaz)/2.0
sfdz(ll)= -zi*(1.0+etaz)
518
sfde(ll)= (1.0-zi*zi)*zy(ll)/2.0
do 519 kk=4,8,4
ziz= zi*zx(kk)
sfdz(kk)= zx(kk)*(1.0-eta*eta)/2.0
sfde(kk)= -eta*(1.0+ziz)
519
sfd(kk)= (1.0+ziz)*(1.0-eta*eta)/2.0
do 520 ii=1,5
do 520 jj=1,nenp
520
bb(ii,jj)=0.0
(continued)