<< Statistics`LinearRegression`

In CRE there are some characteristic graphs, which summarize a large body of information and which are essential (and visually very powerful!) for the analysis and design of chemical reactors. We have seen seven such graphs in Chs. 1-8: 

ln k vs. 1/T, 
(Subscript[C, A]) vs. t, 
kτ vs. (Subscript[X, A])
-1/(Subscript[r, A]) vs. (Subscript[X, A])
(τNMFR/τPFR) vs. 1-(Subscript[X, A])
φ (R/A) vs. (Subscript[C, A]) (for parallel reactions)
(Subscript[C, R]) vs. (Subscript[X, A]) (for consecutive reactions)

We are about to meet another such graph: 

(Subscript[X, A]) vs. T (for non-isothermal reactors)

Let's work through the important Examples 9.2 and 9.3 in Levenspiel.

gasConst = 8.314;
temp1 = 298;
deltaG298 = -14130;
deltaHrxn = -75300;
equilConst298 = Exp[-deltaG298/gasConst/temp1]

equilConst = 
        equilConst298*
          Exp[-deltaHrxn/
                gasConst*(1/temp - 
                  1/temp1)];
xAe = equilConst/(equilConst + 1);
Plot[xAe, {temp, 278, 420}, 
        GridLines -> Automatic];

xAeAt78 = xAe /. temp -> 78 + 273

OK. The above graph is Figure E9.2. Let's proceed now to Example 9.3. We need the complete rate equation for a 1st order reversible reaction, on the basis of which we can construct the conversion vs. temperature graph.

xAf1 = 0.581;
time1 = 1;
temp1 = 65 + 273;
xAe1 = xAe /. temp -> temp1;
k11 = %/time1*(-Log[1 - xAf1/xAe1])

xAf2 = 0.60;
time2 = 10;
temp2 = 25 + 273;
xAe2 = xAe /. temp -> temp2;
k12 = %/time2*(-Log[1 - xAf2/xAe2])

activEnergy1 = gasConst*Log[k12/k11]/(1/temp1 - 1/temp2)

preExp1 = k12/Exp[-activEnergy1/gasConst/temp2]

k1 = 
      preExp1*Exp[-activEnergy1/gasConst/
            temp]
k2 = k1/equilConst

Now we can construct the X vs. T graph. Let's first solve the relevant equations for a space time of 5, 10 and 20 minutes in a PFR (BR), and then let's do the same for an MFR.

tau1 = 10;
tau2 = 5;
tau3 = 20;
tempData = 
      Table[temp[i] = 298 + 10  i, {i, 1, 7}]
xAe1 = xAe /. temp -> temp[1];
xAe2 = xAe /. temp -> temp[2];
xAe3 = xAe /. temp -> temp[3];
xAe4 = xAe /. temp -> temp[4];
xAe5 = xAe /. temp -> temp[5];
xAe6 = xAe /. temp -> temp[6];
xAe7 = xAe /. temp -> temp[7];
k11 = 
        preExp1*Exp[-activEnergy1/gasConst/
              temp[1]];
k12 = 
        preExp1*Exp[-activEnergy1/gasConst/
              temp[2]];
k13 = 
        preExp1*Exp[-activEnergy1/gasConst/
              temp[3]];
k14 = 
        preExp1*Exp[-activEnergy1/gasConst/
              temp[4]];
k15 = 
        preExp1*Exp[-activEnergy1/gasConst/
              temp[5]];
k16 = 
        preExp1*Exp[-activEnergy1/gasConst/
              temp[6]];
k17 = 
        preExp1*Exp[-activEnergy1/gasConst/
              temp[7]];
xAPFR11 = 
        xAe1*(1 - Exp[-k11*tau1/xAe1]);
xAPFR21 = 
        xAe2*(1 - Exp[-k12*tau1/xAe2]);
xAPFR31 = 
        xAe3*(1 - Exp[-k13*tau1/xAe3]);
xAPFR41 = 
        xAe4*(1 - Exp[-k14*tau1/xAe4]);
xAPFR51 = 
        xAe5*(1 - Exp[-k15*tau1/xAe5]);
xAPFR61 = 
        xAe6*(1 - Exp[-k16*tau1/xAe6]);
xAPFR71 = 
        xAe7*(1 - Exp[-k17*tau1/xAe7]);
xAPFR12 = 
        xAe1*(1 - Exp[-k11*tau2/xAe1]);
xAPFR22 = 
        xAe2*(1 - Exp[-k12*tau2/xAe2]);
xAPFR32 = 
        xAe3*(1 - Exp[-k13*tau2/xAe3]);
xAPFR42 = 
        xAe4*(1 - Exp[-k14*tau2/xAe4]);
xAPFR52 = 
        xAe5*(1 - Exp[-k15*tau2/xAe5]);
xAPFR62 = 
        xAe6*(1 - Exp[-k16*tau2/xAe6]);
xAPFR72 = 
        xAe7*(1 - Exp[-k17*tau2/xAe7]);
xAPFR13 = 
        xAe1*(1 - Exp[-k11*tau3/xAe1]);
xAPFR23 = 
        xAe2*(1 - Exp[-k12*tau3/xAe2]);
xAPFR33 = 
        xAe3*(1 - Exp[-k13*tau3/xAe3]);
xAPFR43 = 
        xAe4*(1 - Exp[-k14*tau3/xAe4]);
xAPFR53 = 
        xAe5*(1 - Exp[-k15*tau3/xAe5]);
xAPFR63 = 
        xAe6*(1 - Exp[-k16*tau3/xAe6]);
xAPFR73 = xAe7*(1 - Exp[-k17*tau3/xAe7]);

xAMFR11 = 
        xAe1*k11*tau1/(xAe1 + k11*tau1);
xAMFR21 = 
        xAe2*k12*tau1/(xAe2 + k12*tau1);
xAMFR31 = 
        xAe3*k13*tau1/(xAe3 + k13*tau1);
xAMFR41 = 
        xAe4*k14*tau1/(xAe4 + k14*tau1);
xAMFR51 = 
        xAe5*k15*tau1/(xAe5 + k15*tau1);
xAMFR61 = 
        xAe6*k16*tau1/(xAe6 + k16*tau1);
xAMFR71 = 
        xAe7*k17*tau1/(xAe7 + k17*tau1);
xAMFR12 = 
        xAe1*k11*tau2/(xAe1 + k11*tau2);
xAMFR22 = 
        xAe2*k12*tau2/(xAe1 + k12*tau2);
xAMFR32 = 
        xAe3*k13*tau2/(xAe1 + k13*tau2);
xAMFR42 = 
        xAe4*k14*tau2/(xAe1 + k14*tau2);
xAMFR52 = 
        xAe5*k15*tau2/(xAe1 + k15*tau2);
xAMFR62 = 
        xAe6*k16*tau2/(xAe1 + k16*tau2);
xAMFR72 = 
        xAe7*k17*tau2/(xAe1 + k17*tau2);
xAMFR13 = 
        xAe1*k11*tau3/(xAe1 + k11*tau3);
xAMFR23 = 
        xAe2*k12*tau3/(xAe1 + k12*tau3);
xAMFR33 = 
        xAe3*k13*tau3/(xAe1 + k13*tau3);
xAMFR43 = 
        xAe4*k14*tau3/(xAe1 + k14*tau3);
xAMFR53 = 
        xAe5*k15*tau3/(xAe1 + k15*tau3);
xAMFR63 = 
        xAe6*k16*tau3/(xAe1 + k16*tau3);
xAMFR73 = xAe7*k17*tau3/(xAe1 + k17*tau3);

data11 = 
        List[{temp[1], xAPFR11}, {temp[2], xAPFR21}, {temp[3], 
            xAPFR31}, {temp[4], xAPFR41}, {temp[5], xAPFR51}, {temp[6], 
            xAPFR61}, {temp[7], xAPFR71}];
data12 = 
        List[{temp[1], xAPFR12}, {temp[2], xAPFR22}, {temp[3], 
            xAPFR32}, {temp[4], xAPFR42}, {temp[5], xAPFR52}, {temp[6], 
            xAPFR62}, {temp[7], xAPFR72}];
data13 = 
        List[{temp[1], xAPFR13}, {temp[2], xAPFR23}, {temp[3], 
            xAPFR33}, {temp[4], xAPFR43}, {temp[5], xAPFR53}, {temp[6], 
            xAPFR63}, {temp[7], xAPFR73}];
dataPlot11 = 
        ListPlot[data11, 
          PlotStyle -> PointSize[0.02]];
dataPlot12 = 
        ListPlot[data12, 
          PlotStyle -> PointSize[0.02]];
dataPlot13 = 
        ListPlot[data13, 
          PlotStyle -> PointSize[0.02]];
dataPlot11p = 
      Regress[data11, Table[temp^i, {i, 0, 3}], temp, 
        RegressionReport -> RSquared]
dataPlot12p = 
      Regress[data12, Table[temp^i, {i, 0, 3}], temp, 
        RegressionReport -> RSquared]
dataPlot13p = 
      Regress[data13, Table[temp^i, {i, 0, 3}], temp, 
        RegressionReport -> RSquared]
xAPFRvsT1 = 
        Fit[data11, Table[temp^i, {i, 0, 3}], 
          temp];
Plot[xAPFRvsT1, {temp, 298, 358}];
Show[%, dataPlot11];
xAPFRvsT2 = 
        Fit[data12, Table[temp^i, {i, 0, 3}], 
          temp];
Plot[xAPFRvsT2, {temp, 298, 358}];
Show[%, dataPlot12];
xAPFRvsT3 = 
        Fit[data13, Table[temp^i, {i, 0, 3}], 
          temp];
Plot[xAPFRvsT3, {temp, 298, 358}];
Show[%, dataPlot13];

data21 = 
        List[{temp[1], xAMFR11}, {temp[2], xAMFR21}, {temp[3], 
            xAMFR31}, {temp[4], xAMFR41}, {temp[5], xAMFR51}, {temp[6], 
            xAMFR61}, {temp[7], xAMFR71}];
data22 = 
        List[{temp[1], xAMFR12}, {temp[2], xAMFR22}, {temp[3], 
            xAMFR32}, {temp[4], xAMFR42}, {temp[5], xAMFR52}, {temp[6], 
            xAMFR62}, {temp[7], xAMFR72}];
data23 = 
        List[{temp[1], xAMFR13}, {temp[2], xAMFR23}, {temp[3], 
            xAMFR33}, {temp[4], xAMFR43}, {temp[5], xAMFR53}, {temp[6], 
            xAMFR63}, {temp[7], xAMFR73}];
dataPlot21 = 
        ListPlot[data21, 
          PlotStyle -> PointSize[0.02]];
dataPlot22 = 
        ListPlot[data22, 
          PlotStyle -> PointSize[0.02]];
dataPlot23 = 
        ListPlot[data23, 
          PlotStyle -> PointSize[0.02]];
dataPlot21p = 
      Regress[data21, Table[temp^i, {i, 0, 3}], temp, 
        RegressionReport -> RSquared]
dataPlot22p = 
      Regress[data22, Table[temp^i, {i, 0, 3}], temp, 
        RegressionReport -> RSquared]
dataPlot23p = 
      Regress[data23, Table[temp^i, {i, 0, 3}], temp, 
        RegressionReport -> RSquared]
xAMFRvsT1 = 
        Fit[data21, Table[temp^i, {i, 0, 3}], 
          temp];
Plot[xAMFRvsT1, {temp, 298, 358}];
Show[%, dataPlot21];
xAMFRvsT2 = 
        Fit[data22, Table[temp^i, {i, 0, 3}], 
          temp];
Plot[xAMFRvsT2, {temp, 298, 358}];
Show[%, dataPlot22];
xAMFRvsT3 = 
        Fit[data23, Table[temp^i, {i, 0, 3}], 
          temp];
Plot[xAMFRvsT3, {temp, 298, 358}];
Show[%, dataPlot23];

Plot[{xAPFRvsT1, xAPFRvsT2, xAPFRvsT3, xAMFRvsT1, xAMFRvsT2, 
          xAMFRvsT3, xAe}, {temp, 268, 358}, GridLines -> Automatic, 
        PlotStyle -> {{AbsoluteThickness[1.5], 
              AbsoluteDashing[{10}]}, {AbsoluteThickness[2.5], 
              AbsoluteDashing[{10}]}, {AbsoluteThickness[3.5], 
              AbsoluteDashing[{10}]}, AbsoluteThickness[1.5], 
            AbsoluteThickness[2.5], AbsoluteThickness[3.5], 
            AbsoluteThickness[5]}];

Now compare this graph with Figure E9.3 and figure out why reaction rate is the parameter there, rather than space time (as shown here). Dashed lines correspond to plug flow (or a batch reactor); solid lines correspond to mixed flow. The thickest line is the equilibrium line; the thinnest lines correspond to the base case (10 min).

cA0 = 1;
k160 = k1 /. temp -> 333
xAe60 = xAe /. temp -> 333
k260 = k2 /. temp -> 333

Solve[0.5 == k160*cA0 - k160*cA0*xA - k260*cA0*xA, 
      xA]
Solve[0.3 == k160*cA0 - k160*cA0*xA - k260*cA0*xA, 
      xA]
Solve[0.2 == k160*cA0 - k160*cA0*xA - k260*cA0*xA, 
      xA]
Solve[0.1 == k160*cA0 - k160*cA0*xA - k260*cA0*xA, 
      xA]
Solve[0.05 == k160*cA0 - k160*cA0*xA - k260*cA0*xA, 
      xA]
Solve[0.03 == k160*cA0 - k160*cA0*xA - k260*cA0*xA, 
      xA]
Solve[0.02 == k160*cA0 - k160*cA0*xA - k260*cA0*xA, 
      xA]
Solve[0.01 == k160*cA0 - k160*cA0*xA - k260*cA0*xA, 
      xA]

The above values do correspond quite well with the values for the reaction rate at 60 °C.

k110 = k1 /. temp -> 283
xAe10 = xAe /. temp -> 283
k210 = k2 /. temp -> 283

Solve[0.03 == k110*cA0 - k110*cA0*xA - k210*cA0*xA, 
      xA]
Solve[0.02 == k110*cA0 - k110*cA0*xA - k210*cA0*xA, 
      xA]
Solve[0.01 == k110*cA0 - k110*cA0*xA - k210*cA0*xA, 
      xA]
Solve[0.005 == k110*cA0 - k110*cA0*xA - k210*cA0*xA, 
      xA]
Solve[0.003 == k110*cA0 - k110*cA0*xA - k210*cA0*xA, 
      xA]
Solve[0.002 == k110*cA0 - k110*cA0*xA - k210*cA0*xA, 
      xA]
Solve[0.001 == k110*cA0 - k110*cA0*xA - k210*cA0*xA, 
      xA]

The above values do correspond very well with the values for the reaction rate at 10 °C.

k190 = k1 /. temp -> 90 + 273
xAe90 = xAe /. temp -> 90 + 273
k290 = k2 /. temp -> 90 + 273

Solve[3.  == k190*cA0 - k190*cA0*xA - k290*cA0*xA, 
      xA]
Solve[2.  == k190*cA0 - k190*cA0*xA - k290*cA0*xA, 
      xA]
Solve[1.  == k190*cA0 - k190*cA0*xA - k290*cA0*xA, 
      xA]
Solve[0.5 == k190*cA0 - k190*cA0*xA - k290*cA0*xA, 
      xA]
Solve[0.3 == k190*cA0 - k190*cA0*xA - k290*cA0*xA, 
      xA]
Solve[0.2 == k190*cA0 - k190*cA0*xA - k290*cA0*xA, 
      xA]
Solve[0.1 == k190*cA0 - k190*cA0*xA - k290*cA0*xA, 
      xA]

And finally, the above values do correspond very well with the values for the reaction rate at 90 °C. Can you now prepare a more "elegant" solution to the construction of the graph in Figure E9.3?