Dynamic of a Diabatic Continuous 
Stirred Tank Reactor
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :
Numerical values were obtained from Process Dynamics, Modeling, Analysis and Simulation, B. Wayne Bequette, Prentice Hall, 1998.

Dimensionless Governing Equations and dimensionless groups

Numerical values

δ=150/500.;

γ=11834/(298.15 1.987);

(Subscript[D, a]) = 9703 3600 Exp[-γ];

B=5960 10/500 1/298.15 γ;

(Subscript[x, 1  c]) = 0;

Multiple steady states are found for the numerical values of γ, B, δ, (Subscript[x, 1  c]) and (Subscript[D, a]). pt1 is the low temperature and low conversion stable steady state. pt2 is the high temperature and high conversion stable steady state. pt3 is the intermediate temperature and intermediate conversion unstable steady state

root1 = 
      FindRoot[{-a1 - δ (a1 - (Subscript[x, 1  c])) + 
              B (Subscript[D, a]) (1 - a2) Exp[a1/(1 + a1/γ)] == 
            0, -a2 + 
              (Subscript[D, a]) (1 - a2) Exp[a1/(1 + a1/γ)] == 
            0}, {a1, 0.5}, {a2, 0.5}]

pt1={a2,a1}/.root1

root2 = 
      FindRoot[{-a1 - δ (a1 - (Subscript[x, 1  c])) + 
              B (Subscript[D, a]) (1 - a2) Exp[a1/(1 + a1/γ)] == 
            0, -a2 + 
              (Subscript[D, a]) (1 - a2) Exp[a1/(1 + a1/γ)] == 
            0}, {a1, 5}, {a2, 0.8}]

pt2={a2,a1}/.root2

root3 = 
      FindRoot[{-a1 - δ (a1 - (Subscript[x, 1  c])) + 
              B (Subscript[D, a]) (1 - a2) Exp[a1/(1 + a1/γ)] == 
            0, -a2 + 
              (Subscript[D, a]) (1 - a2) Exp[a1/(1 + a1/γ)] == 
            0}, {a1, 2}, {a2, 0.3}]

pt3={a2,a1}/.root3

Dynamic behavior depends on initial condition. Only steady states one and two are stable. Steady state three is a saddle point.

sys[x10_, 
        x20_] := {D[
            (Subscript[x, 1])[t], {t, 
              1}] == -(Subscript[x, 1])[
                t] - δ ((Subscript[x, 1])[t] - (Subscript[x, 1  c])) + 
            B (Subscript[D, a]) (1 - (Subscript[x, 2])[t]) Exp[
                (Subscript[x, 1])[t]/(1 + (Subscript[x, 1])[t]/γ)], 
        D[(Subscript[x, 2])[t], {t, 1}] == -(Subscript[x, 2])[t] + 
            (Subscript[D, a]) (1 - (Subscript[x, 2])[t]) Exp[(Subscript[x, 1])[t]/(1 + (Subscript[x, 1])[t]/γ)], 
        (Subscript[x, 1])[0] == x10, (Subscript[x, 2])[0] == x20}

sol[a_, b_] := 
      NDSolve[sys[x10, x20] /. x10 -> a /. x20 -> b, {(Subscript[x, 1])[t], 
          (Subscript[x, 2])[t]}, {t, 0, 1000}]

Steady State One: Temperature and Conversion versus time

plt1 = 
      ParametricPlot[Evaluate[{(Subscript[x, 2])[t], (Subscript[x, 1])[t]} /. sol[0, 0]], {t, 0, 1000}, 
        PlotRange -> All]

PLT1 = 
      Plot[Evaluate[298.15/19.9756 (Subscript[x, 1])[t] + 298.15 /. sol[0, 0]], {t, 0, 
          100}, PlotRange -> All, 
        PlotStyle -> RGBColor[1, 0, 0]]

PLT2 = 
      Plot[Evaluate[(Subscript[x, 2])[t] /. sol[0, 0]], {t, 0, 100}, 
        PlotRange -> All, 
        PlotStyle -> RGBColor[1, 0, 0]]

Steady state conversion

(Subscript[x, 2])[t] /. sol[0, 0] /. t -> 100

Steady state temperature

298.15/19.97 (Subscript[x, 1])[t] + 298.15 /. sol[0, 0] /. t -> 100

Steady State Two: Temperature and Conversion versus time.

plt2 = 
      ParametricPlot[
        Evaluate[{(Subscript[x, 2])[t], (Subscript[x, 1])[t]} /. sol[4, 0.8]], {t, 0, 100}, 
        PlotRange -> All]

PLT3 = 
      Plot[Evaluate[298.15/19.9756 (Subscript[x, 1])[t] + 298.15 /. sol[4, 0.8]], {t, 0, 
          100}, PlotRange -> All, 
        PlotStyle -> RGBColor[0, 1, 0]]

PLT4 = 
      Plot[Evaluate[(Subscript[x, 2])[t] /. sol[4, 0.8]], {t, 0, 100}, 
        PlotRange -> All, 
        PlotStyle -> RGBColor[0, 1, 0]]

Steady state conversion

(Subscript[x, 2])[t] /. sol[4, 0.8] /. t -> 100

Steady state temperature

298.15/19.97 (Subscript[x, 1])[t] + 298.15 /. sol[4, 0.8] /. 
      t -> 100

Showing both stable steady state temperature in the same plots

Show[PLT1,PLT3]

Showing both stable steady state conversion in the same plots

Show[PLT2,PLT4]

Phase plane plot. Steady state one, two and three are the red, green and blue dots respectively. One can see clearly that steady states one and two are stable while steady state three is a saddle point. Temperature overshoot is observed. Knowing about possible overshoot behavior is a valuable information since temperature overshoot can lead to safety problems.

For[i = 0, i < 100, i += 1, 
      plt[i] = ParametricPlot[
          Evaluate[{(Subscript[x, 2])[t], (Subscript[x, 1])[t]} /. sol[ i 0.1, 1]], {t, 0, 10}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

For[i = 0, i < 100, i += 1, 
      PLT[i] = ParametricPlot[
          Evaluate[{(Subscript[x, 2])[t], (Subscript[x, 1])[t]} /. sol[i 0.1, 0]], {t, 0, 10}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

For[i = 0, i < 100, i += 1, 
      plot[i] = 
        ParametricPlot[
          Evaluate[{(Subscript[x, 2])[t], (Subscript[x, 1])[t]} /. sol[0, 0.01 i]], {t, 0, 10}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

tbl=Table[plt[i],{i,0,99,3}];

TBL=Table[PLT[i],{i,0,99,2}];

TbL=Table[plot[i],{i,0,99,4}];

Show[Flatten[Join[TBL,tbl,TbL]],Graphics[{{PointSize[0.02],RGBColor[1,0,0],Point[pt1]},{PointSize[0.02],RGBColor[0,1,0],Point[pt2]},{PointSize[0.02],RGBColor[0,0,1],Point[pt3]}}],PlotRange->{{0,1},{0,6}},DisplayFunction->$DisplayFunction]