
C -- start program ammonia --------------------------------------------

C This program steps down the PRESSURE
C    Written by: John Morud (1992-95)
C    Kinetics changed by: G. Juonys (1997)
C    Step down PRESSURE: S. Skogestad (1997)
C The results are stored in the file datasim and can be plotted with MATLAB
C See also file README

      program ammonia
      parameter (Npoint=30,Nbed1=10,Nbed2=10,Nbed3=10)
      real mc, mh, Cp, Cpc, Q1, Q2, Q3,  
     & Tf, cf, dm1, dm2, dm3, dHrx
      common /crap/ mc, mh, Cp, Cpc, Q1, Q2, Q3, 
     & Tf, cf, dm1, dm2, dm3, dHrx
	  common /pres/ p
	  real p
      real T(Npoint), dt,f(Npoint), Tp(Npoint),fp(Npoint)
      real tt, rhoc, vol1, vol2, vol3, time
      real steptime, Tfstep(10), stepindex
      integer i,j,k
	  real Tinit(Npoint)
      data Tinit/369.4310, 376.3844, 383.8233, 391.8714, 400.6834, 
     & 410.4583, 421.4572, 434.0176, 448.5350, 465.2707, 444.9769, 
     & 461.0420, 478.4799, 494.4777, 505.6704, 511.7069, 514.4583, 
     & 515.6106, 516.0755, 516.2603, 498.7332, 506.8531, 509.9740, 
     & 511.0447, 511.3971, 511.5115, 511.5485, 511.5604, 511.5642, 
     & 511.5655/

      open(unit=15,file='datasim.m',status='new')
      
	 Tf=250
	 p=200
C  Final simulation time (min)
	  tsim = 150
C  Time between each step in temperature
      tstep = 20

	  tsims = (tsim*60)/10
	  tsteps = tstep*60
C  Initial parameters:
      steptime=-10.
      stepindex=1
      Tfstep(1)=170.
      Tfstep(2)=150.
      Tfstep(3)=150.
      Tfstep(4)=150.
      Tfstep(5)=150.
      Tfstep(6)=150.
      Tfstep(7)=200.
      Tfstep(8)=Tfstep(7)
      Tfstep(9)=Tfstep(8)
      Tfstep(10)=Tfstep(9)

      mc=127.*1000./3600.
      mh=252.*1000./3600.
      Cp=3500.
      Cpc=1100.
      Q1=58.*1000./3600.
      Q2=35.*1000./3600.
      Q3=32.*1000./3600.
      p=Tfstep(stepindex)
      cf=0.08
      rhoc=2200.
      vol1=6.69
      vol2=9.63
      vol3=15.2
      dHrx=2.7e6
      dm1=vol1/Nbed1*rhoc
      dm2=vol2/Nbed2*rhoc
      dm3=vol3/Nbed3*rhoc
C Start at upper steady-state for Tf=250C
      do 100 i=1,Npoint
        T(i)=Tinit(i)
  100 continue
      dt=0.1
      do 500 k=1,tsims
C Step in temperature:
      tt=100.*(k-1)*dt
      steptime=steptime+10.
      if (steptime .GE. tsteps) then
        steptime=0.
        stepindex=stepindex+1
C Step the PRESSURE
        p=Tfstep(stepindex)
      endif
      write(*,*) tt, T(1), T(Npoint)      
C Data to file 'datasim.m':
      time=(tt)/60.
      write(15,*) 'tsim(',k,')=',time,';' 
      write(15,*) 'TTsim(:,',k,')=['
         do 550 i=1,Npoint
      write(15,*) T(i)              
  550 continue
      write(15,*) '];'
C Start integration:
C Predictor:
      do 200 i=1,100       
        call righthand(T, fp)       
        do 300 j=1,Npoint
          Tp(j)=T(j)+fp(j)*dt
  300   continue
C Corrector:
       call righthand(Tp,f)
       do 400 j=1,Npoint
         T(j)=T(j)+0.5*(fp(j)+f(j))*dt
  400  continue
C End of integration:
  200  continue
  500  continue
       end

C --------------------------------------------------------------------
      real function hx(To)
      parameter (Npoint=30,Nbed1=10,Nbed2=10,Nbed3=10)
      real mc, mh, Cp, Cpc, Q1, Q2, Q3,  
     & Tf, cf, dm1, dm2, dm3, dHrx
      common /crap/ mc, mh, Cp, Cpc, Q1, Q2, Q3, 
     & Tf, cf, dm1, dm2, dm3, dHrx
      real To, UA, cstar, ntu, eps, Q, Tho, Ti
      UA=536.*283. 
      cstar=mc/mh
      ntu=UA/(mc*Cp)
      eps=(1-exp(-ntu*(1-cstar)))/(1-cstar*exp(-ntu*(1-cstar)))
      Q=eps*mc*Cp*(To-Tf)
      Tho=To-Q/(mh*Cp)
      Ti=Tf+Q/(mc*Cp)
      hx=Ti
      end

C --------------------------------------------------------------------
      real function rx(T,c)
      parameter (Npoint=30,Nbed1=10,Nbed2=10,Nbed3=10)
      real mc, mh, Cp, Cpc, Q1, Q2, Q3,  
     & Tf, cf, dm1, dm2, dm3, dHrx
      common /crap/ mc, mh, Cp, Cpc, Q1, Q2, Q3, 
     & Tf, cf, dm1, dm2, dm3, dHrx
      common /pres/ p
      real mNH3,mH2,mN2,nNH3,nH2,nN2,rate,T,c,x
      real R, p, pnh3,pn,ph,k1,k2
      mNH3=c 
      mH2=(1-c)*6./34. 
      mN2=(1-c)*28./34. 
      nNH3=mNH3/17. 
      nH2=mH2/2. 
      nN2=mN2/28. 
      x=nNH3/(nNH3+nH2+nN2) 
      R=8.31
      k1=1.79E+4*exp(-87090./(R*(T+273.))) 
      k2=2.57E+16*exp(-198464./(R*(T+273.)))
      pnh3=x*p  
      pn=(1-x)*0.25*p 
      ph=(1-x)*0.75*p
      rate=k1*pn*(ph**3/pnh3**2)**0.5 - k2*(pnh3**2/ph**3)**0.5 
      rx=rate*34./2200./3600. 
      rx=rx*4.75
      end


C --------------------------------------------------------------------
      subroutine righthand(T,f)
     parameter (Npoint=30,Nbed1=10,Nbed2=10,Nbed3=10)
      real mc, mh, Cp, Cpc, Q1, Q2, Q3,  
     & Tf, cf, dm1, dm2, dm3, dHrx
      common /crap/ mc, mh, Cp, Cpc, Q1, Q2, Q3, 
     & Tf, cf, dm1, dm2, dm3, dHrx
      real T(Npoint), f(Npoint),m1,m2,m3,c(Npoint)
      real cin2,cin3,Tin1,Tin2,Tin3,dTdt(Npoint)
      integer i
      external hx, rx
      real hx, rx
comment: For simplicity we use rx((T(i),c(i-1)) rather than rx((T(i),c(i)) 
C   Mass Flows:
      m1=mc+Q1
      m2=m1+Q2 
      m3=m2+Q3 
C Concentration:
C Conc.  in bed 1:
      c(1)=cf+dm1*rx(T(1),cf)/m1 
      do 100 i=2,Nbed1
        c(i)=c(i-1)+dm1*rx(T(i),c(i-1))/m1 
  100 continue
C Quench for bed 2: 
      cin2=m1/m2*c(Nbed1)+Q2/m2*cf 
C Conc. in bed 2
      c(Nbed1+1)=cin2+dm2*rx(T(Nbed1+1),cin2)/m2 
      do 200 i=Nbed1+2,Nbed1+Nbed2
        c(i)=c(i-1)+dm2*rx(T(i),c(i-1))/m2 
  200 continue
C Quench 3:
      cin3=m2/m3*c(Nbed2+Nbed1)+Q3/m3*cf 
C Conc. in bed 3
      c(Nbed1+Nbed2+1)=cin3+dm3*rx(T(Nbed1+Nbed2+1),cin3)/m3 
      do 300 i=Nbed1+Nbed2+2, Nbed1+Nbed2+Nbed3
        c(i)=c(i-1)+dm3*rx(T(i),c(i-1))/m3 
  300 continue
C Temperature:
C Quench 1:
      Tin1=mc/m1*hx(T(Npoint))+Q1/m1*Tf 
C dTdt in bed 1:
      dTdt(1)=(m1*Cp*(Tin1-T(1))+dm1*rx(T(1),cf)*dHrx)/(dm1*Cpc) 
      do 400 i=2,Nbed1
        dTdt(i)=(m1*Cp*(T(i-1)-T(i))+dm1*rx(T(i),c(i-1))*dHrx)/(dm1*Cpc) 
  400 continue
C Quench 2:
      Tin2=m1/m2*T(Nbed1)+Q2/m2*Tf 
C dTdt in bed 2:
      dTdt(Nbed1+1)=(m2*Cp*(Tin2-T(Nbed1+1))
     &             +dm2*rx(T(Nbed1+1),cin2)*dHrx)/(dm2*Cpc) 
      do 500 i=Nbed1+2, Nbed1+Nbed2
        dTdt(i)=(m2*Cp*(T(i-1)-T(i))
     &         +dm2*rx(T(i),c(i-1))*dHrx)/(dm2*Cpc) 
  500 continue
C Quench 3:
      Tin3=m2/m3*T(Nbed1+Nbed2)+Q3/m3*Tf 
C dTdt i bed 3:
      dTdt(Nbed1+Nbed2+1)=(m3*Cp*(Tin3-T(Nbed1+Nbed2+1))
     &              +dm3*rx(T(Nbed1+Nbed2+1),cin3)*dHrx)/(dm3*Cpc) 
      do 600 i=Nbed1+Nbed2+2,Npoint
        dTdt(i)=(m3*Cp*(T(i-1)-T(i))
     &         +dm3*rx(T(i),c(i-1))*dHrx)/(dm3*Cpc) 
  600 continue
      do 700 i=1,Npoint
        f(i)=dTdt(i)
  700 continue 
      end

C -- end program ammonia --------------------------------------------


