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 --------------------------------------------