#Nonisothermal Reactor Network Synthesis #Constant Density Assumption #CSTRs and CFRs #linear control of sidestream flowrates #quadratic control of CFR temperature with continuous 1st derivatives # # C. A. Schweiger and C. A. Floudas, "Optimization # Framework for the Synthesis of Chemical Reactor # Networks", Ind. Eng. Chem. Res. 1998, submitted # OPTIONS {{ EVTOL = "1e-10"; FTOL = "1e-6"; OTOL = "1e-6"; SNOPT = "Iterations Limit 100000"; SNOPT = "Major Iterations Limit 10000"; MINOS = "Iterations Limit 100000"; MINOS = "Major Iterations Limit 10000"; }} DECLARATIONS {{ INDEX {i,j,k,k',l,l',r,p}; #Sets SET I = |1:4|; #reaction components SET J = |1:4|; #reaction paths SET K = |1:1|; #crossflow reactors (CFRs) SET L = |1:1|; #mixed reactors (CSTRs) SET R = |1|; #feeds SET P = |1|; #products SET CI = |0:9|; #control intervals SET CP = |0:10|; #control points (intervals+1) SET CO = |0:2|; #collocation points SET OR = |0:1|; #order or the collocation equations #Parameters #Gas constant PARAMETER rg = 1.987; #cal/gmol/K #reaction rate constants [1/h,1/h,1/h,1/h] PARAMETER k0(J) = {2.0e13,2.0e13,8.15e17,2.1e5}; #activation energy [cal/gmol] PARAMETER e(J) = {38000,38000,50000,20000}; #heat of reaction (-H/rho/Cp) [K L/gmol] PARAMETER h(J) = {36,129,108,222}; #stoichiometric coefficients PARAMETER nu(I,J) = {-1,-1, 0, 0, 1, 0,-1, 0, 0, 1, 1,-1, 0, 0, 0, 1}; #feed flowrate [L/h] PARAMETER fa(R) = {100}; #feed composition [gmol/L] PARAMETER ca(R,I) = {1,0,0,0}; XVAR {fab(R,K), #flowrate from feed splitters to CFR side mixers fac(R,K), #flowrate from feed splitters to CFR main mixers fad(R,L), #flowrate from feed splitters to CSTR mixers fb(K), #flowrate of CFR side inlets cb(K,I), #composition of CFR side inlets fc(K), #flowrate of CFR main inlets cc(K,I), #composition of CFR main inlets fd(L), #flowrate of CSTR inlets cd(L,I), #composition of CSTR inlets fe(K), #flowrate of CFR main outlets ce(K,I), #composition of CFR main outlets ff(K), #flowrate of CFR side outlets cf(K,I), #composition of CFR side outlets fg(L), #flowrate of CSTR outlets cg(L,I), #composition of CSTR outlets feb(K,K), #flowrate from CFR main splitters to CFR side mixers fec(K,K), #flowrate from CFR main splitters to CFR main mixers fed(K,L), #flowrate from CFR main splitters to CSTR mixers ffb(K,K), #flowrate from CFR side splitters to CFR side mixers ffc(K,K), #flowrate from CFR side splitters to CFR main mixers ffd(K,L), #flowrate from CFR side splitters to CSTR fgb(L,K), #flowrate from CSTR splitters to CFR side mixers fgc(L,K), #flowrate from CSTR splitters to CFR main mixers fgd(L,L), #flowrate from CSTR splitters to CSTR mixers feh(K,P), #flowrate from CFR main splitters to product mixers ffh(K,P), #flowrate from CFR side splitters to product mixers fgh(L,P), #flowrate from CSTR splitters to product mixers fh(P), #flowrate of product ch(P,I), #composition of product vm(L), #CSTR reactor volumes tm(L), #CSTR reactor temperatures rm(L,J), #CSTR reaction rates vt(K), #PFR reactor volumes ks(K,CP), #control parameters for dsdv ko(K,CP), #control parameters for dldv kt(K,CP,OR) #control parameters for dtdv }; #Bounds and starting points LBDS fab(R,K) = {0}; STP fab(R,K) = {0}; UBDS fab(R,K) = {1000}; LBDS fac(R,K) = {0}; STP fac(R,K) = {50}; UBDS fac(R,K) = {1000}; LBDS fad(R,L) = {0}; STP fad(R,L) = {50}; UBDS fad(R,L) = {1000}; LBDS fb(K) = {0}; STP fb(K) = {0}; UBDS fb(K) = {100}; LBDS cb(K,I) = {0}; STP cb(K,I) = {ca(R,I)}; UBDS cb(K,I) = {10}; LBDS fc(K) = {0}; STP fc(K) = {50}; UBDS fc(K) = {1000}; LBDS cc(K,I) = {0}; STP cc(K,I) = {ca(R,I)}; UBDS cc(K,I) = {10}; LBDS fd(L) = {0}; STP fd(L) = {50}; UBDS fd(L) = {1000}; LBDS cd(L,I) = {0}; STP cd(L,I) = {ca(R,I)}; UBDS cd(L,I) = {10}; LBDS fe(K) = {0}; STP fe(K) = {50}; UBDS fe(K) = {1000}; LBDS ce(K,I) = {0}; STP ce(K,I) = {ca(R,I)}; UBDS ce(K,I) = {10}; LBDS ff(K) = {0}; STP ff(K) = {0}; UBDS ff(K) = {1000}; LBDS cf(K,I) = {0}; STP cf(K,I) = {ca(R,I)}; UBDS cf(K,I) = {10}; LBDS fg(L) = {0}; STP fg(L) = {50}; UBDS fg(L) = {1000}; LBDS cg(L,I) = {0}; STP cg(L,I) = {ca(R,I)}; UBDS cg(L,I) = {10}; LBDS feb(K,K) = {0}; STP feb(K,K) = {0}; UBDS feb(K,K) = {1000}; LBDS fec(K,K) = {0}; STP fec(K,K) = {0}; UBDS fec(K,K) = {1000}; LBDS fed(K,L) = {0}; STP fed(K,L) = {0}; UBDS fed(K,L) = {1000}; LBDS ffb(K,K) = {0}; STP ffb(K,K) = {0}; UBDS ffb(K,K) = {1000}; LBDS ffc(K,K) = {0}; STP ffc(K,K) = {0}; UBDS ffc(K,K) = {1000}; LBDS ffd(K,L) = {0}; STP ffd(K,L) = {0}; UBDS ffd(K,L) = {1000}; LBDS fgb(L,K) = {0}; STP fgb(L,K) = {0}; UBDS fgb(L,K) = {1000}; LBDS fgc(L,K) = {0}; STP fgc(L,K) = {0}; UBDS fgc(L,K) = {1000}; LBDS fgd(L,L) = {0}; STP fgd(L,L) = {0}; UBDS fgd(L,L) = {0}; LBDS feh(K,P) = {0}; STP feh(K,P) = {50}; UBDS feh(K,P) = {1000}; LBDS ffh(K,P) = {0}; STP ffh(K,P) = {0}; UBDS ffh(K,P) = {1000}; LBDS fgh(L,P) = {0}; STP fgh(L,P) = {50}; UBDS fgh(L,P) = {1000}; LBDS fh(P) = {0}; STP fh(P) = {100}; UBDS fh(P) = {1000}; LBDS ch(P,I) = {0}; STP ch(P,I) = {ca(R,I)}; UBDS ch(P,I) = {10}; LBDS vm(L) = {0}; STP vm(L) = {1}; UBDS vm(L) = {10000}; LBDS tm(L) = {900}; STP tm(L) = {900}; UBDS tm(L) = {1500}; LBDS rm(L,J) = {0}; STP rm(L,J) = {0}; UBDS rm(L,J) = {10000}; LBDS vt(K) = {0}; STP vt(K) = {1}; UBDS vt(K) = {1000}; LBDS ks(K,CP) = {0}; STP ks(K,CP) = {0}; UBDS ks(K,CP) = {1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,0}; LBDS ko(K,CP) = {0}; STP ko(K,CP) = {0}; UBDS ko(K,CP) = {0,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000}; LBDS kt(K,CP,OR) = {900}; STP kt(K,CP,OR) = {900}; UBDS kt(K,CP,OR) = {1500}; ZVAR {ft(K), #CFR flowrate ct(K,I), #CFR concentration tt(K), #CFR temperature co(K,I), #CFR side exit concentration fs(K), #CFR sidestream flowrate fo(K), #CFR side exit flowrate rt(K,J), #CFR reaction rates dtdv(K), dsdv(K), dodv(K) }; #States which specify the initial condition ISPE {ft(K), ct(K,I), co(K,I)}; #States whose initial conditions are to be optimized ICS {ft(K), ct(K,I)}; #Parameters corresponding to the initial conditions ICP {fc(K), cc(K,I)}; #Initial conditions of remaining z-varaibles STP co(K,I) = {0}; #10 intervals PARA time(CP) = {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}; TIME {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}; #collocation points PARA pt(CI,CO) = {0.00,0.05,0.10, 0.10,0.15,0.20, 0.20,0.25,0.30, 0.30,0.35,0.40, 0.40,0.45,0.50, 0.50,0.55,0.60, 0.60,0.65,0.70, 0.70,0.75,0.80, 0.80,0.85,0.90, 0.90,0.95,1.00}; }} MODEL {{ MIN: -ch(1,3); #feed splitter splita(r E R): fa(r) =e= <> + <> ; #pre-reaction mixers mixbt(k E K): fb(k) =e= <> + <> + <>; mixbc(k E K, i E I): cb(k,i)*fb(k) =e= <> + <> + <>; mixct(k E K): fc(k) =e= <> + <> + <>; mixcc(k E K, i E I): cc(k,i)*fc(k) =e= <> + <> + <>; mixdt(l E L): fd(l) =e= <> + <> + <>; mixdc(l E L, i E I): cd(l,i)*fd(l) =e= <> + <> + <>; #CSTRs cstrt(l E L): fg(l) =e= fd(l); cstrc(l E L, i E I): cg(l,i)*fg(l) =e= cd(l,i)*fd(l) + vm(l)*<>; #CSTR reaction rates rate1(l E L): rm(l,1) =e= k0(1)*exp[-e(1)/rg/tm(l)]*cg(l,1); rate2(l E L): rm(l,2) =e= k0(2)*exp[-e(2)/rg/tm(l)]*cg(l,1); rate3(l E L): rm(l,3) =e= k0(3)*exp[-e(3)/rg/tm(l)]*cg(l,2); rate4(l E L): rm(l,4) =e= k0(4)*exp[-e(4)/rg/tm(l)]*cg(l,3); #CFRs cfrt(k E K): ft'(k) =e= -dsdv(k) - dodv(k); cftc(k E K, i E I): ct'(k,i)*(ft(k)+1e-6) =e= vt(k)*<> - (cb(k,i)-ct(k,i))*dsdv(k); point1(k E K)[10]: fe(k) =e= ft(k); point2(k E K, i E I)[10]: ce(k,i)*fe(k) =e= ct(k,i)*ft(k); control3(k E K): tt(k) =e= interv[i E CI| kt(k,i ,0)*(t-pt(i,1))*(t-pt(i,2))/((pt(i,0)-pt(i,1))*(pt(i,0)-pt(i,2))) + kt(k,i ,1)*(t-pt(i,0))*(t-pt(i,2))/((pt(i,1)-pt(i,0))*(pt(i,1)-pt(i,2))) + kt(k,i+1,0)*(t-pt(i,0))*(t-pt(i,1))/((pt(i,2)-pt(i,0))*(pt(i,2)-pt(i,1))) ]; dcontrol3(k E K): dtdv(k) =e= interv[i E CI| kt(k,i ,0)*((t-pt(i,1))+(t-pt(i,2)))/(pt(i,0)-pt(i,1))/(pt(i,0)-pt(i,2)) + kt(k,i ,1)*((t-pt(i,0))+(t-pt(i,2)))/(pt(i,1)-pt(i,0))/(pt(i,1)-pt(i,2)) + kt(k,i+1,0)*((t-pt(i,0))+(t-pt(i,1)))/(pt(i,2)-pt(i,0))/(pt(i,2)-pt(i,1)) ]; #CFR reaction rates rate1(k E K): rt(k,1) =e= k0(1)*exp[-e(1)/rg/tt(k)]*ct(k,1); rate2(k E K): rt(k,2) =e= k0(2)*exp[-e(2)/rg/tt(k)]*ct(k,1); rate3(k E K): rt(k,3) =e= k0(3)*exp[-e(3)/rg/tt(k)]*ct(k,2); rate4(k E K): rt(k,4) =e= k0(4)*exp[-e(4)/rg/tt(k)]*ct(k,3); #CFR side feed mat1(k E K): fb(k) =e= ks(k,0); control1(k E K): fs(k) =e= interv[i E CI| ks(k,i )*(t-pt(i,2))/(pt(i,0)-pt(i,2)) + ks(k,i+1)*(t-pt(i,0))/(pt(i,2)-pt(i,0)) ]; dcontrol1(k E K): dsdv(k) =e= interv[i E CI| ks(k,i )/(pt(i,0)-pt(i,2)) + ks(k,i+1)/(pt(i,2)-pt(i,0)) ]; #CFR side exit exitc(k E K, i E I): co'(k,i)*(fo(k)+1e-6) =e= (ct(k,i)-co(k,i))*dodv(k); mat2(k E K): ff(k) =e= ko(k,10); point3(k E K, i E I)[10]: cf(k,i)*ff(k) =e= co(k,i)*fo(k); control2(k E K): fo(k) =e= interv[i E CI| ko(k,i )*(t-pt(i,2))/(pt(i,0)-pt(i,2)) + ko(k,i+1)*(t-pt(i,0))/(pt(i,2)-pt(i,0)) ]; dcontrol2(k E K): dodv(k) =e= interv[i E CI| ko(k,i )/(pt(i,0)-pt(i,2)) + ko(k,i+1)/(pt(i,2)-pt(i,0)) ]; #post-reactor splitters splite(k E K): fe(k) =e= <> + <> + <

>; splitf(k E K): ff(k) =e= <> + <> + <

>; splitg(l E L): fg(l) =e= <> + <> + <

>; #product mixer mixht(p E P): fh(p) =e= <> + <>; mixhc(p E P, i E I): ch(p,i)*fh(p) =e= <> + <>; #continuity of derivatives for temperature control cont1(k E K, i E CI & i>0): kt(k,i-1,0)*(pt(i-1,2)-pt(i-1,1))/ (pt(i-1,0)-pt(i-1,1))/(pt(i-1,0)-pt(i-1,2)) + kt(k,i-1,1)*(pt(i-1,2)-pt(i-1,0))/ (pt(i-1,1)-pt(i-1,0))/(pt(i-1,1)-pt(i-1,2)) + kt(k,i ,0)*(2*pt(i-1,2)-pt(i-1,1)-pt(i-1,0))/ (pt(i-1,2)-pt(i-1,0))/(pt(i-1,2)-pt(i-1,1)) - kt(k,i ,0)*(2*pt(i,0)-pt(i,2)-pt(i,1))/ (pt(i,0)-pt(i,1))/(pt(i,0)-pt(i,2)) - kt(k,i ,1)*(pt(i,0)-pt(i,2))/ (pt(i,1)-pt(i,0))/(pt(i,1)-pt(i,2)) - kt(k,i+1,0)*(pt(i,0)-pt(i,1))/ (pt(i,2)-pt(i,0))/(pt(i,2)-pt(i,1)) =e= 0; #logical constraints log1(k E K, i E CI): ks(k,i+1) =l= ks(k,i); log2(k E K, i E CI): ko(k,i+1) =g= ko(k,i); }}