Analysis of a batch reactor

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

Mathematica example

We wish to solve illustration 8.2 in Hill and Root.  This involves some complex kinetics to produce “zeolite A” in a batch reactor

The mass balance for a batch reactor is:

Rate of change of species A in the reactor = rate at which A is consumed by reaction

V ((d (Subscript[C, A]))/dt) = - (Subscript[r, A]) V,

where V is the reactor volume and (Subscript[C, A])is the concentration of the reactant, A.

To integrate this expression, we will need “(Subscript[r, A])” expressed in terms of the concentration of A.

The texts gives this expression

From this we need to express (Subscript[C, Z]) and [(SuperMinus[OH])] in terms of (Subscript[C, A]) using the reaction stoichiometry.  

We see that the equation above can be rearranged to be 

 dt = -  ((d (Subscript[C, A]))/(Subscript[r, A])) = ((d (Subscript[C, A]))/(F[(Subscript[C, A])])), where F[CA] is the right side of the complete rate expression.
 
 H&R suggest that this problem is solved by considering (Subscript[f, A]), the fractional conversion of A
 
 (Subscript[f, A]) ≡(((Subscript[C, A0])-(Subscript[C, A]))/(Subscript[C, A0]))

We can solve for CA

ans1 = Solve[ (Subscript[f, A]) == ((Subscript[C, A0])-(Subscript[C, A]))/(Subscript[C, A0]), (Subscript[C, A])]

which is equation (C) on page 228.  The “ans1” is the “name” of the result we have obtained from solving the simple equation.  We will use it to substitute for (Subscript[C, A]) below.  Also note that a single “=” assigns the variable on the left to be replaced by whatever is on the right.  In the command a “==” is used to denote an “equation” that can be solved as we wish.  These are crucial distinctions in Mathematica (and more generally if you think  about it!)

The zeolite concentration, (Subscript[C, Z]) , can be written in terms of how much A has reacted away as the fractional conversion of A gives the increase in the relative concentration of Z.

ans2 = Solve[ (Subscript[f, A]) == ((Subscript[C, Z])-(Subscript[C, Z0]))/(Subscript[C, A0]), (Subscript[C, Z])]

The [OH-]  remains constant.

Instead of substituting numbers, let’s get a full analytic answer that allows any values of the concentrations and rate constants.

reactor equation in terms of (Subscript[f, A])

dt = -  ((d (Subscript[C, A]))/(Subscript[r, A])) =-(Subscript[C, A0]) ((d (Subscript[f, A]))/(Subscript[r, A])),

Now let’s solve this.

Here is the rate expression .

arg1= ((Subscript[k, 2])( OH^(Subscript[k, 1])) (Subscript[C, A]) (Subscript[C, Z]))/((Subscript[k, 3]) (Subscript[C, Z])+ ( 1 + (Subscript[k, 3])) (Subscript[C, A]))

Here is how we substitute for (Subscript[C, A]) and (Subscript[C, Z]).  
The “/.” is the substitution command. 
The -> in ans2 and ans1 gives  the replacement that will be done for variables already in the equation.   
The [[1]] is a Mathematica “designator” indicating that part 1 of the expression is to be use for the substitution. 
Since there a two substitutions, we do the operation twice.

arg2 =(arg1/.ans1[[1]])/.ans2[[1]]

This is easily integrated as an indefinite integral

ans3 =(Subscript[C, A0]) Integrate[ 1/arg2,(Subscript[f, A])]

This gives the time as for a final fractional conversion of (Subscript[f, Af]) staring at (Subscript[f, A])=0

reactiontime =(ans3/.(Subscript[f, A])-> (Subscript[f, Af]) )- (ans3/.(Subscript[f, A])-> 0)

We can put in numbers to check H&R

numbs ={(Subscript[k, 1])->2.36, (Subscript[k, 2])->.625, OH->1.5, (Subscript[C, A0])-> 24, (Subscript[k, 3])->.36, (Subscript[C, Z0])->1}

time1 =reactiontime/.numbs

For a conversion of .98

finalanswer = time1/.{(Subscript[f, Af])->.98}

This is in ksec,  Since the remnant from imaginary  parts of the integral is not exactly 0, we can “Chop”

Chop[finalanswer]

Which matches H&R.    Our answer works for any values of the constants.
In fact, why not put the correct units in with the numbers!

numbs2 ={(Subscript[k, 1])->2.36, (Subscript[k, 2])->.625 /ks, OH->1.5 kmol/m^3, (Subscript[C, A0])-> 24 kg/m^3, (Subscript[k, 3])->.36, (Subscript[C, Z0])->1 kg/m^3}

time2 =reactiontime/.numbs2

finalanswer2 = time2/.{(Subscript[f, Af])->.98}

FullSimplify[%]

We see the answer from above. 
However the result is rather interesting and reveals that the rate expression (A) as given in the text, is dimensionally inconsistent.  We see this now as the quantity, (Subscript[k, 2])( OH^(Subscript[k, 1])), should together have the dimension of 1/time.  But for this to occur, the k2 would need a different dimension.

How about a definite integral?

Mathematica can calculate definite integrals, but it will try to give the most general result for any possible values of the constants.  In this attempt, the time to get a result could get very excessive!  Here I show how I make the time less by designating the constants to be positive and hence Real by Mathematica convention.)

arg2

With the “assumptions” as designated the time to get an answer is a little more than 1 sec on my machine

Timing[(Subscript[C, A0])Integrate[ 1/arg2,{(Subscript[f, A]),0, (Subscript[f, Af])},Assumptions->{(Subscript[k, 1])>0,(Subscript[k, 2])>0,(Subscript[k, 3])>0,OH>0,(Subscript[C, A0])>0, (Subscript[C, Z0])>0,(Subscript[f, A])>0, (Subscript[f, Af])>0}]]

You can see that Mathematica is giving an answer if the fractional conversion is less than 1.  Presumably the expression is singular if (Subscript[f, Af]), plus, we know that in a physical sense, the reaction time would be infinite anyway.

Let’s check the answer.  Here I use a different designation for the expression I want to work on,  “%” tells me to pick the last answer.  the “[[2]]” is the 2nd part, which is the integral, not the timing!

%[[2]]/.numbs2

We pick (Subscript[f, Af])=.98

%/.(Subscript[f, Af])->.98

Which is the expected result.  

Suppose that I just drop one key assumption!, k2 can be negative or imaginary.

Timing[(Subscript[C, A0])Integrate[ 1/arg2,{(Subscript[f, A]),0, (Subscript[f, Af])},Assumptions->{(Subscript[k, 1])>0,(Subscript[k, 3])>0,OH>0,(Subscript[C, A0])>0, (Subscript[C, Z0])>0,(Subscript[f, A])>0, (Subscript[f, Af])>0}]]

Which is the expected result.  

If I don’t restrict the rate constants to be real.  The time increases a bit

Timing[(Subscript[C, A0])Integrate[ 1/arg2,{(Subscript[f, A]),0, (Subscript[f, Af])},Assumptions->{(Subscript[C, A0])>0, (Subscript[C, Z0])>0,(Subscript[f, A])>0, (Subscript[f, Af])>0}]]

Actually if you re-run these, Mathematica learns and saves some answers.  So the time goes down.  This started at >15 sec.  You can also see that MMA has given a condition that must be true for this to be the correct answer.

Timing[(Subscript[C, A0])Integrate[ 1/arg2,{(Subscript[f, A]),0, (Subscript[f, Af])},Assumptions->{(Subscript[C, A0])>0, (Subscript[f, A])>0, (Subscript[f, Af])>0}]]

Finally, making only 1 assumption the time is much longer, but  the solution has many more restrictions

Timing[(Subscript[C, A0])Integrate[ 1/arg2,{(Subscript[f, A]),0, (Subscript[f, Af])},Assumptions->{(Subscript[f, A])>0}]]

Here are the restrictions.  But since we know what we needed for the variables,  we could avoid this situation.

%149[[2,2]]

But  this answer is fine for our problem!

FullSimplify[(%149[[2,1]]/.numbs2 )/.(Subscript[f, Af])->.98]

The 2nd part of this example is very straight forward!  However to get more insight into an optimal batch time for a reactor that has downtime to empty, clean and refill see another Mathematica notebook:

production_maximum_batch_reactor.nb