Briggs-Rauscher Mechanism
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reactions involved in the BR mechanism:
2H++I-+IO3- ⇀ HOI+HIO2
H++I-+HIO2 ⇀ 2HOI
HOI+ H++I- ⇌  I2+H2O
IO3-+HIO2+H+   ⇌  2IO2+H2O
2 HIO2   ⇀   HOI+ IO3-+ H+
IO2+ H2O+Mn2+  ⇌   HIO2 +MnOH2+
MnOH2++ H2O2  ⇀  HO2 +Mn2++ H2O
2HO2 ⇀  H2O2+O2
I2+MA  ⇀  IMA+ H++I-
HOI+ H2O2  ⇀  O2+ H++I-+ H2O
where MA is Malonic Acid
Reference:
P. DeKepper and I. R. Epstein, American Chemical Society, Vol. 104, 1982, N. 1, p.49-55

Symbols in the program:
x1, x2, x3,x4, x5,x6, x7,x8, x9 and x10 in the Mathematica code are the 
compositions of  IO3- , I- , I2  , HIO2  , HOI  , IO2  , MnOH2+ ,  HO2  ,  MA  ,  H2O2  . 

Concentration of H+ is taken as constant equal to h. 

The symbols C9, Ki and ki are the different rate constants and constants 
appearing in the rate laws (see DeKepper's reference for details). 

Xi0 are the 10 concentration of the various components in the feed stream of this open reactor. 

The residence time equal to 1/k0=156 s. 

A is a constant equal to [Mn2+ ]+[MnOH2+].

Bistability:
The values of X10 and X30 are chosen in order to be in the region of the phase diagram where there is bistability.
We get two different steady states for different initial conditions:

Initial condition 1:
x1[0]==10^-2, x2[0]==10^-8, x3[0]==6 10^-7,
x4[0]==10^-10, x5[0]==10^-10, x6[0]==10^-10,
x7[0]==10^-13, x8[0]==0, x9[0]==10^-3, x10[0]== 0.

Initial Condition 2:
x1[0]==10^-2, x2[0]==10^-9, x3[0]==6 10^-4,
x4[0]==10^-6, x5[0]==10^-6, x6[0]==10^-7,
x7[0]==10^-9, x8[0]==0, x9[0]==10^-3, x10[0]== 0

Click on link to see steady states values for the different components: Steady state 1 and steady state 2.

Concentration Profiles are given by the links below: 
Component \(TraditionalForm\`\(I\^-\)\) : figure 1
Component (Subscript[I, 2]) : figure 2
Component (Subscript[HIO, 2]) : figure 3
Component HOI : figure 4
Component (Subscript[IO, 2]) : figure 5
Component \(TraditionalForm\`\(MnOH\^+\)\) : figure 6
Component MA : figure 7

X10 = 0.04; X20 = 0; X30 = 2.5 10^-6;
X40 = 0; X50 = 0; X60 = 0; X70 = 0; X80 = 0; X90 = 0.0015;  
    X100 = 0.33; h = 0.056; A = 0.004; k1 = 1.43 10^3; k2 = 2 10^10;  
    k3 = 3.1 10^12; K3 = 2.2; k4 = 7.3 10^3; K4 = 1.7 10^7; k5 = 6 10^5; 
    k6 = 1 10^4;  k7 = 3.2 10^4;  k8 = 7.5 10^5; k9 = 40;  C9 = 10^4;  
    k0 = 1/156; k10 = 37;

Here is solution for the first initial condition

\(\(\(sol\)\(=\)\(NDSolve[{
\(x1'\)[
              t] == \(-k1\)\ h^2\ x2[t]\ x1[t] - k4\ h\ x1[t]\ x4[t] + 
              K4\ x6[t]^2 + k5\ x4[t]^2 + 
              k0\ \((X10 - x1[t])\), 
\(x2'\)[
              t] == \ \(-k1\)\ h^2\ x2[t]\ x1[t] - 
              k2\ h\ x4[t]\ x2[t] - k3\ h\ x5[t]\ x2[t] + K3\ x3[t] + 
              k9\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + k10\ x10[t]\ x5[t] + 
              k0\ \((X20 - x2[t])\), 
\(x3'\)[t] == \ 
            k3\ h\ x5[t]\ x2[t] - K3\ x3[t] - 
              k9\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + 
              k0\ \((X30 - x3[t])\), 
\(x4'\)[t] == \ 
            k1\ h^2\ x2[t]\ x1[t] - k2\ h\ x4[t]\ x2[t] - 
              k4\ h\ x1[t]\ x4[t] + K4\ x6[t]^2 - 2  k5\ x4[t]^2 + 
              k6\ x6[t]\ \((A - x7[t])\) + 
              k0\ \((X40 - x4[t])\), 
\(x5'\)[t] == 
            k1\ h^2\ x2[t]\ x1[t] + 2\ k2\ h\ x4[t]\ x2[t]\  - 
              k3\ h\ x5[t]\ x2[t] + K3\ x3[t] + k5\ x4[t]^2 - 
              k10\ x10[t]\ x5[t]
 + 
              k0\ \((X50 - x5[t])\), 
\(x6'\)[t] == \ 
            2\ k4\ h\ x1[t]\ x4[t] - 2\ K4\ x6[t]^2 - 
              k6\ x6[t]\ \((A - x7[t])\) + 
              k0\ \((X60 - x6[t])\), 
\(x7'\)[t] == \ 
            k6\ x6[t]\ \((A - x7[t])\) - k7\ x7[t]\ x10[t] + 
              k0\ \((X70 - x7[t])\), 
\(x8'\)[t] == \ 
            k7\ x7[t]\ x10[t] - 2\ k8\ x8[t]^2 + 
              k0\ \((X80 - x8[t])\), 
\(x9'\)[
              t] == \ \(-k9\)\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + 
              k0\ \((X90 - x9[t])\), 
\(x10'\)[
              t] == \ \(-k7\)\ x7[t]\ x10[t] + \ k8\ x8[t]^2 - 
              k10\ x10[t]\ x5[t] + 
              k0\ \((X100 - x10[t])\), 
x1[0] == 
            10^\(-2\), x2[0] == 10^\(-8\), 
          x3[0] == 6\ 10^\(-7\), 
x4[0] == 
            10^\(-10\), x5[0] == 10^\(-10\), 
          x6[0] == 10^\(-10\), 
x7[0] == 
            10^\(-13\), x8[0] == 0, x9[0] == 10^\(-3\), 
          x10[0] == \ 0}, {x1, x2, x3, x4, x5, x6, x7, x8, x9, 
          x10}, {t, 0, 10000}, WorkingPrecision -> 15, 
        AccuracyGoal -> 10, 
\ \ \ PrecisionGoal -> 10, 
        MaxSteps -> Infinity]\)\(
\)\)\)

Plot[Evaluate[x2[t] /. sol], {t, 0, 1000}, 
      PlotRange -> All]

Plot[Evaluate[x3[t] /. sol], {t, 0, 2000}, 
      PlotRange -> All]

Plot[Evaluate[x4[t] /. sol], {t, 0, 2000}, 
      PlotRange -> All]

Plot[Evaluate[x5[t] /. sol], {t, 0, 2000}, 
      PlotRange -> All]

Plot[Evaluate[x6[t] /. sol], {t, 0, 2000}, 
      PlotRange -> All]

Plot[Evaluate[x7[t] /. sol], {t, 0, 2000}, 
      PlotRange -> All]

Plot[Evaluate[x9[t] /. sol], {t, 0, 2000}, 
      PlotRange -> All]

Steady State Values for the first initial condition

Evaluate[x2[t] /. sol] /. t -> 10000

Evaluate[x3[t] /. sol] /. t -> 10000

Evaluate[x4[t] /. sol] /. t -> 10000

Evaluate[x5[t] /. sol] /. t -> 10000

Evaluate[x6[t] /. sol] /. t -> 10000

Evaluate[x7[t] /. sol] /. t -> 10000

Evaluate[x9[t] /. sol] /. t -> 10000

Here is the solution for the second initial condition

\(sol1 = 
      NDSolve[{
\(x1'\)[
              t] == \(-k1\)\ h^2\ x2[t]\ x1[t] - k4\ h\ x1[t]\ x4[t] + 
              K4\ x6[t]^2 + k5\ x4[t]^2 + 
              k0\ \((X10 - x1[t])\), 
\(x2'\)[
              t] == \ \(-k1\)\ h^2\ x2[t]\ x1[t] - 
              k2\ h\ x4[t]\ x2[t] - k3\ h\ x5[t]\ x2[t] + K3\ x3[t] + 
              k9\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + k10\ x10[t]\ x5[t] + 
              k0\ \((X20 - x2[t])\), 
\(x3'\)[t] == \ 
            k3\ h\ x5[t]\ x2[t] - K3\ x3[t] - 
              k9\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + 
              k0\ \((X30 - x3[t])\), 
\(x4'\)[t] == \ 
            k1\ h^2\ x2[t]\ x1[t] - k2\ h\ x4[t]\ x2[t] - 
              k4\ h\ x1[t]\ x4[t] + K4\ x6[t]^2 - 2  k5\ x4[t]^2 + 
              k6\ x6[t]\ \((A - x7[t])\) + 
              k0\ \((X40 - x4[t])\), 
\(x5'\)[t] == 
            k1\ h^2\ x2[t]\ x1[t] + 2\ k2\ h\ x4[t]\ x2[t]\  - 
              k3\ h\ x5[t]\ x2[t] + K3\ x3[t] + k5\ x4[t]^2 - 
              k10\ x10[t]\ x5[t]
 + 
              k0\ \((X50 - x5[t])\), 
\(x6'\)[t] == \ 
            2\ k4\ h\ x1[t]\ x4[t] - 2\ K4\ x6[t]^2 - 
              k6\ x6[t]\ \((A - x7[t])\) + 
              k0\ \((X60 - x6[t])\), 
\(x7'\)[t] == \ 
            k6\ x6[t]\ \((A - x7[t])\) - k7\ x7[t]\ x10[t] + 
              k0\ \((X70 - x7[t])\), 
\(x8'\)[t] == \ 
            k7\ x7[t]\ x10[t] - 2\ k8\ x8[t]^2 + 
              k0\ \((X80 - x8[t])\), 
\(x9'\)[
              t] == \ \(-k9\)\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + 
              k0\ \((X90 - x9[t])\), 
\(x10'\)[
              t] == \ \(-k7\)\ x7[t]\ x10[t] + \ k8\ x8[t]^2 - 
              k10\ x10[t]\ x5[t] + 
              k0\ \((X100 - x10[t])\), 
x1[0] == 
            10^\(-2\), x2[0] == 10^\(-9\), 
          x3[0] == 6\ 10^\(-4\), 
x4[0] == 
            10^\(-6\), x5[0] == 10^\(-6\), 
          x6[0] == 10^\(-7\), 
x7[0] == 
            10^\(-9\), x8[0] == 0, x9[0] == 10^\(-3\), 
          x10[0] == \ 0}, {x1, x2, x3, x4, x5, x6, x7, x8, x9, 
          x10}, {t, 0, 10000}, WorkingPrecision -> 15, 
        AccuracyGoal -> 10, 
\ \ \ PrecisionGoal -> 10, 
        MaxSteps -> Infinity]\)

Plot[Evaluate[x2[t] /. sol1], {t, 0, 1000}, 
      PlotRange -> All]

Plot[Evaluate[x3[t] /. sol1], {t, 0, 1000}, 
      PlotRange -> All]

Plot[Evaluate[x4[t] /. sol1], {t, 0, 1000}, 
      PlotRange -> All]

Plot[Evaluate[x5[t] /. sol1], {t, 0, 1000}, 
      PlotRange -> All]

Plot[Evaluate[x6[t] /. sol1], {t, 0, 1000}, 
      PlotRange -> All]

Plot[Evaluate[x7[t] /. sol1], {t, 0, 1000}, 
      PlotRange -> All]

Plot[Evaluate[x9[t] /. sol1], {t, 0, 1000}, 
      PlotRange -> All]

Steady State Values for the second initial condition

Evaluate[x2[t] /. sol1] /. t -> 10000

Evaluate[x3[t] /. sol1] /. t -> 10000

Evaluate[x4[t] /. sol1] /. t -> 10000

Evaluate[x5[t] /. sol1] /. t -> 10000

Evaluate[x6[t] /. sol1] /. t -> 10000

Evaluate[x7[t] /. sol1] /. t -> 10000

Evaluate[x9[t] /. sol1] /. t -> 10000

Comparison of the concentration profiles for the first and second initial condition.

Plot[Evaluate[{x2[t] /. sol, x2[t] /. sol1}], {t, 0, 1000}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

Plot[Evaluate[{1000 x3[t] /. sol, x3[t] /. sol1}], {t, 0, 1000}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

Plot[Evaluate[{1000 x4[t] /. sol, x4[t] /. sol1}], {t, 0, 2000}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

Plot[Evaluate[{1000 x5[t] /. sol, x5[t] /. sol1}], {t, 0, 2000}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

Plot[Evaluate[{1000 x6[t] /. sol, x6[t] /. sol1}], {t, 0, 2000}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

Plot[Evaluate[{1000 x7[t] /. sol, x7[t] /. sol1}], {t, 0, 2000}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

Plot[Evaluate[{x9[t] /. sol, x9[t] /. sol1}], {t, 0, 2000}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], 
          RGBColor[1, 0, 0]}, 
AxesLabel -> {t, 
          x9}]

Getting Oscillations in the composition profile and a limit cyle:
Here we change the values of X10 and X30 to move into the region of the phase diagram where the concentrations of the components will oscillate.
Click to see Oscillations.
Click to see Limit Cycle.

X10 = 0.035; X20 = 0; X30 = 10^-6;
X40 = 0; X50 = 0; X60 = 0; X70 = 0; X80 = 0; X90 = 0.0015;  
    X100 = 0.33; h = 0.056; A = 0.004; k1 = 1.43 10^3; k2 = 2 10^10;  
    k3 = 3.1 10^12; K3 = 2.2; k4 = 7.3 10^3; K4 = 1.7 10^7; k5 = 6 10^5; 
    k6 = 1 10^4;  k7 = 3.2 10^4;  k8 = 7.5 10^5; k9 = 40;  C9 = 10^4;  
    k0 = 1/156; k10 = 37;

\(sol2 = 
      NDSolve[{
\(x1'\)[
              t] == \(-k1\)\ h^2\ x2[t]\ x1[t] - k4\ h\ x1[t]\ x4[t] + 
              K4\ x6[t]^2 + k5\ x4[t]^2 + 
              k0\ \((X10 - x1[t])\), 
\(x2'\)[
              t] == \ \(-k1\)\ h^2\ x2[t]\ x1[t] - 
              k2\ h\ x4[t]\ x2[t] - k3\ h\ x5[t]\ x2[t] + K3\ x3[t] + 
              k9\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + k10\ x10[t]\ x5[t] + 
              k0\ \((X20 - x2[t])\), 
\(x3'\)[t] == \ 
            k3\ h\ x5[t]\ x2[t] - K3\ x3[t] - 
              k9\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + 
              k0\ \((X30 - x3[t])\), 
\(x4'\)[t] == \ 
            k1\ h^2\ x2[t]\ x1[t] - k2\ h\ x4[t]\ x2[t] - 
              k4\ h\ x1[t]\ x4[t] + K4\ x6[t]^2 - 2  k5\ x4[t]^2 + 
              k6\ x6[t]\ \((A - x7[t])\) + 
              k0\ \((X40 - x4[t])\), 
\(x5'\)[t] == 
            k1\ h^2\ x2[t]\ x1[t] + 2\ k2\ h\ x4[t]\ x2[t]\  - 
              k3\ h\ x5[t]\ x2[t] + K3\ x3[t] + k5\ x4[t]^2 - 
              k10\ x10[t]\ x5[t]
 + 
              k0\ \((X50 - x5[t])\), 
\(x6'\)[t] == \ 
            2\ k4\ h\ x1[t]\ x4[t] - 2\ K4\ x6[t]^2 - 
              k6\ x6[t]\ \((A - x7[t])\) + 
              k0\ \((X60 - x6[t])\), 
\(x7'\)[t] == \ 
            k6\ x6[t]\ \((A - x7[t])\) - k7\ x7[t]\ x10[t] + 
              k0\ \((X70 - x7[t])\), 
\(x8'\)[t] == \ 
            k7\ x7[t]\ x10[t] - 2\ k8\ x8[t]^2 + 
              k0\ \((X80 - x8[t])\), 
\(x9'\)[
              t] == \ \(-k9\)\ x9[t]\ x3[t]/\((1 + C9\ x3[t])\) + 
              k0\ \((X90 - x9[t])\), 
\(x10'\)[
              t] == \ \(-k7\)\ x7[t]\ x10[t] + \ k8\ x8[t]^2 - 
              k10\ x10[t]\ x5[t] + 
              k0\ \((X100 - x10[t])\), 
x1[0] == 
            10^\(-2\), x2[0] == 10^\(-8\), 
          x3[0] == 6\ 10^\(-7\), 
x4[0] == 
            10^\(-10\), x5[0] == 10^\(-10\), 
          x6[0] == 10^\(-10\), 
x7[0] == 
            10^\(-13\), x8[0] == 0, x9[0] == 10^\(-3\), 
          x10[0] == \ 0}, {x1, x2, x3, x4, x5, x6, x7, x8, x9, 
          x10}, {t, 0, 10000}, WorkingPrecision -> 15, 
        AccuracyGoal -> 10, 
\ \ \ PrecisionGoal -> 10, 
        MaxSteps -> Infinity]\)

Plot[Evaluate[Log[10, x3[t]] /. sol2], {t, 0, 2000}, 
      PlotRange -> All]

Plot[Evaluate[Log[10, x2[t]] /. sol2], {t, 0, 2000}, 
      PlotRange -> All]

Plot[Evaluate[{Log[10, x5[t]], Log[10, x4[t]]} /. sol2], {t, 0, 4000}, 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

ParametricPlot[
      Evaluate[{-Log[10, x2[t]], 10000 x3[t]} /. sol2], {t, 0, 2000}, 
      PlotRange -> {{5, 12}, {-1, 8}}, 
      PlotStyle -> RGBColor[1, 0, 1]]

Plot[
      Evaluate[{Log[10, x5[t]], Log[10, x4[t]], Log[10, x2[t]], 
            Log[10, x3[t]]} /. sol2], {t, 0, 1000}, 
      PlotRange -> {-10, -3}, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0], 
          RGBColor[0, 1, 0], RGBColor[0, 0, 0]}]