plas.for

一个组合程序

**********************************
*       Program elas_plas.f
*       Simulates the NON-LINEAR RESPONSE of a SDOF using
*       an ELASTO-PLASTIC hysteresis loop to model the 
*       Spring Resistance Function, Rm
**********************************
*       The program uses the Newmark-B integration scheme
*       to evaluate the response at each discrete time step
**********************************
*       Program Author : Sanjoy Chakraborty
**********************************

	program resp_nonlin
	implicit real *8 (a-h,o-z)
	real *8 k,m
	character *80 header
	dimension x(50000),x1(50000),x2(50000),ft(20),tt(20),
     1  p(50000)

	write(*,10)
10      format(//,' Elasto-Plastic Response of a SDOF',//,
     1  ' Program Author : Sanjoy Chakraborty, Auburn University',//)


	call input(k,m,c,rm,tt,ft,nfor,tend,h,header)

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

	call resp(x,x1,x2,m,c,k,h,tend,rm,p,ic)

	call output(m,c,k,rm,tt,ft,nfor,tend,h,ic,x,x1,x2,header)

	stop
	end

**********************************
	subroutine input(k,m,c,rm,tt,ft,nfor,tend,h,header)
**********************************
*       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 *22 fin
	character *80 header
	dimension tt(20),ft(20)

	write(*,'(//a)')' Input file :'
	read(*,'(a)')fin
	open(1,file=fin,status='old')

	read(1,'(a80)')header
	read(1,*)k,m,c,rm,tend,h
	read(1,*)nfor

	do 110 i=1,nfor
	  read(1,*)tt(i),ft(i)
110     continue

	close(1)

	return
	end

**********************************
	subroutine inter_f(nf,tt,ft,p,h,tend)
**********************************
*       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(20),ft(20),p(50000)

	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

	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(50000),x1(50000),x2(50000),p(50000)

	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(50000),x(50000),x1(50000),x2(50000)

	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,header)
**********************************
*       Outputs the computed results to 2 files :
*       'fout' -> system parameters and maximum response values
*       x.out -> time vs. displacement data
**********************************

	implicit real *8 (a-h,o-z)
	real *8 m,k
	character *22 fout
	character *80 header
	dimension tt(20),ft(20),x(50000),x1(50000),x2(50000)

	write(*,'(//a)')' Output file :'
	read(*,'(a)')fout
	open(2,file=fout,status='unknown')
	open(3,file='x.out',status='unknown')

	do 601 i=1,ic-1
	  t=float(i-1)*h
	  write(3,609)t,x(i),x1(i),x2(i)
609       format(4e16.8)
601     continue

	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

	write(2,603)header,m,c,k,rm,h,tend
603     format(/,10x,
     1  ' NON-LINEAR / ELASTO-PLASTIC RESPONSE OF A SDOF',/,10x,
     1  '   Direct Integration using Newmark-B Method',/,10x,
     1  '      Program Author : Sanjoy Chakraborty',//,
     1  //,a80,//,' Mass      :',e16.6,/,
     2            ' Damping   :',e16.6,/,
     3            ' Stiffness :',e16.6,/,
     4            ' Rm        :',e16.6,/,
     5            ' Time step :',e16.6,//,
     5  ' Response was calculated from : 0.0  to  ',f5.2,' secs')


	write(2,604)xm,txm,x1m,tx1m,x2m,tx2m
604     format(//,' Values of Maximum response :',//,
     1  ' Maxm. Displacement  =',e15.6,'  at t=',e14.6,/,
     1  ' Maxm. Velocity      =',e15.6,'  at t=',e14.6,/,
     1  ' Maxm. Acceleration  =',e15.6,'  at t=',e14.6)



	return
	end