Maximum production from a batch reactor

Mark J McCready
Professor of Chemical and Biomolecular Engineering
University of Notre Dame
mjm@nd.edu
chemeprof.com
8/27/21

Reaction time plus downtime are considered

When you run a batch reactor, there is reaction time and down time necessary for filling, emptying and cleaning.

The optimal reaction time for maximum production will  be a function of the reaction kinetics and the downtime.

In this notebook, we investigate this issue for different kinetic expressions and find that indeed, there is an optimal time for running the reaction for some reaction orders.

Conclusions

For a zero order reaction, the maximum production (total moles of product/ total time) will occur if the reaction is allowed to reach completion no matter what the downtime is.  This is because the reactor is always running at maximum rate.

For a first order reaction, if the reaction is slow enough that k*downtime is less than order ~1,  the optimal reaction time is approximately Sqrt[ 2/ (downtime*k)], where k is the first order reaction rate constant.   For faster reactions the optimal time decreases faster with k than this relation, but there does not appear to be a simple scaling relation.

For a 2nd order reaction, the optimal reaction time is (Sqrt[(1/(c0*tfixed*k))])essentially analogous to the first order, slow reaction result.

zero order kinetics

A-- > M
r = k

The production (moles/time) for a batch reactor is
Cmfinal V/(treaction + tdowntime)

We will call the downtime:  tfixed.
We will pick the reaction time as: α* tfixed.
Thus the final product concentration is Cmf = k α tfixed

Production/volume/feed concentration

Production0  =k α tfixed / (1 + α)/tfixed/Cf

Let's see what this function looks like

Plot[Production0/.{k->1,Cf->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","production/volume"}]

We can see that there is no max.  The function will asymptotically approach k/Cf

Limit[Production0,α->Infinity]

So for 0 order kinetics, you would let the reaction go to completion.  Of course no reaction will remain at 0 order as the reactants are exhausted.

first order kinetics

A-- > M
r = k Ca

The normalized (for unit initial concentration)production (moles/time/Volume/Caf) for a batch reactor is
Cmfinal  /(treaction + tdowntime)/Caf

We will call the downtime:  tfixed.
The reaction time is α*tfixed

This is dimensionless production, in a (total) time (α*tfixed+tfixed).  The maximum production should occur for a specific value of α that is a function of the rate constant.   To keep this as simple as possible, tfixed is always kept the same.

Production1  = Caf(1- Exp[-k α tfixed])  V/(1 + α)/tfixed/V/Caf/k

The concentration at the optimal time can be obtained from this expression.

Concentration  = (1- Exp[-k α tfixed])

Let's see what this function looks like.  You would never need to wait past about α tfixed =~ 3/k, to get most of the way to completion.  
We shall see that we might want to cut off the reaction much sooner though!

Plot[Concentration/.{k->1, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless concentration"},PlotStyle->Blue, PlotLegends->{"k*tfixed =1"}]

We see that we would get more total production of “M” if we stop the  reaction at α*tfixed = tfixed, that is, α =1.  In doing this we are taking advantage of the fast reaction that is occurring for shorter times.  (This does increase the separation costs however.)

Plot[Production1/.{k->1, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless production rate"},PlotStyle->Blue, PlotLegends->{"k*tfixed =1"}]

Increasing values of the rate constant.  As k is increased, the optimal time decreases and even with a faster reaction, the optimal final concentration is decreasing.   

While the dimensionless production is decreasing, the dimensional production will increase for a faster reaction because the rate constant is used to make the production dimensionless!

Plot[Concentration/.{k->2, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless concentration"},PlotStyle->Red, PlotLegends->{"k*tfixed =2"}]

Plot[Production1/.{k->2, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless production rate"},PlotStyle->Red, PlotLegends->{"k*tfixed =2"}]

Plot[Concentration/.{k->5, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless concentration"},PlotStyle->Orange, PlotLegends->{"k*tfixed =5"}]

Plot[Production1/.{k->5, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless production rate"},PlotStyle->Orange, PlotLegends->{"k*tfixed =5"}]

Plot[Concentration/.{k->10, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless concentration"},PlotStyle->Purple, PlotLegends->{"k*tfixed =10"}]

Plot[Production1/.{k->10, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless production rate"},PlotStyle->Purple, PlotLegends->{"k*tfixed =10"}]

Show[%175,%173,%171,%169]

Now for decreasing values of k.

Plot[Concentration/.{k->1, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless concentration"},PlotStyle->Blue, PlotLegends->{"k*tfixed =1"}]

Plot[Production1/.{k->1, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless production rate"},PlotStyle->Blue, PlotLegends->{"k*tfixed =1"}]

Plot[Concentration/.{k->1, Caf->1,tfixed->2},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless concentration"},PlotStyle->Red, PlotLegends->{"k*tfixed =1/2"}]

Plot[Production1/.{k->1, Caf->1,tfixed->2},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless production rate"},PlotStyle->Red, PlotLegends->{"k*tfixed =1/2"}]

Plot[Concentration/.{k->1/5, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless concentration"},PlotStyle->Orange, PlotLegends->{"k*tfixed =1/5"}]

Plot[Production1/.{k->1/5, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless production rate"},PlotStyle->Orange, PlotLegends->{"k*tfixed =1/5"}]

Plot[Concentration/.{k->1/10, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless concentration"},PlotStyle->Purple, PlotLegends->{"k*tfixed =1/10"}]

Plot[Production1/.{k->1/10, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","dimensionless production rate"},PlotStyle->Purple, PlotLegends->{"k*tfixed =1/10"}]

Show[%183,%185,%188,%190]

So if the  reaction rate is slower, the optimal time needs to increase.

It requires some thought to explain why for a very slow reaction, you want to stop at α~4.  For a slow reaction, the downtime is irrelevant.  Thus once the reaction starts to slow down very much at all, it is best to dump the contents and start with a higher concentration reactants.

rate1 = D[Concentration,α]

Plot[Concentration/.{k->1/10, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","final concentration"},PlotStyle->Purple, PlotLegends->{"k*tfixed =1"},PlotLabel->{"first order reaction"}]

Here is the reaction rate:

Plot[rate1/.{k->1/10, Caf->1,tfixed->1},{α,0,10},PlotRange->All, AxesLabel->{"t/tfixed","reaction rate"},PlotStyle->Green, PlotLegends->{"k*tfixed =1"},PlotLabel->{"first order reaction"}]

Let’s directly get the optimal values of production and concentration

We will use the function “Maximize”.  You can find the description of this function using “??Maximize”

We need to restrict the range of α over which the function will search for a maximum.

𝒟=Interval[{0,20}]

For the Red curve on the page above.  The first entry in the function is what is being maximized, the second is the restriction for the range of the variable, α.  Then outside of the bracket is (again), α which is the variable that is being adjusted to find the maximum.

Maximize[{Production1/.{k->1/2, Caf->1,tfixed->1},{α}∈𝒟},α]

Oh, you want numbers....  You can see that it works.

N[%]

Here is the entire range:

LogLogPlot[N[Maximize[{Production1/.{ Caf->1,tfixed->1},{α}∈𝒟},α]][[1]],{k,.01,20},PlotRange->All,AxesLabel->{" k* tfixed","maximum dimensionless production"}, PlotLabel->{"First order reaction with downtime"}]

Here is the dimensional production.  It does increase as the reaction rate is increased.

LogLogPlot[k N[Maximize[{Production1/.{ Caf->1,tfixed->1},{α}∈𝒟},α]][[1]],{k,.01,20},PlotRange->All,AxesLabel->{" k* tfixed","maximum dimensional production"}, PlotLabel->{"First order reaction with downtime"}]

Here is the optimal reaction time  as a function of the rate constant,

LogLogPlot[α/.N[Maximize[{Production1/.{ Caf->1,tfixed->1},{α}∈𝒟},α]][[2]],{k,.01,20},PlotRange->All,AxesLabel->{" k* tfixed","optimal reaction time"}, PlotLabel->{"First order reaction with downtime"},PlotLegends->{"full solution"}]

Perhaps this is interesting.  This is the concentration at the optimal batch time.

LogLinearPlot[(Concentration/.tfixed->1)/.N[Maximize[{Production1/.{ Caf->1,tfixed->1},{α}∈𝒟},α]][[2]],{k,.01,20},PlotRange->All,AxesLabel->{" k* tfixed","optimal concentration"}]

Simple scaling for 1st order reaction maximum.

Let's see if we can find the maximum analytically

diff1 =D[Production1,α]

time1max =Solve[diff1==0,α]

Ok use Reduce

Reduce[diff1==0,α]

??Reduce

This answer is not useful.   We don’t have a way to determine (C[1]).  If it is 0, then expression is always negative.  

We need to actually do some work to get a reasonable answer that will work in particular ranges

It may not be obvious, but this answer does not work.  It is always negative.

We need to actually do some work to get a reasonable answer that will work in particular ranges

Together[diff1]

So the middle expression needs to be 0 for the function to be a maximum.  It is apparent why there is a problem,

exp1 = (-1+(E^(k t))-k t-k tfixed)

The problem is that t is in both an exponential and a linear term.  Let's consider the case of k t being large enough so that Exp[kt] >> kt.

The equation to solve is thus:

-1 + Exp[ k t] - k tfixed ==0

ans11=Solve[%,t]

This answer doesn’t match the real result above!

t1max = (t/.ans11[[1]])/.(C[1])->0

t1max/.{k->10, tfixed->1.}

Maximize[{Production1/.{k->10, Caf->1,tfixed->1},{α}∈𝒟},α]//N

slow reaction

do an expansion around k=0

ex1 =Together[diff1]

ex2 = ex1/.k-> K /tfixed

ex3 = Series[ex2, {K,0,2}]

ex4 = Normal[Series[ex3, {K,0,1}]]

alfans  = Solve[ex4 ==0,α]

alfmax = α/.alfans[[2]]

Plot[alfmax,{K,0.001,1}]

If we realize that the numerator leaves only the Sqrt[2 K], We see that the maximum production batch time scales as Sqrt[2/K],  K = tfixed*k.

Plot[-1 + Sqrt[2/k],{k,.1,1}, PlotStyle->Red]

Show[%,%%]

Let’s see how well the scaling works.  I adjusted the constant to make this fit

LogLogPlot[Sqrt[ 2*.92/k],{k,.01,20},PlotStyle->Red,PlotLegends->{"1/\!\(\*SqrtBox[\(tfixed\\ k\)]\) scaling"}]

Show[%7,%]

This simplification works well for slow reactions

Even if we keep more of the approximate expression, we don’t do any better!

LogLogPlot[(1.05 alfmax)/.K-> k ,{k,.01,20},PlotStyle->Purple,PlotLegends->{""((-K+(Sqrt[2 K+(K^2)]))/K) "scaling"}]

Show[%66,%]

We can show that the original simplification did not give a useful result:

LogLogPlot[((Log[1+k *10])/k),{k,.01,20},PlotStyle->Green,PlotLegends->{"\!\(\*FractionBox[\(Log[1 + k\\ *tfixed]\), \(k\)]\)"}]

We see that for slow reactions, the limiting scaling works pretty well.  However above order 1, the attempt to find a relation did not work.  Note that this total range exceeds the limits of practicality.  If the reaction time is very short compared to the downtime, then you would probably try to run the reactor continuously.

Show[%78,%]

2nd order kinetics

So we consider A + B -> M. 
r = k Ca^2

First solve the reactor equation for this case:

ans12nd=DSolve[{D[ca[t],t] == - k ca[t] ca[t],ca[0]==caf},ca[t],t]

To get the product we need  Caf - Ca

Production2 =( (caf-ca[t]/.ans12nd[[1]])V/caf/V/(1 + α)/tfixed/(k caf))/.t-> α tfixed

Plot[Production2 /.{caf->1,tfixed->1,k->.1},{α,0,10}]

FullSimplify[Production2]

LogLogPlot[ N[Maximize[{Production2/.{ caf->1,tfixed->1},{α}∈𝒟},α]][[1]],{k,.01,20},PlotRange->All,AxesLabel->{" k* tfixed","maximum dimensionless production"}, PlotLabel->{"2nd order reaction with downtime"}]

Dimensional production

LogLogPlot[k  N[Maximize[{Production2/.{ caf->1,tfixed->1},{α}∈𝒟},α]][[1]],{k,.01,20},PlotRange->All,AxesLabel->{" k* tfixed","maximum dimensional production"}, PlotLabel->{"2nd order reaction with downtime"}]

Here is the optimal reaction time  as a function of the rate constant,

secondorderplot =LogLogPlot[α/.N[Maximize[{Production2/.{ caf->1,tfixed->1},{α}∈𝒟},α]][[2]],{k,.01,20},PlotRange->All,AxesLabel->{" k* tfixed","optimal reaction time"}, PlotLabel->{"Second order reaction with downtime"},PlotLegends->{"full solution"}]

This looks simple...

D[Production2,α]

Yes, it is.  This is the maximum.

Solve[%==0,α]

Let’s check

LogLogPlot[Sqrt[ 1/k],{k,.01,20},PlotStyle->Red,PlotLegends->{"1/(Sqrt[tfixed*k]) scaling"}]

Show[secondorderplot,%]