Description
Description
Purpose
Purpose
This notebook computes the steady state of a Reactive Distillation (RD) column that decomposes methyl tert-butyl ether (MTBE) into iso-butene (IBUT) and methanol (MEOH), MTBE ⇄ IBUT + MEOH. It reproduces the 8(5)/6(4)/5 design of Huang & Wang (Ind. Eng. Chem. Res. 2007, 46, 2508) and plots the resulting profiles. Stage numbering follows the paper: stage 1 is the total condenser, stages 2-4 rectify, stages 5-15 are reactive, stages 16-20 strip, stage 21 is the partial reboiler, and pure MTBE is fed on stage 14.
Acknowledgement
Acknowledgement
AI was used interactively in order to code and solve this RD problem.
Model
Model
Each stage is an equilibrium (MESH) stage: material balances for IBUT and MEOH, a total mole balance, an energy balance, and the equilibrium condition that the vapour mole fractions sum to one. Reactive stages also satisfy chemical equilibrium with ln Keq = 16.33 - 6820/T. Vapour-liquid equilibrium uses Wilson activity coefficients and Antoine vapour pressures with the constants of the original notebook. The condenser is total, the reboiler partial, and pressure is a constant 11 atm.
The switch activityInKeq
The switch activityInKeq
In the first code block, activityInKeq = False applies Keq to mole fractions (Keq = (x_IBUT x_MEOH)/x_MTBE). activityInKeq = True puts the Wilson activity coefficients into the reaction equilibrium (Keq g_MTBE x_MTBE = g_IBUT g_MEOH x_IBUT x_MEOH). Using the paper's Fig. 6 compositions and temperatures back through the equilibrium relation shows that the paper's profile needs g_IBUT g_MEOH / g_MTBE close to 1, so activityInKeq = False reproduces the composition and temperature profiles much better. The analysis given at the end of the notebook support the assumption: activityInKeq = False.
Unknowns and specifications
Unknowns and specifications
The unknowns are the IBUT and MEOH liquid fractions, temperature, liquid flow and vapour flow on stages 2-21, the reaction extent on stages 5-15, the reflux ratio R, the reboiler duty Qreb and the drum temperature Td. There are 114 unknowns in total. There are 114 equations: 5 per stage (100), 11 chemical-equilibrium conditions, the distillate bubble-point condition, and two product specifications (94 mol% IBUT in the distillate, 94 mol% MEOH in the bottoms).
Energy balance
Energy balance
Enthalpies are measured from liquid at 298 K. A formation-enthalpy offset makes the decomposition absorb 37.7 kJ/mol. Vapour enthalpy is liquid enthalpy plus a latent heat that follows the Watson correlation from the 298 K values of the paper's Table 1. The liquid heat capacities (130, 80 and 190 J/mol/K for IBUT, MEOH and MTBE) are estimates so the duties are approximate.
Expected results
Expected results
For activityInKeq = False: R = 1.68, Qreb = 7.40 MW, Qcon = 4.74 MW, D = 100.24, B = 93.99, conversion 94.2%. Paper (Table 2): R = 2.87, Qreb = 7.39 MW, Qcon = 6.98 MW, D = 100.03, B = 94.00, conversion 94.03%. With False the compositions (RMS error about 0.03 against Fig. 6) and the reboiler duty match but the reflux and condenser duty are lower (the Wilson IBUT–MeOH equilibrium could be the likely cause). Using AI and UNIQUAC data (see Ref. 2) did not help either.
Known limitations
Known limitations
Activities are Wilson, not the UNIQUAC model of the paper and pressure is a constant 1114 kPa while the paper's stated 50 kPa drop per stage. Vapour is an ideal gas without mixing enthalpy. The model is steady state only.
Reference:
[1] K. Huang and S. J. Wang, Design and Control of a Methyl Tertiary Butyl Ether (MTBE) Decomposition Reactive Distillation Column, Ind. Eng. Chem. Res., 46(8), 2007, pp. 2508-2519.
[2] A. Rehfinger and U. Hoffmann, Kinetics of methyl tertiary butyl ether liquid phase synthesis catalyzed by ion exchange resin—I. Intrinsic rate expression in liquid phase activities. Chemical Engineering Science 1990, 45 (6), 1605–1617.
[1] K. Huang and S. J. Wang, Design and Control of a Methyl Tertiary Butyl Ether (MTBE) Decomposition Reactive Distillation Column, Ind. Eng. Chem. Res., 46(8), 2007, pp. 2508-2519.
[2] A. Rehfinger and U. Hoffmann, Kinetics of methyl tertiary butyl ether liquid phase synthesis catalyzed by ion exchange resin—I. Intrinsic rate expression in liquid phase activities. Chemical Engineering Science 1990, 45 (6), 1605–1617.
Code
Code
Parameters, switch and thermodynamic functions
Parameters, switch and thermodynamic functions
In[]:=
activityInKeq=False;(*SWITCH:False=Keqappliedtomolefractions(bestmatchtopaperFig6compositions);True=WilsonactivitycoefficientsinsideKeq*)
In[]:=
aij={{0,169.9953,-30.2477},{2576.8532,0,1483.2478},{271.5669,-406.3902,0}};(*Wilsonenergyparametersa_ijincal/mol;row/columnorder=IBUT,MEOH,MTBE*)Vm={93.33,44.44,118.8};(*liquidmolarvolumesincm3/mol(IBUT,MEOH,MTBE),usedintheWilsonmodel*)
In[]:=
Rg=1.987;(*gasconstantincal/mol/K,matchestheunitsofaij*)
In[]:=
P=11*1.013*^5;(*columnpressure:11atmexpressedinPa,takenasconstantoneverystage*)
In[]:=
Aa={20.64559,23.49989,20.71616};(*AntoineconstantA(lnofvapourpressureinPa)*)Ba={-2125.74886,-3643.31362,-2571.58460};(*AntoineconstantB*)Ca={-33.160,-33.434,-48.406};(*AntoineconstantC(K)*)
In[]:=
Tc={417.9,512.6,497.1};(*criticaltemperaturesinK,usedforthelatent-heattemperaturecorrection*)
In[]:=
lam298={20270.,38020.,29730.};(*latentheatsofvaporisationat298KinJ/mol(paperTable1)*)
In[]:=
cpL={130.,80.,190.};(*liquidheatcapacitiesinJ/mol/K--ESTIMATES,notfromthepaper*)
In[]:=
dHf={0.,37700.,0.};(*formation-enthalpyoffsetsinJ/molsothatMTBE->IBUT+MEOHabsorbs+37.7kJ/mol*)
In[]:=
F=100.;zF={0,0,1};TF=426.73;(*feed:100mol/sofpureMTBE,saturatedliquidat426.73K*)
In[]:=
jf=14;reac=Range[5,15];(*feedstagenumberandlistofreactivestages(papernumbering)*)
In[]:=
gam[x_List,T_]:=Module[{A=Table[Vm[[j]]/Vm[[i]]Exp[-aij[[i,j]]/(RgT)],{i,3},{j,3}],s},(*Wilsonactivitycoefficients:buildmatrixA_ij=Vm_j/Vm_iexp(-a_ij/RT)*)s=A.x;(*s_i=sum_jA_ijx_j*)Exp[1-Log[s]-Table[Sum[x[[k]]A[[k,i]]/s[[k]],{k,3}],{i,3}]]];(*ln(gamma_i)=1-lns_i-sum_kx_kA_ki/s_k*)
In[]:=
psat[T_]:=Exp[Aa+Ba/(T+Ca)];(*AntoinevapourpressuresofthethreespeciesattemperatureT,inPa*)
In[]:=
Keq[T_]:=Exp[16.33-6820/T];(*reactionequilibriumconstant,papereq4*)
In[]:=
hL[x_List,T_]:=x.(dHf+cpL(T-298.));(*liquidmolarenthalpyofamixture(J/mol)relativetoliquidat298K*)
In[]:=
lam[T_]:=lam298((1-T/Tc)/(1-298./Tc))^0.38;(*Watsoncorrelation:latentheatatTfromthe298Kvalue*)
In[]:=
HV[y_List,T_]:=y.(dHf+cpL(T-298.)+lam[T]);(*vapourmolarenthalpy=liquidenthalpy+latentheat(idealgas,nomixingeffects)*)
Stage equations
Stage equations
In[]:=
xv[j_]:={x1[j],x2[j],1-x1[j]-x2[j]};(*liquidcompositionofstagejas{IBUT,MEOH,MTBE};MTBEbydifference*)yv[j_]:=xv[j]*gam[xv[j],T[j]]*psat[T[j]]/P;(*equilibriumvapourcompositiononstagej:y_i=x_igamma_iPsat_i/P*)ydv=yv[2]/Total[yv[2]];(*totalcondenser:distillateandrefluxhavethecompositionofthevapourleavingstage2*)Dflow=V[2]/(1+R);L1=RDflow;(*distillateflowfromV2=D(1+R),andrefluxflowL1=RD*)Lin[j_]:=If[j==2,L1,L[j-1]];(*liquidenteringstagej:refluxforstage2,otherwisetheliquidfromthestageabove*)xin[j_]:=If[j==2,ydv,xv[j-1]];(*compositionofthatenteringliquid*)hin[j_]:=If[j==2,hL[ydv,Td],hL[xv[j-1],T[j-1]]];(*enthalpyofthatenteringliquid(refluxsitsatthedrumtemperatureTd)*)Vin[j_]:=If[j==21,0,V[j+1]];(*vapourenteringstagejfromthestagebelow;nothingentersthereboiler(stage21)*)yin[j_]:=If[j==21,{0,0,0},yv[j+1]];(*compositionofthatenteringvapour*)Hin[j_]:=If[j==21,0,HV[yv[j+1],T[j+1]]];(*enthalpyofthatenteringvapour*)Fs[j_]:=If[j==jf,F,0];(*feedflow,nonzeroonlyonthefeedstage*)xiF[j_]:=If[5<=j<=15,xi[j],0];(*reactionextent(mol/sofMTBEdecomposed),zeroonnon-reactivestages*)stageEqs[j_]:=Module[{en},(*thefiveequationswrittenforstagej*)en=Lin[j]hin[j]+Vin[j]Hin[j]+Fs[j]hL[zF,TF](*energyin:enteringliquid+enteringvapour+feed*)-L[j]hL[xv[j],T[j]]-V[j]HV[yv[j],T[j]](*energyout:leavingliquidandvapour*)+If[j==21,Qreb,0];(*reboilerheatdutyaddedonthelaststage(21)only*){(Lin[j]xin[j][[1]]+Vin[j]yin[j][[1]]+Fs[j]zF[[1]]+xiF[j]-L[j]xv[j][[1]]-V[j]yv[j][[1]])/F==0,(*IBUTbalance:in+generated-out,scaledbyF*)(Lin[j]xin[j][[2]]+Vin[j]yin[j][[2]]+Fs[j]zF[[2]]+xiF[j]-L[j]xv[j][[2]]-V[j]yv[j][[2]])/F==0,(*MEOHbalance(MEOHisalsogenerated1:1)*)(Lin[j]+Vin[j]+Fs[j]+xiF[j]-L[j]-V[j])/F==0,(*totalmolebalance(reactionaddsonemolepermolereacted)*)en/(F3.*^4)==0,(*energybalance,scaledtobeO(1)*)Total[yv[j]]-1==0}];(*equilibriumclosure:vapourmolefractionssumto1,whichfixesthestagetemperature*)gg[j_]:=If[activityInKeq,gam[xv[j],T[j]],{1,1,1}];(*activitycoefficientsusedinthereactionequilibrium;all1whentheswitchisFalse*)reacEq[j_]:=Keq[T[j]]gg[j][[3]]xv[j][[3]]-gg[j][[1]]gg[j][[2]]x1[j]x2[j]==0;(*chemicalequilibriumonareactivestage:Keqg_MTBEx_MTBE=g_IBUTg_MEOHx_IBUTx_MEOH*)bubbleEq=Total[ydv*gam[ydv,Td]*psat[Td]/P]-1==0;(*drumtemperatureTdisthebubblepointofthedistillate*)specs={ydv[[1]]==0.94,x2[21]==0.94};(*productspecs:94mol%IBUTdistillate,94mol%MEOHbottoms*)eqs=Flatten[{Table[stageEqs[j],{j,2,21}],reacEq/@reac,bubbleEq,specs}];(*collectall114equationsintoonelist*)
Unknowns and starting values
Unknowns and starting values
In[]:=
startA={{0.9686260378,0.9669481816,0.9270819782,0.8138762959,0.6539000384,0.5031767584,0.3991679221,0.3324765312,0.2868100222,0.2524316873,0.2242278816,0.199418591,0.1764298977,0.0795602532,0.01179697636,0.001542877628,0.0001853660704,2.00717385*^-05,1.916158283*^-06,1.581259821*^-07},{0.02489098254,0.01172607685,0.009557864783,0.01196412076,0.03520600356,0.08085407303,0.1383481247,0.1939669001,0.2441725864,0.290344405,0.3347469337,0.3795263592,0.4265601456,0.7044619881,0.7618379747,0.7883666929,0.8186426666,0.8573341421,0.9002941225,0.94},{347.59192,348.8170679,350.6520742,354.9493925,360.2169237,365.589966,369.543307,372.194153,374.0737333,375.5357714,376.7785195,377.9191332,379.0354568,385.7688251,403.682569,407.4860151,408.4252738,409.19481,410.2147587,411.420275},{173.6665535,175.0043874,170.8448916,141.8931202,116.3657873,102.2838011,97.51174694,96.81413207,97.62532378,99.04632267,100.784821,102.7526258,253.6096794,357.2613657,365.0043034,365.8946967,364.0943525,361.5673128,359.087201,93.98534431},{268.4695363,273.9110626,275.2488966,271.0894007,247.9996748,227.2000371,215.2568558,210.3487976,208.5633084,207.9930735,207.9331214,208.1194707,208.4419277,219.4933921,263.2760214,271.0189591,271.9093524,270.1090082,267.5819685,265.1018567},{-5.862045423,-4.72769521,-2.138804924,0.1360040231,1.087874361,1.381426668,1.480950955,1.552148973,1.645347891,39.80558915,59.86905701},1.678147,7402747.3,346.692};(*startvalues{x1,x2,T,L,V,xi,R,Qreb,Td}foractivityInKeq=False:convergedPythonsolution,stages2-21*){x1s,x2s,Ts,Ls,Vs,xis,R0,Q0,Td0}=startA;(*pickthestartsetthatbelongstothechosenswitchvalue*)vars=Join[(*listof{unknown,startingvalue}pairshandedtoFindRoot*)Table[{x1[j],x1s[[j-1]]},{j,2,21}],(*IBUTfractions*)Table[{x2[j],x2s[[j-1]]},{j,2,21}],(*MEOHfractions*)Table[{T[j],Ts[[j-1]]},{j,2,21}],(*temperatures*)Table[{L[j],Ls[[j-1]]},{j,2,21}],(*liquidflows*)Table[{V[j],Vs[[j-1]]},{j,2,21}],(*vapourflows*)Table[{xi[reac[[k]]],xis[[k]]},{k,Length[reac]}],(*reactionextents:k-threactivestagepairedwithk-thstartvalue*){{R,R0},{Qreb,Q0},{Td,Td0}}];(*refluxratio,reboilerdutyinW,drumtemperatureinK*)
Solve
Solve
Results
Results
Profile plots
Profile plots
Composition comparison with Huang & Wang Fig. 6
Composition comparison with Huang & Wang Fig. 6
Temperature comparison with Huang & Wang Fig. 6
Temperature comparison with Huang & Wang Fig. 6
Flow comparison with Huang & Wang Fig. 6
Flow comparison with Huang & Wang Fig. 6
Here, we observe the same discrepancy as the reflux ratio (1.68 in the model against the paper' s 2.87) . Columns above the feed run at a higher internal reflux in the paper than in the Mathematica model . Together with the top temperature being 4 K high, it points to the Wilson IBUT–MeOH equilibrium as the likely cause.