Successive First-Order Reversible Reactions

Authors: Benjamin J. McCoy and Housam BINOUS

Benjamin J. McCoy
Department of Chemical Engineering 
Louisiana State University Baton Rouge, La. 70803 
Email: bjmccoy@lsu.edu 

Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reaction Networks

We consider successive first-order reversible reactions:

	A1 \(TraditionalForm\`\( ⇌ \+\(k\_21\)\%\(k\_12\)\)\) A2 \(TraditionalForm\`\( ⇌ \+\(k\_32\)\%\(k\_23\)\)\) A3 \(TraditionalForm\`\( ⇌ \+\(k\_43\)\%\(k\_34\)\)\)... \(TraditionalForm\`\( ⇌ \+\(k\_\(n\ \ n - 
                1\)\)\%\(k\_\(n - 1\ \ n\)\)\)\) An

One can write the kinetic expressions where Ai represents the molar concentration:

	dA1/dt =  -k12 A1 + k21 A2
	dA2/dt =  k12 A1 - (k21 + k23)A2  + k32 A3
	dA3/dt =  k23 A2 - (k32 + k34)A3  + k43 A4
	:
	dAn/dt =  kn-1 n A n-1 - kn n-1An  
	
Adding the equations provides the differential mass balance:
	d(A1 + A2 + A3 +...+ An)/dt = 0.

By setting the left hand side of the differential equations equal to zero, one can solve the simultaneous linear algebraic equations and get the steady state solutions.

 With initial conditions on A1, A2, A3,... An, one can solve the simultaneous ordinary differential equations and get the transient solution. Here, A1(t=0) = 1.0 and other initial concentrations are zero. Values of the rate constants are chosen as 1.0 for forward and 0.9 for reverse reactions.

Final equilibrium solution for n reactions 

We plot the steady state solutions.

<< Graphics`Graphics`

Clear[soleq, eqns]; n = 1000; A1o = 1. ; A[0] = 0; A[n + 1] = 0; 
    A[n] = (A1o - Sum[A[j], {j, 1, n - 1}]);
k[0, 1] = 0; k[1, 0] = 0; k[n, n + 1] = 0; k[n + 1, n] = 0;
Table[k[i, (i + 1)] = 1. , {i, 1, n}]; 
    Table[k[(i + 1), i] =  .9, {i, 1, n}];

eqns = 
        Table[0 == 
            k[j - 1, j] A[j - 1] - (k[j, j - 1] + k[j, j + 1] ) A[j] + 
              k[j + 1, j] A[j + 1], {j, 1, n - 1}];

soleq = 
        NSolve[
eqns, 
Table[
            A[j], {j, 1, n - 1}]];

A[n] = (A1o - Sum[A[j], {j, 1, n - 1}]) /. soleq

Flatten[Table[A[j], {j, 1, n}] /. soleq]

LogLogListPlot[
        Flatten[Table[A[j], {j, 1, n}] /. soleq]];

Transient solution for n reactions

We plot the concentration profile for species:  (Subscript[A, 100]), (Subscript[A, 200]), (Subscript[A, 300]), (Subscript[A, 400]), (Subscript[A, 500]), (Subscript[A, 600]), (Subscript[A, 700]), (Subscript[A, 800]) and (Subscript[A, 900]).

We also plot the concentration profile for (Subscript[A, 1]), and finally we plot the concentration profile for (Subscript[A, 1000]).

Clear[soleq, eqns]; tf = 50000;
n = 1000; A1o = 1. ; (Subscript[A, 0])[t] = 0; (Subscript[A, n + 1])[t] = 0; 
    (Subscript[A, n])[t] = (A1o - 
          Sum[(Subscript[A, j])[t], {j, 1, n - 1}]);
k[0, 1] = 0; k[1, 0] = 0; k[n, n + 1] = 0; 
    k[n + 1, n] = 0;
Table[k[i, (i + 1)] = 1. , {i, 1, n}]; 
    Table[k[(i + 1), i] =  .9, {i, 1, n}];

eqns[j_] = 
        D[(Subscript[A, j])[t], {t, 1}] == 
          k[j - 1, j] (Subscript[A, j - 1])[t] - (k[j, j - 1] + k[j, j + 1] ) 
              (Subscript[A, j])[t] + k[j + 1, j] (Subscript[A, j + 1])[t];

soleq = 
        NDSolve[
Flatten[{Table[eqns[j], {j, 1, n - 1}], 
              (Subscript[A, 1])[0] == 1, 
Table[
                (Subscript[A, j])[0] == 0, {j, 2, 
                  n - 1}]}], 
Table[
            (Subscript[A, j])[t], {j, 1, n - 1}], {t, 0, tf}];

Plot[
      Evaluate[Table[(Subscript[A, j])[t], {j, 100, n - 1, 100}] /. soleq], {t, 0, tf}, 
      Frame -> True, 
      PlotStyle -> Table[Hue[100 j/( n)], {j, 1, n}], 
      DefaultFont -> {"Symbol-Bold", 10}, 
      FrameLabel -> {"t", "Concentration"}, 
      PlotRange -> {{0, 15000}, {0, 0.01}}]

Plot[Evaluate[(Subscript[A, 1])[t] /. soleq], {t, 0, tf}, Frame -> True, 
      PlotStyle -> Hue[ 1/( n)], 
      DefaultFont -> {"Symbol-Bold", 10}, 
      FrameLabel -> {"t", "Concentration"}, 
      PlotRange -> {{0, 400}, {0, 1}}, 
      PlotPoints -> 900]

Plot[Evaluate[(Subscript[A, 1000])[t] /. soleq], {t, 0, tf}, Frame -> True, 
      PlotStyle -> Hue[ 1000/( n)], 
      DefaultFont -> {"Symbol-Bold", 10}, 
      FrameLabel -> {"t", "Concentration"}, 
      PlotRange -> {{0, 15000}, {0, 0.15}}]