Optimization of a reactor-separator train

This notebook has been written in Mathematica by 

Mark J. McCready
Professor of Chemical and Biomolecular Engineering
University of Notre Dame
Notre Dame IN 46556
USA

Mark.J.McCready.1@nd.edu
ButtonBox[RowBox[{RowBox[{"http",":"}],"//",RowBox[{"chemeprof",".","com"}]}],BaseStyle->"Hyperlink",ButtonData->{"chemeprof.com",None}]


It is copyrighted to the extent allowed by whatever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy, that this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  6/3/21

Problem statement

We want to find an analytical expression(s) describing the costs of a reactor followed by a distillation column.  From this we will see how physical properties affect the outcome and how decisions about the residence time of the reactor, affect the size and energy costs of a distillation column.  
We will use a A--> B reaction with first order kinetics

Solution

We will consider the reactor to be a ideal CFSTR.  We will  feed pure A and there will be a recycle stream that is (essentially) pure A.  

Let us first solve the problem we know well.

First order kinetics, r = k (Subscript[C, A])

The mass balance for component A

eqa1 =  q (Subscript[C, A0 ]) - q (Subscript[C, A ])- k (Subscript[C, A ])V

eqa2= ExpandAll[eqa1/((Subscript[C, A ]) q)]

We introduce the conversion, X, (fraction of feed that has reacted) and the residence time, τ.

eqa3=eqa2/.{V->τ q, (Subscript[C, A ]) -> (Subscript[C, A0 ])(1-X)}

Here is an expression for the conversion.

ans1=Solve[eqa3==0,X]

For simplicity, I introduce ξ = k τ so that we can vary the dimensionless residence time for comparison.

ans2=ans1/.{k-> 1/τ ξ}

Recall the behavior of the conversion function.

Plot[Evaluate[X/.ans2[[1]]],{ξ,0,10},PlotRange->{0,1}]

When we solve a distillation problem, the minimum information about the chemistry that we need is the equilibrium vapor and liquid compositions.  This information is commonly given in the form of y-x diagrams.  These are plots of vapor composition, y versus liquid composition, x.

Recall relation for VLE from CBE 30367  (or last week).  

	(Subscript[x, i]) (Subscript[γ, i]) (Subscript[P, i]^vap) = (Subscript[y, i]) P

For mixtures of nonpolar, chemically similar molecules of close molecular weights (i.e., hydrocarbons), the mixtures are often "ideal".  This means that no chemical or size effects influence the compositions.  For these systems, the activity coefficient may be very close to 1.  Thus we have: 
		   (Subscript[y, i])   = (Subscript[x, i])   (Subscript[P, i]^vap)/ P.

The composition is only a function of temperature.  (If it were not we might write: y = K x and the "K" is found in table or nomograph.)   

For an ideal mixture we might find that the ratio

		α ≡  (( (Subscript[y, 1])  /(Subscript[x, 1]))/((Subscript[y, 2])  /(Subscript[x, 2]))) = (( (Subscript[P, 1]^vap)  /P)/((Subscript[P, 2]^vap)  /P))= (( (Subscript[P, 1]^vap) )/((Subscript[P, 2]^vap) ))
		
 is almost constant.  This is called "constant relative volatility" and α is the relative volatility.

Since we have a binary mixture, (Subscript[y, 2]) = 1 - (Subscript[y, 1]) (etc.)  thus:

alfaeq=α - ((y/(1-y))/(x/(1-x)))

alfans=Solve[alfaeq==0,y]

whyeq=y/.alfans[[1]]

Let us see what this looks like:

alfaplot=Plot[{x,whyeq/.α->3},{x,0,1},AspectRatio->1,Frame->True,PlotRange->{{0,1},{0,1}}]

The 45 degree line is commonly plotted for reference, for further construction, and undoubtedly for visual appeal!

It turns out that we will want the intersection of a line with the equilibrium curve at x = (Subscript[x, F]), the feed composition to the distillation column.

yfeed = whyeq/.{ x->(Subscript[x, F])}

I need to make a line out of two points, so I define this function:

lineqfun[ x0_,y0_,x1_,y1_ ]:=   ((y1-y0)/(x1-x0))(x-x0) + y0;

The two points that I need a line through are x=y=(Subscript[x, D]),and which is the final product composition. and (Subscript[x, F]),yfeedas defined above.  This assumes a saturated liquid feed.

refluxline = lineqfun[ (Subscript[x, D]),(Subscript[x, D]) ,(Subscript[x, F]),yfeed]

refplot=Plot[refluxline/.{(Subscript[x, D])->0.9`,(Subscript[x, F])->0.4`,α->3},{x,0,.9},AspectRatio->1,
PlotRange->{{0,1},{0,1}},PlotStyle->Red]

Here is the construction, part of the McCabe  - Thiele method for distillation design.

Show[refplot,alfaplot]

Now we can find a quantity called the reflux ratio, L/D which the ratio of liquid sent back down the column to the product distillate that is recycled back down the tower (so that it will work!).  We find the information that we need at the intercept.

temp1 =((Subscript[x, D])/(Rmin+1)) == refluxline/.{(Subscript[x, D])->xd,(Subscript[x, F])->xf,x->0}/.{xf->(Subscript[x, F]),xd->(Subscript[x, D]) }

temp2 = Solve[temp1,Rmin]

rmin=Rmin/.temp2[[1]]

We have calculated the smallest possible value for the reflux ratio, L/D.  Note that this is a form of the Underwood equation p 447, King (Separation Processes) 1st edition -- which is for L/F.  Both assume a saturated liquid feed.  Note that even if very high purity is desired, the minimum reflux ratio is still bounded as 2/(α-1).  For difficult separations increasing the reflux ratio will not produce the rerquired concentrations.  The trays must be increased.

The vapor flowrate in the column, which is constant and equal to the amount of vapor that we need to produce in the reboiler (and hence pay energy costs for) is G (aka, gee).  From a mass balance around the top of the column we would find:

gee = D (Rmin*β+1)/.temp2[[1]]

We have used Rmin with a factor β, to turn the minimum reflux into the value we use for reflux.

Since G is constant in the column, the costs to maintain a vaporization rate equal to G is

energycosts = gee (Subscript[C, E]) (Subscript[H, v])

where we have used (Subscript[H, v])for heat of vaporization and (Subscript[C, E]) for the cost per unit of energy.

energycosts/.{(Subscript[x, D])->.99, D->1000,(Subscript[C, E]) -> .01,(Subscript[H, v])->1000,α->3,β->1.5}

Let's make a plot to see how the energy costs change as we change the feed composition to the column.

We are making the same amount of distillate, D, and the composition is the same, (Subscript[x, D]), however if the feed is more dilute, we need to boil a lot more liquid.

Plot[energycosts/.{(Subscript[x, D])->0.95`,D->1000,(Subscript[C, E])->0.01`,(Subscript[H, v])->1000,β->1.2`,α->3},{(Subscript[x, F]),0.1`,0.9`},PlotRange->{0,70000}]

Translation of notation for the reactor into the distillation column

When we analyze the reactor, we have volumetric flowrates and molar compositions.  When we do a distillation column, we need molar flowrates, which is fortunately the product, q (Subscript[C, A]).  We then use mole fractions for compositions.

What is a value of (Subscript[C, A]) that might occur?  If the density of the liquid is 0.8 g/ml and the molecular weight is 75 g/gmol, we would have

N[800 g/l/(75 g/gmol)]

I can calculate the mole fraction of the feed stream to the distillation column: (and so can you!) as moles/time of B, divided by total mole flow rate.

exf =(( (Subscript[C, A0 ]) X)/(X (Subscript[C, A0 ]) + (1-X) (Subscript[C, A0 ])))

This is (Subscript[x, F]).  And it simplifies (for this special case) to:

Simplify[exf]

We can now write the energy costs in terms of reactor variables and output variables

energycosts

dist1energy = energycosts/.{(Subscript[x, F])->X}

dist2energy =dist1energy/.ans2[[1]]

Let’s check energy needed for separation versus the amount of time for the reaction.

Plot[dist2energy/.{(Subscript[x, D])->0.95`,D->1000,(Subscript[C, E])->0.01`,(Subscript[H, v])->1000,β->1.2`,α->3},{ξ,0.5`,10},AxesLabel->{"reaction time","energy cost"},PlotRange->{0,27000}]

The reactor cost will depend on reactor volume.  The tank is pretty simple, and the relation I found was

reactorcost = (Subscript[C, R]) (V^(.75))

The factor (Subscript[C, R])is installed cost per volume of reactor.

To get the reactor volume, we need the real flow rate, q.  This stream directly goes into the distillation column as F.  Thus we need to determine how much bigger F is than D.  Since the mole fraction in the feed is X, and X F is essentially D, then F = D/X.   

Now we need to make sure of the conversion between F and q.  F is moles/time.  q is liters/time.  We decided that no change of moles occurs, no change in density occurs so (Subscript[C, a0]) q = moles/time = F.  So that q= D/X/(Subscript[C, a0])

If we look back at ξ, which is kτ and recall that τ = V/q.

We can find a way to substitute for V in terms of ξ and

temp1= k V/q  - ξ

temp2=temp1/.q->D/ X/Ca0

temp3=temp2/.X->ξ/(1+ξ)

temp4 =Solve[%==0,V]

reactorcost = (Subscript[C, R]) (V^(.75))/.temp4[[1]]

We need the cost of the distillation column.  It depends on the flowrate and the number of trays.  We have the flowrate dependence from a chart in Peters and Timmerhaus.  However, we need some sort of relation between the number of trays and the reflux ratio, to finally complete this.

distcolumncost =( (Subscript[C, Area]) (F^(.8))) numbertrays

distcolumncost=distcolumncost/.F->D/X

distcolumncost=distcolumncost/.X->ξ/(1+ξ)

We need to determine how the number of trays varies.  The Fenske equation gives minimum trays -- corresponding to total reflux.

Fenske = Log[ ( xd/(1-xd) (1-xb)/xb)]/Log[α]


We expect the number of trays to decrease as the reflux ratio increases.  Here I pick a function.  (I could do better with some scaling theory!)

numbtray1=Log[((Subscript[x, D])/(1-(Subscript[x, D])) (1-(Subscript[x, B]))/(Subscript[x, B]))]/Log[α](Rmin/(((R-Rmin))^(.5)))

numbtray2 =numbtray1/.{Rmin->rmin,R->rmin β}

We need to get rid of xf

numbtray3 =numbtray2/.(Subscript[x, F])->X

Get this in terms of ξ.

numbtray4=numbtray3/.X->ξ/(1+ξ)

Finally do the substitutions

distcolumncost1=distcolumncost/.numbertrays->numbtray4

Equations for the total cost for this scheme

Totalcostfunction = distcolumncost1+reactorcost+dist2energy

reactorcostratio = reactorcost/(distcolumncost1)

distcostratio =  dist2energy/  distcolumncost1

capcosts =   distcolumncost1+reactorcost

distcosts = distcolumncost1+dist2energy

Plots of the cost functions

Let’s make some plots!

β is fractional increase of the reflux ratio above the minimum.

numbs2={(Subscript[x, D])->.95,(Subscript[x, B])->.05, D->1000,(Subscript[C, E]) -> .03,(Subscript[H, v])->1000,α->2,k->.001,Ca0->10,(Subscript[C, Area])->1,(Subscript[C, R])->2.5};

Plot[Evaluate[(Totalcostfunction/.numbs2)/.β->{1.05,1.1,1.2,1.5}],{ξ,1`,9},AxesLabel->{" k τ","total cost"},PlotLegends->{"β->1.05","β->1.1","β->1.2","β->1.5"}]

FindMinimum[{Totalcostfunction/.numbs2,ξ>0, β>0},{{β,1.3},{ξ,3}}]

More expensive energy.  You want more conversion.

numbs2={(Subscript[x, D])->.95,(Subscript[x, B])->.05, D->1000,(Subscript[C, E]) -> .04,(Subscript[H, v])->1000,α->2,k->.001,Ca0->10,(Subscript[C, Area])->1,(Subscript[C, R])->2.5};

Plot[Evaluate[(Totalcostfunction/.numbs2)/.β->{1.05,1.1,1.2,1.5}],{ξ,1`,9},AxesLabel->{" k τ","total cost"},PlotLegends->{"β->1.05","β->1.1","β->1.2","β->1.5"}]

FindMinimum[{Totalcostfunction/.numbs2,ξ>0, 5>β>0},{{β,1.3},{ξ,3}}]

Cheaper energy

numbs2={(Subscript[x, D])->.95,(Subscript[x, B])->.05, D->1000,(Subscript[C, E]) -> .025,(Subscript[H, v])->1000,α->2,k->.001,Ca0->10,(Subscript[C, Area])->1,(Subscript[C, R])->2.5};

Plot[Evaluate[(Totalcostfunction/.numbs2)/.β->{1.1,1.2,1.3,1.5}],{ξ,0.5`,9},AxesLabel->{" k τ","total cost"},PlotLegends->{"β->1.1","β->1.2","β->1.3","β->1.5"}]

FindMinimum[{Totalcostfunction/.numbs2,ξ>0, 5>β>0},{{β,1.3},{ξ,3}}]

Harder distillation, decrease α

numbs2={(Subscript[x, D])->.95,(Subscript[x, B])->.05, D->1000,(Subscript[C, E]) -> .03,(Subscript[H, v])->1000,α->1.7,k->.001,Ca0->10,(Subscript[C, Area])->1,(Subscript[C, R])->2.5};

Plot[Evaluate[(Totalcostfunction/.numbs2)/.β->{1.1,1.2,1.3,1.5}],{ξ,0.5`,10},AxesLabel->{" k τ","total cost"},PlotLegends->{"β->1.1","β->1.2","β->1.3","β->1.5"}]

Ok.

FindMinimum[{Totalcostfunction/.numbs2,ξ>0, 5>β>0},{{β,1.5},{ξ,6}}]

Harder distillation, decrease α.  We see that we have to make the residence time longer in an attempt to reduce the amount of distillation that we need to do.

numbs2={(Subscript[x, D])->.95,(Subscript[x, B])->.05, D->1000,(Subscript[C, E]) -> .03,(Subscript[H, v])->1000,α->1.6,k->.001,Ca0->10,(Subscript[C, Area])->1,(Subscript[C, R])->2.5};

Plot[Evaluate[(Totalcostfunction/.numbs2)/.β->{1.1,1.2,1.3,1.5}],{ξ,0.5`,10},AxesLabel->{" k τ","total cost"},PlotLegends->{"β->1.1","β->1.2","β->1.3","β->1.5"}]

Ok.

FindMinimum[{Totalcostfunction/.numbs2,ξ>0, 5>β>0},{{β,1.5},{ξ,12}}]

What if the reaction rate is faster?

numbs2={(Subscript[x, D])->.95,(Subscript[x, B])->.05, D->1000,(Subscript[C, E]) -> .03,(Subscript[H, v])->1000,α->2,k->.0012,Ca0->10,(Subscript[C, Area])->1,(Subscript[C, R])->2.5};

Plot[Evaluate[(Totalcostfunction/.numbs2)/.β->{1.1,1.3,1.5,1.7}],{ξ,0.5`,10},AxesLabel->{" k τ","total cost"},PlotLegends->{"β->1.1","β->1.3","β->1.5","β->1.7"}]

Ok.

FindMinimum[{Totalcostfunction/.numbs2,ξ>0, 5>β>0},{{β,1.5},{ξ,6}}]

We see that we will let k τ get bigger because this cuts distillation costs while not needing a bigger reactor,

Here is distillation costs to see the effect of α

numbs3={(Subscript[x, D])->.95,(Subscript[x, B])->.05, D->1000,(Subscript[C, E]) -> .03,(Subscript[H, v])->1000,ξ->4,k->.001,Ca0->10,(Subscript[C, Area])->1,(Subscript[C, R])->2.5};

Here we see that if the distillation is made easier, the total separation costs go down pretty quickly and the reflux ratio is further above the minimum.

Plot[Evaluate[distcosts/.numbs3/.α->{1.6,1.8,2,2.2}],{β,1.02`,2},AxesLabel->{" reflux factor","distillation costs"}]

Table[{α,FindMinimum[{((distcosts/.numbs3)), 5>β>0},{{β,1.5}}]},{α,1.6,2.2,.1}]

Here is distillation costs to see the effect of kτ

numbs3={(Subscript[x, D])->.95, (Subscript[x, B])->.05,D->1000,(Subscript[C, E]) -> .03,(Subscript[H, v])->1000,k->.001,α->2,Ca0->10,(Subscript[C, Area])->1,(Subscript[C, R])->2.5};

Here we see that if the reaction goes further to completion, the total separation costs go down pretty quickly and the reflux ratio is further above the minimum.

Plot[Evaluate[distcosts/.numbs3/.ξ->{2,3,4,5}],{β,1.02`,2},AxesLabel->{" reflux factor","distillation costs"},PlotLegends->{"ξ->2","ξ->3","ξ->4","ξ->5"}]

Table[{ξ,FindMinimum[{((distcosts/.numbs3)), 5>β>0},{{β,1.5}}]},{ξ,1.5,8,.5}]

Comments

We find that the optimum dimensionless residence time is in the range of 2-4.  This is a little less than what the first order reaction conversion function shows, the conversion is increasing rapidly and then begins to slow down at 4.

The effect of reaction rate constant is profound.  If the reaction is quicker, we can get closer to completion in a smaller, cheaper reactor which reduces the distillation costs.  Thus we confirm that the cost increase at higher kτ, is that V is getting too big -- not that k is too big.

If the relative volatility of the components strongly influences the total cost of the process.  As the volatility gets closer to unity, the distillation gets harder and we are willing to use a bigger reactor to get a longer residence time.

If a separate optimization calculation is done for the reflux ratio factor, which is the main parameter that is adjusted for finding the optimal balance between number of trays and energy costs, the optimum is in the range of 1.1.  If energy costs go up, this goes down.