Civil Engineering Reference
In-Depth Information
√
0
.
004931
=
0
.
070 as compared with 0
.
068 previously computed for a similar problem
with lumped mass.
Program 10.4 Eigenvalue analysis of an elastic solid in plane strain using 4-node
rectangular quadrilaterals. Lanczos Method. Lumped mass. Element-by-element for-
mulation. Mesh numbered in
x
-or
y
-direction.
PROGRAM p104
!-------------------------------------------------------------------------
! Program 10.4: Eigenvalue analysis of an elastic solid in plane strain
! using 4-node rectangular quadrilaterals. Lanczos Method.
! Lumped mass. Element by element formulation.
! Mesh numbered in x- or y-direction.
!-------------------------------------------------------------------------
USE main; USE geom; iMPLICIT NONE
INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
INTEGER::i,iel,iflag=-1,iters,jflag,k,lalfa,leig,lp=6,lx,lz,ndim=2,ndof, &
neig=0,nels,neq,nip=4,nmodes,nn,nod=4,nodof=2,nprops=3,np_types,nr,
&
nst=3,nxe,nye
REAL(iwp)::acc,det,el,er,one=1.0_iwp,zero=0.0_iwp
CHARACTER(LEN=15)::element='quadrilateral'
!--------------------------- dynamic arrays-------------------------------
INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),g_num(:,:),jeig(:,:),nf(:,:),&
nu(:),num(:)
REAL(iwp),ALLOCATABLE::alfa(:),bee(:,:),beta(:),coord(:,:),dee(:,:), &
del(:),der(:,:),deriv(:,:),diag(:),ecm(:,:),eig(:),fun(:),g_coord(:,:),&
jac(:,:),km(:,:),mm(:,:),points(:,:),prop(:,:),ua(:),udiag(:),va(:),
&
vdiag(:),v_store(:,:),weights(:),w1(:),x(:),x_coords(:),y(:,:),
&
y_coords(:),z(:,:)
!----------------------input and initialisation---------------------------
OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
READ(10,*)nxe,nye,nod,np_types
CALL mesh_size(element,nod,nels,nn,nxe,nye); ndof=nod*nodof
ALLOCATE(nf(nodof,nn),points(nip,ndim),dee(nst,nst),g_coord(ndim,nn), &
coord(nod,ndim),fun(nod),jac(ndim,ndim),weights(nip),g_num(nod,nels), &
der(ndim,nod),deriv(ndim,nod),bee(nst,ndof),num(nod),km(ndof,ndof),
&
g(ndof),g_g(ndof,nels),mm(ndof,ndof),ecm(ndof,ndof),eig(leig),x(lx),
&
del(lx),nu(lx),jeig(2,leig),alfa(lalfa),beta(lalfa),z(lz,leig),
&
prop(nprops,np_types),x_coords(nxe+1),y_coords(nye+1),etype(nels))
READ(10,*)prop; etype=1; IF(np_types>1)READ(10,*)etype
READ(10,*)x_coords,y_coords
nf=1; READ(10,*)nr,(k,nf(:,k),i=1,nr); CALL formnf(nf); neq=MAXVAL(nf)
READ(10,*)nmodes,el,er,lalfa,leig,lx,lz,acc
!------------------loop the elements to set up global arrays--------------
elements_1: DO iel=1,nels
CALL geom_rect(element,iel,x_coords,y_coords,coord,num,'y')
CALL num_to_g(num,nf,g); g_num(:,iel)=num
g_coord(:,num)=TRANSPOSE(coord); g_g(:,iel)=g
END do elements_1
WRITE(11,'(A,I5,A)')" There are",neq," equations"
ALLOCATE(ua(0:neq),va(0:neq),vdiag(0:neq),diag(0:neq),udiag(0:neq),
&
v_store(0:neq,lalfa),w1(0:neq),y(0:neq,leig))
ua=zero; va=zero; eig=zero; jeig=0; x=zero; del=zero; nu=0; alfa=zero
beta=zero; diag=zero; udiag=zero; w1=zero; y=zero; z=zero
CALL sample(element,points,weights)
!--------------- element stiffness integration and assembly---------------
elements_2: DO iel=1,nels
