%this code is meant to function exactly as the original fortran cpdheat %code % % % C THIS IS THE CPD-NLG MODEL % % % C % % % C THIS MODEL WAS DEVELOPED BY SANDIA NATIONAL LABORATORIES UNDER % % % C FWP 0709 FOR THE DEPARTMENT OF ENERGY'S PITTSBURGH ENERGY % % % C TECHNOLOGY CENTER AND THE DOE DIVISION OF ENGINEERING AND GEOSCIENCES % % % C THROUGH THE OFFICE OF BASIC ENERGY SCIENCES; % % % C AND BY THE UNIVERSITY OF UTAH THROUGH FUNDING FROM % % % C THE ADVANCED COMBUSTION ENGINEERING RESEARCH CENTER (ACERC), WHICH % % % C IS PRINCIPALLY SPONSORED BY THE NATIONAL SCIENCE FOUNDATION, THE % % % C STATE OF UTAH, AND BY A CONSORTIUM OF INDUSTRIAL COMPANIES. % % % C THE CODE WILL NOT BE FORMALLY LICENSED. NEITHER THE U.S. OR THE % % % C DOE, NOR ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, % % % C OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY, % % % C COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT, OR % % % C PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD INFRINGE PRIVATELY % % % C OWNED RIGHTS. % % % c % % % c The CPD model is intended to solve pyrolysis rate equations % % % c based on a percolative bond-breaking scheme. This version includes the % % % c flash distillation program to distinguish between tar and metaplast. % % % c This program also includes a crosslinking scheme. (January, 1991) % % % c % % % c % % % c Most recent modifications to this model include a (a) nitrogen release model % % % c and (b) a model to break the light gas into species based on a correlation. % % % c These modifications were made by Dominic Genetti in his M.S. Thesis work at % % % c BYU (1999). % % % c % % % c This is the black liquor version of the cpd model that was made for a % % % c constant heating rate. The self-adjusting time step was modified % % % c (about January, 2004). % % % c function log_data=cpdheat(Input_matrix) format long %code in the inputs nmax=10; nt=1; intar=false; idiff=false; ftar0=0; fgas0=0; fchar0=1; ip=30; zero=0; ftar=0; fchar=0; fmet=0; fvol=0; fvolold=0; lib=0; %these next arrays are possibly in need of modification for enhanced %flexibility y=zeros(4,1); yp=zeros(4,1); ypp=zeros(4,1); ypred=zeros(4,1); tim=zeros(50,1); tem=zeros(50,1); ft=zeros(35,1); mt=zeros(35,1); u0=0; kkk=1; umax=0; yield=0; g0=0; g1=0; g2=0; % % % % % % % % % %read in the situation specific inputs % % % % % p0=Input_matrix(1);% !p0 c0=Input_matrix(2);% !c0 sigp1=Input_matrix(3);% !sig+1 mw1=Input_matrix(4);% !mw mdel=Input_matrix(5);% !mdel (7 will be subtracted internally to the CPD model) fnit=0.0119;% !fnit (daf mass fraction of nitrogen in unreacted coal) fst=0.50;% !fst (fraction of nitrogen remaining in the char) fhyd=0.0595;% !fhyd (daf mass fraction of hydrogen in unreacted coal) fcar=0.7406;% !fcar (daf mass fraction of carbon in unreacted coal) foxy=0.1841;% !foxy (daf mass fraction of oxygen in unreacted coal) ab=Input_matrix(12);% !ab eb0=Input_matrix(13);% !eb ebsig=Input_matrix(14);% !ebsig ac=Input_matrix(15);% !ac=rho ec0=Input_matrix(16);% !ec ag=Input_matrix(17);% !ag eg0=Input_matrix(18);% !eg egsig=Input_matrix(19);% !egsig acr=Input_matrix(20);% !Acr (pre-exponential factor for crosslinking rate) ecr=Input_matrix(21);% !Ecr (Activation energy for crosslinking rate) an=Input_matrix(22);% !an (Pre-exponential factor for HCN release) en0=Input_matrix(23);% !en (Activation energy for HCN release) ensig=Input_matrix(24);% !ensig (deviation bound for dist en for HCN release) press=1;% !pressure (atm) % %calculate mdel and other structure coefficiencets, values that should not be changed % mdel=mdel-7; l0=p0-c0; mb=2*mdel; ma=mw1-sigp1*mdel; beta = ma; sig = sigp1-1; finf = 1/(2*ma/(mb*sigp1*(1-c0))+1); rba = mb/ma; fnca0 = fnit*mw1/(mw1-sigp1*mdel); %c calculate O/C and H/C ratios for light gas model xoc = (foxy/16)/(fcar/12); yhc = fhyd/(fcar/12); %c test to see if p0 is above percolation threshhold pthresh = 1/sig; if p0<=pthresh error('Error: value of p0 is below the percolation threshold') end tim(1)=0; % %resume situation specific parameters % tem(1)=300;% !initial temperature (K) ntim=2;% !number of time points; this needs to be an integer % %calculate the number of time spoints % ntim=ntim+1; % %resume situation specific parameters % %The number of values below depend on the number of time points. %They should begin at 2, and one value should be put in for each time point %beyond 2. If the points become numerous and repetetive, this section can %be recoded as a for loop. heat(2)=Input_matrix(6);%heating rate tem(2)=Input_matrix(7);%temperature tim(2)=Input_matrix(8); % %check the value of heat, and recalculate tim % if heat(2)~=0 tim(2)=(tem(2)-tem(2-1))/heat(2)+tim(2-1); end % %resume situation specific inputs % heat(3)=Input_matrix(9); tem(3)=Input_matrix(10); tim(3)=Input_matrix(11); % %check the value of heat, and recalculate tim % if heat(3)~=0 tim(3)=(tem(3)-tem(3-1))/heat(3)+tim(3-1); end %initialize the nt+1 value beyond the final input heat value heat=heat'; heat=[heat;0]; % %resume situation specific inputs % dt0=5e-4;%!dt (s) % %input the print increment % iprint=1; % %resume situation specific parameters % dtmax=1; timax=600;% !timax (maximum residence time [s] for calculations) nmax=20;% !nmax (maximum number of mers for tar molecular wt) % %input a logic value % inside=true; %!true for Light Gas Calculation, otherwise false % % % % % % % % % %end of situation specific parameters % % % % % %initialize the data blocks ftold=zeros(35,1); metold=zeros(35,1); x3=.2; x2=.3; xmat=zeros(36,1); ymat=zeros(36,1); k=zeros(36,1); lmat=zeros(36,1); v=zeros(36,1); xmw=zeros(36,1); tarold=zeros(35,1); f=zeros(36,1); z=zeros(36,1); pv=zeros(36,1); fgas0=0; ip=30; ftart=0; fx=0; l=0; kc=0; kb=0; kg=0; kn=0; kp=0; ft=zeros(35,1); mt=zeros(35,1); mtot=0; sumu=zeros(200,1); VoL=0; Ltot=0; Vtot=0; xmwtot=0; mgas=0; fgas=0; fchar=0; fgas0old=0; % initialize some more variables breakcounter=0; dtold=dt0; itest=0; sumy=0; yygas=zeros(5,1); y(1) = l0; y(2) = 2.*(1.-c0-l0); y(3) = c0; y(4) = fnca0; fnca=y(4); aind0 = l0 + (1.-c0-l0); siginv = 1./sig; pstar = 0.5*siginv; ynhcn = 0.0; yntar = 0.0; yytar = 0.0; yf = 0.0; yyyy=false; %start the calculation loop Rg=1.987;%(cal/mol) iii=0;%initialize the count fgas0 = 0; fcross = 0.0; fntar = 0.0; fntd = 0.0; fnt = 0.0; cmw1 = 0.0; smw1 = 0.0; time=0; if heat(2)~=0 dt0=tim(2)/1000; end dt=min(dt0,dtmax); ntmax=100*(timax/dt+1); for iii=1:ntmax time=time+dt; % CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC % C % C PREDICTOR % C % CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC while 0==0 %this formulation is to translate the "go to" command in fortran if time<=tim(nt+1) Tp=tem(nt)+(time-tim(nt))*(tem(nt+1)-tem(nt))/(tim(nt+1)-tim(nt)); break elseif nt4000 print1='warning, gas temp is too high' Tp end %C-- DEVOLATILIZATION RATES while 0==0%this formulation allows a more direct translation from a 'goto' %statement in the fortran version [ypp,Tp,l0,l,c0,g0,ma,rba,finf,sig,siginv,nmax,pstar,... ab,eb0,ebsig,ac,ec0,ag,eg0,egsig,Rg,... fst,fnca0,fnca,an,en0,ensig,fx,kg,kn,kp,ft,mt,kb]... =perks(y,Tp,l0,l,c0,g0,ma,rba,finf,sig,siginv,nmax,pstar,... ab,eb0,ebsig,ac,ec0,ag,eg0,egsig,Rg,... fst,fnca0,fnca,an,en0,ensig,fx,kg,kn,kp,ft,mt,kb,iii);%call for j=1:4 ypred(j)=y(j)+dt*ypp(j); ypred(j)=max(ypred(j),zero); end [yp,Tp,l0,l,c0,g0,ma,rba,finf,sig,siginv,nmax,pstar,... ab,eb0,ebsig,ac,ec0,ag,eg0,egsig,Rg,... fst,fnca0,fnca,an,en0,ensig,fx,kg,kn,kp,ft,mt,kb]... =perks(ypred,Tp,l0,l,c0,g0,ma,rba,finf,sig,siginv,nmax,pstar,... ab,eb0,ebsig,ac,ec0,ag,eg0,egsig,Rg,... fst,fnca0,fnca,an,en0,ensig,fx,kg,kn,kp,ft,mt,kb,iii);%call dtold=dt; dymax=0; dymin=1; for j=1:4 if y(j)>5e-3 dy1=abs(dt*0.5*(yp(j)+ypp(j))); dymax=max(dymax,dy1); end end dm=fvol-fvolold; bigdm=max(dymax,dm); if bigdm<.0001 itest=1; dt=dt*1.2; elseif bigdm>.0005 itest=2; dt=dtold/2; dy1 end dt=min(dt,dtmax); if dtold==dt | itest>=1 break end end fvolold=fvol; dtn=dt; for j=1:4 y(j) = y(j) + dt*0.5*(yp(j)+ypp(j)); y(j) = max(zero,y(j)); end fracr=1; if (fmet>1e-5)&&(acr>0) ratecr = acr*exp(-ecr/Rg/Tp)*fmet*dt; fracr = 1.-ratecr/fmet; fmet = fmet-ratecr; fcross = fcross+ratecr; if fmet<0 fcross = fcross-ratecr+fmet+ratecr; fmet = 0.0; fracr = 0; end end %c print data if y(1)>1e-5 intar=true; end [mgas,ftar,ftart,fgas,fchar,ft,mt,intar,l0,c0,g0,ma,rba,finf,sig,siginv,nmax,pstar,... ab,eb0,ebsig,ac,ec0,ag,eg0,egsig,Rg,... u0,beta,kkk,umax,yield,press,l,mb,kb,kc,kg,kp,mtot,sumu]... =perkp(y,mgas,ftar,ftart,fgas,fchar,intar,l0,c0,g0,ma,rba,finf,sig,siginv,nmax,pstar,... ab,eb0,ebsig,ac,ec0,ag,eg0,egsig,Rg,... u0,beta,kkk,umax,yield,press,l,mb,kb,kc,kg,kp,ft,mt,mtot,sumu,iii);%call intar=false; tms=time*1e3; if idiff==true fdift = ftar-ftar0; ftar0 = ftar; fdifg = fgas-fgas0; fdifc = -(fchar-fchar0); fchar0 = fchar; else gasmw=rba*ma/2; if fgas>=1e-5 fgasd=fgas-fgas0; fgas0=fgas; imolw=false; if mod(iii,iprint)==0 imolw=true; end [fgas,gasmw,ft,mt,fracr,ftar,fmet,... Tp,press,nmax,imolw,sumy,tms,time,ftold,metold,x3,x2,... xmat,ymat,k,lmat,v,xmw,tarold,f,z,pv,fgas0old,ip,VoL,Ltot,... Vtot,xmwtot... ]=flash(fgas,gasmw,ft,mt,fracr,ftar,fmet,... Tp,press,nmax,imolw,sumy,tms,time,... ftold,metold,x3,x2,xmat,ymat,k,lmat,v,xmw,tarold,f,z,... pv,fgas0old,ip,VoL,Ltot,Vtot,xmwtot,iii); elseif fgas<1e-5 fmet=ftart; ftar=0; end intar=false; fvol=fgas+ftar; fchar=1-fvol; if fvol>.995 return end % if(iii.eq.1)then % write(6,201) % write(2,202) % 201 format(' time(ms) temp fcross labile ftar ', % x'fgas fsolid ftot fmet') % 202 format('c time(ms) temp fcross labile ftar ', % x'fgas fsolid ftot fmet') % endif end l = y(1); del = y(2); c = y(3); fnca = y(4); p = l+c; tarfac = 1; g1 = (2.*(1-p)-del); g2 = 2.*(c-c0); g = g1+g2; del2 = del/2; aind = del2 + l; g12 = g1/2; g22 = g2/2; gtot = g/2; % c % c % C NITROGEN RELEASE CALCULATIONS % % c new average molecular weight cmw1 = ma + (c0 + l + del2)*sigp1*mdel; % c nitrogen content of char and tar fnt = y(4)*ma/cmw1; % c differential mass of nitrogen released with tar yntd = fnt*sumy; % c total nitrogen released with tar yntar = yntar + yntd; % c nitrogen remaining in char ynchar = fchar*fnt; % c differential mass of nitrogen released as light gas ynhcnd = fchar*dt*(-0.5)*(yp(4)+ypp(4))*ma/cmw1; % c total mass of nitrogen released as light gas ynhcn = ynhcn + ynhcnd; % c mass of nitrogen released with tar by difference yntar = fnit - ynhcn - ynchar; % c fraction of original nitrogen remaining in char fnchar = ynchar/fnit; % c fraction of original nitrogen released as tar fntar = yntar/fnit; % c fraction of original nitrogen released as light gas fnhcn = ynhcn/fnit; % c total fractional release of nitrogen fntot = (fnit - fnt*fchar)/fnit; % % % c DISTRIBUTE LIGHT GAS INTO H20, CO2, CO, CH4, & other HC's % % c yf is a CPD indicator of the fraction of total light gas % c that has been released. The look up table on light gas % c composition is based on yf. if inside==true yf=1-aind/aind0; [yygas,inside,lib,yyyy,yf,xoc,yhc]=... lightgas(yygas,inside,lib,yyyy,yf,xoc,yhc,iii);%call % c calculate fraction of total mass release that is h2o, co2, ch4, % c co, and other light gases for ik=1:5 ffgas(ik)=fgas*yygas(ik); end end % c PRINT RESULTS IN OUTPUT FILES % c % if(mod(iii-1,iprint).eq.0)then % write(20,220)tms,l,c,del2,g12,g22,gtot,p % write(21,221)tms,tp,cmw1,y(4),fnt,fnchar,fntar,fnhcn,fntot % if(inside) % x write(22,222)tms,fgas,ffgas(1),ffgas(2),ffgas(3),ffgas(4), % x ffgas(5),yygas(1),yygas(2),yygas(3),yygas(4), % x yygas(5),yf % endif % 220 format(' ',1pe10.3,2x,7(0pf7.5,2x),0pf12.9,0pf9.7,0pf10.5,0pf12.9) % 221 format(' ',2(1pe10.3),0pf8.2,6(0pf10.6)) % 222 format(' ',1pe10.3,0pf9.4,11(0pf9.4)) log_data(iii,1)=ftar; log_data(iii,2)=fvol; log_data(iii,3)=time; log_data(iii,4)=Tp; end intar=false; fvol=fgas+ftar; fchar=1-fvol; if fvol>.995 return end end