modes.for

一个组合程序

******************************
*       Program to generate the response of an N-DOF structure
*       to any arbitrary applied loading using the 
*       MODE SUPERPOSITION METHOD
******************************
*       Program Author : Sanjoy Chakraborty
*                        Auburn University
******************************

	program mode_s
	use msimsl
	implicit real *8 (a-h,o-z)
	real *8 mr,kr
	character *12 flin,fout
	dimension mr(10,1),cr(10,1),kr(10,1),phi(10,10),tt(20,10),
     1  ft(20,10),fr(6000,10),q(6000,10),q1(6000,10),q2(6000,10),
     2  qt(10,1),q1t(10,1),q2t(10,1),nfor(10),xmin(10,1),vmin(10,1)
	data ndim/10/


*******************************
*       Assign I/O Devices and read in system parameters
*******************************

	write(*,50)
50      format(///,' Program for Eigensolution and Mode Superposition'
     1  ,//,' Program Author : Sanjoy Chakraborty, Auburn University'
     2  ,//////,' Enter Input File Name : ')
	read(*,'(a)')flin
	open(1,file=flin,status='old')
	write(*,'(//a)')' Enter Output File Name : '
	read(*,'(a)')fout
	open(5,file=fout,status='unknown')

	call input1(mr,cr,kr,phi,tt,ft,nd,nfor,ndim,delt,ttot,xmin,vmin)

	open(2,file='x.out',status='unknown')
	open(3,file='x1.out',status='unknown')
	open(4,file='x2.out',status='unknown')

*******************************
*       Generate Modal Force Vectors at discrete time intervals
*       fr(i,j) contains the modal force vector Fr at time step 'i'
*               j=1,nd corresponds to the DOF at which fr is applied
*       {fr} = [Phi]'.{ft}
*       where, {ft} = physical force vector at time step 'i'
*              [Phi] = modal transformation matrix for given system
*              [Phi]' = transpose of [Phi]
*******************************
	call force1(tt,ft,delt,phi,nd,fr,nfor,ndim,ttot)

*******************************
*       Generate MODAL time history for all DOF's
*       q(i,j), q1(i,j), q2(i,j) contain the modal
*       displacements, velocities and accelerations resp.
*         i --> corresponds to time step number
*         j --> corresponds to a particular DOF
*******************************
*       Direct Integration of the uncoupled equations of motion 
*       are carried out using the Newmark-B method
*******************************
	do 10 i=1,nd
	 call newb(i,mr(i,1),cr(i,1),kr(i,1),fr,q,q1,q2,delt,ndim,ttot,
     1   xmin(i,1),vmin(i,1))
10      continue

*******************************
*       Transform MODAL time histories back to PHYSICAL Coords.
*       {x(t)} = [Phi].{q(t)}  ..etc.
*       Print Output
*******************************

	ic=1
	do 20 t=0.,ttot,delt

	  do 30 i=1,nd
		qt(i,1)=q(ic,i)
		q1t(i,1)=q1(ic,i)
		q2t(i,1)=q2(ic,i)
30        continue

	  call matmul(phi,qt,nd,mr,ndim)
	  write(2,43)t,(mr(nn,1),nn=1,nd)
43        format(f8.4,10e14.6)
	  do 40 i=1,nd
	    q(ic,i)=mr(i,1)
40        continue

	  call matmul(phi,q1t,nd,mr,ndim)
	  write(3,43)t,(mr(nn,1),nn=1,nd)
	  do 41 i=1,nd
	    q1(ic,i)=mr(i,1)
41        continue

	  call matmul(phi,q2t,nd,mr,ndim)
	  write(4,43)t,(mr(nn,1),nn=1,nd)
	  do 42 i=1,nd
	    q2(ic,i)=mr(i,1)
42        continue

	ic=ic+1
20      continue

*************************
	call res_spec(ic-1,q,q1,q2,delt,ndim,nd)
*************************
*       Generate values of MAXIMUM RESPONSE
*       @ each DOF
*************************

	stop
	end

*************************
	subroutine input1(m,c,k,phi,t,f,nd,nfor,ndim,delt,ttot,xmin,
     1  vmin)
*************************
*       [m1] is the system mass matrix (input)
*       [k1] is the system stiffness matrix (input)
*       [m] contains the modal mass parameters (generated)
*       [c] contains the modal damping parameters (generated)
*       [k] contains the modal stiffness parameters (generated)
*       alpha, beta = parameters to generate proportional damping matrix
*         [c] = alpha[m] + beta[k]
*       [phi] contains the modal transformation matrix
*       {t} contains the times at which the load vector
*       {f} is specified
*       nd = the number of degrees of freedom
*       nsol = 1 for EIGENSOLUTION ONLY
*       nsol = 2 for EIGENSOLUTION + MODE SUPERPOSITION ANALYSIS
*       delt = step size
*       ttot = total time for which response will be evaluated
*************************
*	This program uses the IMSL Math Subroutine DGVCSP to evaluate
*	the eigenvalues and eigenvectors of the dynamic system
*************************

	implicit real *8 (a-h,o-z)
	real *8 m,k,m1,k1
	external dgvcsp
	dimension m(ndim,1),c(ndim,1),k(ndim,1),phi(ndim,ndim),t(20,10)
	dimension nfor(10),f(20,10),m1(10,10),c1(10,10)
	dimension k1(10,10),omega(10)
	dimension phit(10,10),eva(10),evec(10,10)
	dimension temp1(10,10),temp2(10,10),temp3(10,10),xin(10,1),
     1  vin(10,1),xmin(10,1),vmin(10,1)

	read(1,*)nsol,nd       
	
	do 90005 i=1,nd
	    read(1,*)(m1(i,j),j=1,nd)
90005   continue

	do 90006 i=1,nd
	    read(1,*)(k1(i,j),j=1,nd)
90006   continue

********************
	do 825 i=1,nd
	  do 825 j=1,nd
	    temp1(i,j)=m1(i,j)
	    temp2(i,j)=k1(i,j)
825     continue
	
c	Call the IMSL routine DGVCSP to evaluate the eigenvalues and
c	eigenvectors of the system.

	call dgvcsp(nd,temp2,ndim,temp1,ndim,eva,evec,ndim)

	do 903 i=1,nd
	  omega(i)=sqrt(eva(nd-i+1))
903     continue

c	write(*,*)(omega(i),i=1,nd)
c	pause
c	stop

	do 904 i=1,nd
	  do 904 j=1,nd
		phi(i,j)=evec(i,nd-j+1)
904     continue

******************************
*       Normalize Eigenvectors wrt 1st row values
******************************

	do 826 j=1,nd
	  temp=phi(1,j)
	  do 826 i=1,nd
	    phi(i,j)=phi(i,j)/temp
826     continue

******************************
*       Print EIGEN OUTPUT
*       Terminate if NSOL=1
******************************

	write(5,905)
905     format(/,' E I G E N S O L U T I O N',//,
     1  ' Mode #        Omega(rad/sec)',/,
     2  ' ----------------------------')
	do 906 i=1,nd
	  write(5,907)i,omega(i)
907       format(i6,e22.8)
906     continue

	write(5,912)
912     format(//,' Mode Shapes',/)
	do 908 i=1,nd
	  write(5,'(10E14.6)')(phi(i,j),j=1,nd)
908     continue

	if(nsol.eq.1)stop

********************
********************

	read(1,*)delt,ttot
	read(1,*)alpha,beta

	call trans(phi,phit,nd,ndim)
	call matmul1(phit,m1,temp1,nd,ndim)
	call matmul1(temp1,phi,temp2,nd,ndim)
	do 909 i=1,nd
	  m(i,1)=temp2(i,i)
909     continue

	call matmul1(phit,k1,temp1,nd,ndim)
	call matmul1(temp1,phi,temp2,nd,ndim)
	do 910 i=1,nd
	  k(i,1)=temp2(i,i)
910     continue

	do 911 i=1,nd
	  zeta=0.5*(alpha/omega(i)+beta*omega(i))
	  c(i,1)=2.0*m(i,1)*omega(i)*zeta
911     continue

	do 921 i=1,nd
	  do 921 j=1,nd
	    c1(i,j)=alpha*m1(i,j)+beta*k1(i,j)
921     continue

	write(5,922)nd
922     format(//,' S Y S T E M  P R O P E R T I E S',
     1  //,' Degrees of Freedom =',i3,//,' Mass Matrix',/)
	do 913 i=1,nd
	  write(5,914)(m1(i,j),j=1,nd)
914       format(10e15.6)
913     continue

	write(5,915)
915     format(//,' Stiffness Matrix',/)
	do 916 i=1,nd
	  write(5,914)(k1(i,j),j=1,nd)
916     continue


	write(5,917) 
917     format(//,' Damping Matrix',/)
	do 918 i=1,nd
	  write(5,914)(c1(i,j),j=1,nd)
918     continue


	write(5,923)(m(i,1),i=1,nd)
923     format(//,' Modal Masses',//,10e15.6)

	write(5,924)(k(i,1),i=1,nd)
924     format(//,' Modal Stiffnesses',//,10e15.6)

	write(5,925)(c(i,1),i=1,nd)
925     format(//,' Modal Damping',//,10e15.6)

	read(1,*)(xin(i,1),i=1,nd)
	read(1,*)(vin(i,1),i=1,nd)

	do 851 i=1,nd
	  do 851 j=1,nd
	    temp1(i,j)=0.0
	    temp2(i,j)=0.0
	    temp3(i,j)=0.0
851     continue

	do 852 i=1,nd
	  temp1(i,i)=1.0/m(i,1)
852     continue

	call matmul1(temp1,phit,temp2,nd,ndim)
	call matmul1(temp2,m1,temp3,nd,ndim)
	call matmul(temp3,xin,nd,xmin,ndim)
	call matmul(temp3,vin,nd,vmin,ndim)

*********************************
*       Input of forcing Function
*********************************

	do 90011 i=1,nd
	  nfor(i)=2
	  t(1,i)=0.0
	  t(2,i)=ttot+100.
	  f(1,i)=0.0
	  f(2,i)=0.0
90011   continue

90014   read(1,*,err=90012)i,nfor(i)
	do 90013 j=1,nfor(i)
	  read(1,*)t(j,i),f(j,i)
90013   continue
	goto 90014
90012   continue

	return
	end

**************************
	subroutine force1(tt,f,dt,p,nd,fr,nfor,ndim,ttot)
**************************

	implicit real *8 (a-h,o-z)
	dimension tt(20,10),f(20,10),fr(6000,ndim),p(ndim,ndim),
     1  f1(10,1),pt(10,10),nfor(ndim),ft(10,1)

	do 90030 i=1,nd
	  do 90030 j=1,nd
		pt(i,j)=p(j,i)
90030   continue

	ic=1
	do 90020 t=0.,ttot,dt

*******************************
*       Interpolate to obtain value of forcing function Ft
*       at time = t
*******************************
	  do 90022 nf=1,nd
	   do 90024 n=1,nfor(nf)-1
	    if(t.ge.tt(n,nf).and.t.le.tt(n+1,nf))then
		ft(nf,1)=f(n,nf)+(f(n+1,nf)-f(n,nf))*
     1          (t-tt(n,nf))/(tt(n+1,nf)-tt(n,nf))
		goto 90022
	    endif
90024      continue
90022     continue

*******************************
*       Generate MODAL force vector at time = t
*******************************
	  call matmul(pt,ft,nd,f1,ndim)
	  do 90026 ii=1,nd
	    fr(ic,ii)=f1(ii,1)
90026     continue

	ic=ic+1
90020   continue

	return
	end

*******************************
	subroutine newb(ii,m,c,k,f,x,x1,x2,del,ndim,ttot,xin,vin)
*******************************
*       Direct Integration of the uncoupled equations of motion
*       using the NEWMARK-b method
******************************
	implicit real *8 (a-h,o-z)
	real *8 m,k,ks
	dimension f(6000,ndim),x(6000,ndim),x1(6000,ndim),x2(6000,ndim)

******************************
*       Assign Initial Conditions (@ t=0.0)
******************************
	x(1,ii)=xin
	x1(1,ii)=vin
	x2(1,ii)=(f(1,ii)-c*x1(1,ii)-k*x(1,ii))/m
	ks=k+2.0*c/del+4.0*m/del**2

******************************
*       Obtain time history at discrete time intervals of del
******************************
	ic=1
	do 101 t=0.,ttot,del
	  dps=f(ic+1,ii)+m*(4.0*x(ic,ii)/del**2+4.0*x1(ic,ii)
     1    /del+x2(ic,ii))+c*(2.0*x(ic,ii)/del+x1(ic,ii))
	  x(ic+1,ii)=dps/ks
	  x2(ic+1,ii)=(4.0/del**2)*(x(ic+1,ii)-x(ic,ii)-
     1    del*x1(ic,ii))-x2(ic,ii)
	  x1(ic+1,ii)=(2.0/del)*(x(ic+1,ii)-x(ic,ii))-x1(ic,ii)

c       write(5,105)ic,t,f(ic+1,ii),f(ic,ii),dps,dx
c105    format('**',i6,e17.7,/,5x,2f15.9,/,5x,2e20.10)

	ic=ic+1
101     continue

	return
	end

*****************************
	subroutine matmul(a,b,nd,c,ndim)
******************************
*       Used to multiply two matrices
*       [C] = [A].{B}
*****************************
	implicit real *8 (a-h,o-z)
	dimension a(ndim,ndim),b(ndim,1),c(ndim,1)

	do 9030 i=1,nd
	  do 9030 j=1,1
	    c(i,j)=0.0
	    do 9030 ii=1,nd
		c(i,j)=c(i,j)+a(i,ii)*b(ii,j)
9030      continue

	return
	end

*****************************
	subroutine res_spec(n,x,x1,x2,del,ndim,nd)
*****************************
*       Obtain the maximum value of the response
*       (displacement, velocity or acceleration)
*       at each DOF for the time span considered
****************************

	implicit real *8 (a-h,o-z)
	dimension x(6000,ndim),x1(6000,ndim),x2(6000,ndim)
	dimension xl(10),x1l(10),x2l(10),t(10),t1(10),t2(10)

	do 9070 i=1,nd
	  xl(i)=x(1,i)
	  x1l(i)=x1(1,i)
	  x2l(i)=x2(1,i)
	  t(i)=0.
	  t1(i)=0.
	  t2(i)=0.
9070    continue

	do 9060 i=2,n

	  do 9060 j=1,nd

	    if(abs(x(i,j)).gt.abs(xl(j)))then
		xl(j)=x(i,j)
		t(j)=float(i-1)*del
	    endif

	    if(abs(x1(i,j)).gt.abs(x1l(j)))then
		x1l(j)=x1(i,j)
		t1(j)=float(i-1)*del
	    endif

	    if(abs(x2(i,j)).gt.abs(x2l(j)))then
		x2l(j)=x2(i,j)
		t2(j)=float(i-1)*del
	    endif

9060    continue

	write(5,9068)
9068    format(///,'    R E S P O N S E   S P E C T R A',///,
     1  '   DOF         Maxm.    @ t          Maxm.    @ t      ',
     2  '   Maxm.    @ t',/,
     3  '                 X                   Vel.          ',
     4  '       Accln.',/,
     5  '   ----------------------------------------------------',
     6  '----------------')

	do 9065 i=1,nd
	  write(5,9066)i,xl(i),t(i),x1l(i),t1(i),x2l(i),t2(i)
9066      format(i5,e15.7,f7.3,e15.7,f7.3,e15.7,f7.3)
9065    continue

	return
	end

*************************************
	subroutine matmul1(a,b,c,nd,ndim)
*************************************

	implicit real*8 (a-h,o-z)
	dimension a(10,10),b(10,10),c(10,10)

	do 80010 i=1,nd
	  do 80010 j=1,nd
	    c(i,j)=0.0
	    do 80010 k=1,nd
	      c(i,j)=c(i,j)+a(i,k)*b(k,j)
80010   continue

	return
	end

***************************************
	subroutine trans(p,pt,nd,ndim)
***************************************

	implicit real*8 (a-h,o-z)
	dimension p(10,10),pt(10,10)

	do 80001 i=1,nd
	  do 80001 j=1,nd
	    pt(i,j)=p(j,i)
80001   continue

	return
	end
***********************************************