Material and Energy Balances in a Reactor with Heat Exchange

Manipulate[
Module[{Shomate,H,Tin,Fin,CpAcetylene, CpEthylene, ν,Fout,Hin,Hout,sol,tout,ext,result,maxExtent,reactor,components},

(*Moved the reactor from initialization to here so that everything is contained within a single evaluation cell. Personal preference, doesn't make a difference really.*)
reactor=Graphics[{
{EdgeForm@Thick,FaceForm@RGBColor[0.6,1,1],Rectangle[{1,-1/3},{2,1/3}]},
{Thick,Arrow[{{0,0},{1,0}}],Arrow[{{2,0},{3,0}}],Line[{{1.25,-0.5},{1.25,-0.15}}],Arrow[{{1.75,-0.15},{1.75,-0.55}}],
Line[{{1.25+#*0.05,-0.15},{1.25+(1+#)*0.05,-0.05},{1.25+(2+#)*0.05,-0.15}}]&/@Range[0,8,2]},
Text[Style["A + B ⟶ C",22],{1.5,0},{0,-1}]
},ImageSize->{555,350},PlotRange->{All,0.7}];

Tin=25;(*Inlet temperature, degrees C*)
Fin={100,nB0,0,nI}; (*Inlet flow rates of { acetylene, hydrogen, ethylene, and nitrogen } *)
maxExtent=Fin[[2]];(*maximum extent of reaction is necessary to calculate the maximum temperature possible*)

If[extent>maxExtent,extent=maxExtent];(*Cannot have extent of reaction greater than the mols of the limiting reactant.*)

(*
Shomate is a list of parameters of the Shomate equation. In the previous version of this file, this was a function rather than a list. I converted it to a list because it consists purely of constants, so it doesn't make sense as a function.
*)
Shomate[T_]:={
If[
T<1100,(*Acetylene*)
{40.68697,40.73279,-16.17840,3.669741,-0.658411,210.7067,226.7314},
{67.47244,11.75110,-2.021470,0.136195,-9.806418,185.4550,253.5337,226.7314}
],
Which[
T<1000,(*Hydrogen*)
{33.066178,-11.363417,11.432816,-2.772874,-0.158558,-9.980797, 0},
T<2500,
{18.563083,12.257357,-2.859786,0.268238,1.977990,-1.147438,156.288133,0.0},
T>=2500,
{43.413560,-4.293079,1.272428,-0.096876,-20.533862,-38.515158,162.081354,0.0}
],
If[
T<1200,(*Ethylene*)
{-6.387880,184.4019,-112.9718,28.49593,0.315540,48.17332,52.46694},
{106.5104,13.73260,-2.628481,0.174595,-26.14469,-35.36237,275.0424,52.46694}
],
Which[
T<500,(*Nitrogen*)
{28.98641,1.853978,-9.647459,16.63537,0.000117,-8.671914,0},
T<2000,
{19.50583,19.88705,-8.598535,1.369784,0.527601,-4.935202,212.3900,0.0},
T>=2000,
{35.51872,1.128728,-0.196103,0.014662,-4.553760,-18.97091,224.9810,0.0}
]
};

(*stoichiometric coefficients -{acetylene, hydrogen, ethylene, nitrogen}. Also converted to a list, was previously a function*)
ν={-1,-1,1,0};


(*Enthalpy from Shomate equation, using paramaters taken from NIST.gov. You'll notice that I made this equation much more concise.*)
H[n_][T_?NumberQ]:=With[{t=(((T+273.15)/1000))},(Sum[((Shomate[T+273.15][[n]][[i]]*(t^i))/i),{i,1,4,1}]-((Shomate[T+273.15][[n]][[5]])/t)+Shomate[T+273.15][[n]][[6]])*1000];

(*Inlet enthalpy*)

Hin[n_]:=H[n][Tin];

(*Outlet molar flow rate, mass balance*)

Fout[n_][ζ_]:=Fin[[n]]+ν[[n]]*ζ;

(*"sol" is solution to the extent of reaction if the user selects "outlet temperature", and sol is the solution to outlet temperature should the user select "extent of reaction". This is a simple material and energy balance, where the enthalpy of the components entering the reactor must equal the enthalpy of components leaving the reactor plus heat removed (Q)*)
sol=Quiet@If[npr==1,
result=ζ/.FindRoot[(∑_(i=1)^4)(Fin[[i]]*Hin[i])==(∑_(i=1)^4)(Fout[i][ζ]*H[i][T])+Q*1000,{ζ,Fin[[2]],0,Fin[[2]]}][[1]];
If[result>Fin[[2]],result=Fin[[2]]];
result,
Tout/.FindRoot[(∑_(i=1)^4)(Fin[[i]]*Hin[i])==(∑_(i=1)^4)(Fout[i][extent]*H[i][Tout])+Q*1000,{Tout,Tin,Tin,5000}][[1]]
];


tout=If[npr==1,T,sol];
ext=If[npr==1,sol,extent];

(*The text next to the Q slider should have maximum two digits after the decimal place, and the maximum Q slider value should be equal to or less than the maximum amount of heat that can be removed, so we use a floor function in this if statement to remove 3 or more digits after the decimal place while ensuring you can never remove more heat than maxHeat. We could use Round[ ] instead of Floor[ ], but that would sometimes result in Q being higher than maxHeat, which we don't want. Floor guarantees we only set it to a value equal or lower than the maximum possible value*)
If[npr==1,
maxHeat=- ((∑_(i=1)^4)(Fout[i][nB0]*H[i][tout])-(∑_(i=1)^4)(Fin[[i]]*Hin[i]))/1000,
maxHeat =- ( (∑_(i=1)^4)(Fout[i][extent]*H[i][Tin]) -(∑_(i=1)^4)(Fin[[i]]*Hin[i]))/1000 
];
(*If maxHeat is above 100, we only need 0 digits after the decimal place, not 1. If max heat is above 10 but below 100, put only 1 digit after the decimal place. We set the value of precision so that it only runs this code once.*)
Which[
maxHeat>=100&&maxHeat<=Q,maxHeat=Floor[maxHeat],
maxHeat>=10&&maxHeat<=Q,maxHeat=Floor[maxHeat,0.1],
maxHeat>=0&&maxHeat<=Q,maxHeat=Floor[maxHeat,0.01]
];

(*If the Q slider is higher than the maximum amount of heat that can be removed, then set it to the maximum amount of heat that can be removed (rounded down).
Q!=0 is included so that this statement is not evaluated again after Q has already been set to zero*)
If[Q>maxHeat&&Q!=0,Q=Floor[maxHeat]];
If[maxHeat<=0,
Q=0;(*If heat needs to be added to maintain a given temperature, make sure Q is not a negative value, as we cannot add heat in this simulation*)
];
If[npr==1,
(*Since Q=0 in the calculation of maximum temperature possible, we omitted Q from the equation below. It is otherwise the same energy balance as the FindRoot statement above.
This calculates the temperature achieved if extent of reaction is its maximum possible value, corresponding to the maximum possible temperature increase. If T is higher than that temperature (maxTemp), then the maximum value of the temperature slider is set to to maxTemp. Floor function is used for significant figures again*)
maxTemp=Quiet[Tout/.FindRoot[(∑_(i=1)^4)(Fin[[i]]*Hin[i])==(∑_(i=1)^4)(Fout[i][maxExtent]*H[i][Tout])+Q*1000,{Tout,Tin,Tin,5000}][[1]]];
];
(*If the maximum temperature value is somehow higher than what is possible, set T to the maximum possible temperature.*)
If[T>maxTemp,T=Floor[maxTemp]];

components={"A","B","C","I"};

Show[reactor,Graphics@{
Text[Style["Sliders may adjust their ranges to prevent impossible scenarios.",18,Red],{1.55,0.65},{0,0}],
Text[Style[Row@{NumberForm[Tin,{4,0}]," °C"},20,Black],{0.5,0},{0,-3}],
Text[Style[Row@{NumberForm[tout,{4,0}]," °C"},20,If[tout≤0∨tout≥1000,RGBColor[0.8,0,0],Black]],{2.5,0},{0,-3}],
Text[Grid[
(*I changed the three nearly identical statements into a table. Making the code concise.*)
Table[{Style[NumberForm[Fin[[n]],{5,1}],Black,20],Style[Row@{"mol ",components[[n]],"/h"},Black,18]},{n,{1,2,4}}]
],{0.5,-0.25}],
Text[Grid[
(*Same thing here. Trimming the fat off of the code.*)
Table[{Style[NumberForm[Round[Fout[n][ext],0.1],{5,1}],Black,20],Style[Row[{"mol ",components[[n]],"/h"}],Black,18]},{n,1,4}]
],{2.5,-0.3}],
Text[Style[Column[{
Row[{NumberForm[Q,{10,0}]," kW heat removed" }],
Row[{"conversion = ",NumberForm[N@100*((Fin[[1]]-Fout[1][ext])/(Fin[[1]])),{4,1}],"% of A"}]
},Center,0.75],20],{1.5,-0.5},{0,1.65}]
}]
],

Grid[{
{Control[{{npr,1,"change conversion with:"},{1->" outlet temperature ",2->" extent of reaction "},Setter}]},
{PaneSelector[{
1->Control[{{T,1000,"outlet temperature (°C)"},25,maxTemp,1,Appearance->"Labeled"}],
2->Control[{{extent,50,"extent of reaction (mol/h)"},0,nB0,1,Appearance->"Labeled"}]
},Dynamic@npr],SpanFromLeft},
{Control[{{Q,1000,"heat removal rate (kW)"},0,maxHeat,1}],Dynamic@NumberForm[Q,{10,0}]},
{Control[{{nB0,75,Row@{"feed rate ",Subscript["H",2]," (mol B/h)"}},10,100,1,Appearance->"Labeled",ImageSize->Small}],
Control[{{nI,50,"feed rate inert (mol/h)"},0,100,1,Appearance->"Labeled",ImageSize->Small}],

(*We need dynamic variables that can be referenced outside of the Module[] function, because we need to dynamically update the slider limits in order to prevent impossible scenarios from happening. These "controls" are not visible and just act as dynamic variables.*)
Control[{{maxHeat,10000},0,1000,None}],
Control[{{maxTemp,79},10,300,None}]
}
},Alignment->Left],
SaveDefinitions->True]

This Demonstration calculates material and energy balances for a reactor with heat transfer. The reaction A+B->C is acetylene (A) hydrogenation to ethylene (C), where B is hydrogen. The feed rate of A is fixed at 100 mol/h; use sliders to change the feed rates of B and an inert (I). Use either the extent of reaction or the outlet temperature to change the conversion of A. Select "extent of reaction" to set the extent of reaction and heat removal rate with sliders to change the outlet temperature. Select "outlet temperature" to change the outlet temperature and heat removal rate to determine the conversion, and the heat removed is fixed to the value that results in the selected outlet temperature. If the extent of reaction is restricted by the limiting reactant B, the process indicates that the limiting conversion is reached by displaying the message "limiting reactant depleted". If the outlet temperature is below 0 °C or above 1,000 °C, the message "operating conditions unrealistic" appears.

The reaction A+B->C takes place in a reactor with heat exchange. Reactants A and B and inert I enter the reactor. The molar flow rates of each component exiting the reactor are calculated using the extent of reaction:

(Subscript[F, i,out])=(Subscript[F, i,in])+(Subscript[ν, i]) ζ,

where (Subscript[F, i,out]) and (Subscript[F, i,in]) are the molar flow rates (mol/h) of component i=A,B,C,I exiting and entering the reactor, (Subscript[ν, i]) is the stoichiometric coefficient, which is positive for products and negative for reactants, and ζ is the extent of reaction (mol/h).

The extent of reaction is:

ζ=-(((∑_i)(Subscript[F, i,in]) (Subscript[C, p,i])((Subscript[T, out])-(Subscript[T, in]))+Q)/((∑_i)(Subscript[ν, i]) (Subscript[ΔH, f,i]))),

where Δ(Subscript[H, f,i]^o) is the heat of formation (kJ/mol) of component i with a reference state at T=25°C, (Subscript[Cp, i]) is the constant heat capacity (kJ/[mol °C])of component i, (Subscript[T, in]) is the feed temperature (°C), (Subscript[T, out]) is the outlet temperature (°C) and Q is the heat added (kJ/h). The heat removed is -Q.

The outlet temperature is:

(Subscript[T, out])=(Subscript[T, in])-(((∑_i)((Subscript[ν, i]) (Subscript[ΔH, f,i]))*ζ-Q)/((∑_i)(Subscript[F, i,out]) (Subscript[Cp, i]))).

The percent conversion of A is

% conversion = (((Subscript[F, A,in]) -(Subscript[F, A,out]))/(Subscript[F, A,in])) * 100.

Note that constant heat capacities, (Subscript[Cp, i]), were assumed for simplicity. For large temperature changes, this is not a good assumption.

This screencast video [1] show calculations for an energy balance on a reactor with heat transfer. This screencast [2] explores limiting reagents.

References

[1] Energy Balances with Unknown Outlet Conditions [Video]. (Dec 6, 2012) www.youtube.com/watch?v=FPDPxAp-Ms0.

[2] Limiting Reagent (Interactive) [Video]. (Nov 12, 2014) www.youtube.com/watch?v=7H2EnOzGxD0&feature=youtu.be&list=PLD4476BAFA5A65111.

Contributed by: Michael Wrobel and Rachael L. Baumann

Additional contributions by: Garret D. Nicodemus, John L. Falconer, Neil Hendren, and Nick Bongiardina

(University of Colorado Boulder, Department of Chemical and Biological Engineering)