lresp.for

一个组合程序

**********************************
*       Program res_spec.f
*       Generates the RESPONSE SPECTRA of a LINEAR
*       Undamped SDOF System to any arbitrary applied load
**********************************
*       The program uses the Newmark-B integration scheme
*       to evaluate the response at each discrete time step
**********************************
*       Program Author : Sanjoy Chakraborty : 09CE
**********************************

	program res_spec
	implicit real *8 (a-h,o-z)
	real *8 k,m
	character *12 fout
	dimension x(200000),x1(200000),x2(200000),ft(100),tt(100),
     1  p(200000)
	data pi/3.141592654/

	call input(m,c,rm,tt,ft,nfor,tend,h,td,
     1  trat1,trat2,trath)

	call inter_f(nfor,tt,ft,p,h,tend,pmax)

**********************************
*       Loop for computing the value of maximum DISPLACEMENT 
*       response for successive values of system TIME-Period
**********************************
	
	write(*,'(//a)')' Output File Name :'
	read(*,'(a)')fout
	open(2,file=fout,status='unknown')

	do 701 trat=trat1,trat2,trath
	  t=td/trat
	  k=m*(2.*pi/t)**2
	  cc=c*sqrt(2.*k*m)
	  call resp(x,x1,x2,m,cc,k,h,tend,rm,p,ic)
	  call output(m,cc,k,rm,tt,ft,nfor,tend,h,ic,x,x1,x2,
     1    trat,t,pmax)
701     continue

	stop
	end

**********************************
	subroutine input(m,c,rm,tt,ft,nfor,tend,h,
     1  td,trat1,trat2,trath)
**********************************
*       Used to input the system parameters and forcing 
*       function from the input file fin
**********************************

	implicit real *8 (a-h,o-z)
	real *8 k,m
	character *12 fin
	dimension tt(100),ft(100)
	rm=0.1e12

	write(*,115)
115     format(//,' Program to Generate Load Response Spectra',
     1  //,' Program Author : Sanjoy Chakraborty, Auburn University',
     2  ///,' Input File Name :')
	read(*,'(a)')fin
	open(1,file=fin,status='old')
	
	read(1,*)m,c,tend,h
	read(1,*)td,trat1,trat2,trath

	nfor=1
110     continue
	  read(1,*,err=111)tt(nfor),ft(nfor)
	  nfor=nfor+1
	goto 110
111     nfor=nfor-1

	close(1)

	return
	end

**********************************
	subroutine inter_f(nf,tt,ft,p,h,tend,pmax)
**********************************
*       Used to interpolate for the forcing function values
*       at the times at which the solution is sought
**********************************

	implicit real *8 (a-h,o-z)
	dimension tt(100),ft(100),p(200000)

	ic=1
	do 1001 t=0.,tend,h
	  do 1002 i=1,nf-1
	    if(t.ge.tt(i).and.t.le.tt(i+1))then
		p(ic)=ft(i)+(ft(i+1)-ft(i))*(t-tt(i))/(tt(i+1)-tt(i))
		goto 1003
	    endif
1002      continue
1003    ic=ic+1
1001    continue

	pmax=abs(ft(1))
	do 1005 i=2,nf
	  if(abs(ft(i)).gt.pmax)pmax=abs(ft(i))
1005    continue

	return
	end

**********************************
	subroutine resp(x,x1,x2,m,c,k,h,tend,rm,p,ic)
**********************************
*       Computes the system response (displacement, velocity
*       and acceleration) at each time step
**********************************
*       The loop (=i) parameter is used to determine the 
*       exact location of the system on the hysteretic loop
*       describing the resistance function Rm
*       loop =1  -->  elastic loading stage
*       loop =2  -->  plastic loading stage
*       loop =3  -->  elastic rebound stage
*       loop =4  -->  plastic rebound stage
*       xmax, xmin, and xlim are parameters used to control
*       the physical route of the resistance function along
*       the hysteretic loop
*       xmax = max. +ve plastic deformation (loop=2)
*       xmin = max. -ve plastic deformation (loop=4)
**********************************

	implicit real *8 (a-h,o-z)
	real *8 k,m,kelas,kplas
	dimension x(200000),x1(200000),x2(200000),p(200000)

	x(1)=0.
	x1(1)=0.
	x2(1)=p(1)/m

	xel=rm/k
	a1=3./h
	a2=6./h
	a3=h/2.
	a4=6./h**2
	kelas=k+a4*m+a1*c
	kplas=a4*m+a1*c
	xlim=xel
	xmin=-xel
	loop=1
	ic=2

	do 501 t=h,tend,h

	  if(loop.eq.1)then
	    call res1(kelas,p,x,x1,x2,m,c,ic,a2,a3,a1)
	    r=-rm-(xmin-x(ic))*k
	    x2(ic)=(p(ic)-c*x1(ic)-r)/m
	    if(x(ic).ge.xlim)then
	      loop=2
	    endif
	    ic=ic+1
	    goto 501

	  elseif(loop.eq.2)then
	    call res1(kplas,p,x,x1,x2,m,c,ic,a2,a3,a1)
	    r=rm
	    x2(ic)=(p(ic)-c*x1(ic)-r)/m
	    if(x1(ic).le.0.)then
	      loop=3
	      xmax=x(ic)
	      xlim=x(ic)-2.*xel
	    endif
	    ic=ic+1
	    goto 501

	  elseif(loop.eq.3)then
	    call res1(kelas,p,x,x1,x2,m,c,ic,a2,a3,a1)
	    r=rm-(xmax-x(ic))*k
	    x2(ic)=(p(ic)-c*x1(ic)-r)/m
	    if(x(ic).le.xlim)then
	      loop=4
	    endif
	    ic=ic+1
	    goto 501

	  elseif(loop.eq.4)then
	    call res1(kplas,p,x,x1,x2,m,c,ic,a2,a3,a1)
	    r=-rm
	    x2(ic)=(p(ic)-c*x1(ic)-r)/m
	    if(x1(ic).ge.0.)then
	      loop=1
	      xlim=x(ic)+2.*xel
	      xmin=x(ic)
	    endif
	    ic=ic+1
	    goto 501

	  endif

501     continue

	return
	end

**********************************
	subroutine res1(k,p,x,x1,x2,m,c,ic,a2,a3,a1)
**********************************
*       Used to compute displacement and velocity
*       at any time step
**********************************

	implicit real *8 (a-h,o-z)
	real *8 k,m
	dimension p(200000),x(200000),x1(200000),x2(200000)

	dps=p(ic)-p(ic-1)+x1(ic-1)*(a2*m+3.*c)+x2(ic-1)*
     1      (3.*m+a3*c)
	dx=dps/k
	dx1=a1*dx-3.*x1(ic-1)-a3*x2(ic-1)
	x(ic)=x(ic-1)+dx
	x1(ic)=x1(ic-1)+dx1

	return
	end

**********************************
	subroutine output(m,c,k,rm,tt,ft,nfor,tend,h,ic,
     1  x,x1,x2,trat,t,pmax)
**********************************
*       Computes the value of Maximum Response at each value
*       of TIME PERIOD of the system, then writes the output
*       to Unit:3 in the format Td/T vs. (DLF)max
**********************************

	implicit real *8 (a-h,o-z)
	real *8 m,k
	dimension tt(100),ft(100),x(200000),x1(200000),x2(200000)

	xm=x(1)
	txm=0.
	x1m=x1(1)
	tx1m=0.
	x2m=x2(1)
	tx2m=0.

	do 605 i=2,ic-1
	  if(abs(x(i)).gt.abs(xm))then
	    xm=x(i)
	    txm=float(i-1)*h
	  endif
	  if(abs(x1(i)).gt.abs(x1m))then
	    x1m=x1(i)
	    tx1m=float(i-1)*h
	  endif
	  if(abs(x2(i)).gt.abs(x2m))then
	    x2m=x2(i)
	    tx2m=float(i-1)*h
	  endif
605     continue

	xstat=pmax/k
	dlf=abs(xm)/xstat
	write(2,'(2e14.7)')trat,dlf

	return
	end