=Paper= {{Paper |id=Vol-1638/Paper87 |storemode=property |title=Study of oscillatory processes in the one model of electrochemical reactor |pdfUrl=https://ceur-ws.org/Vol-1638/Paper87.pdf |volume=Vol-1638 |authors=Natalia Firstova,Elena Shchepakina }} ==Study of oscillatory processes in the one model of electrochemical reactor == https://ceur-ws.org/Vol-1638/Paper87.pdf
Mathematical Modelling




 STUDY OF OSCILLATORY PROCESSES IN
THE ONE MODEL OF ELECTROCHEMICAL
             REACTOR

                         N. Firstova, E. Shchepakina
              Samara National Research University, Samara, Russia


     Abstract. The paper deals with the study of the features of oscilla-
     tory processes in the electrocatalytic reaction mechanism underlying
     the electrochemical reactor. The conditions for the existence of several
     types of oscillations (small amplitude harmonic oscillations, relaxation
     oscillations and oscillations of transitional type) is obtained.

     Keywords: singular perturbations, critical phenomena, canard,
     Andronov-Hopf bifurcation, electrocatalytic reaction

     Citation: Firstova N, Shchepakina E. Study of oscillatory processes
     in the one model of electrochemical reactor. CEUR Workshop Pro-
     ceedings, 2016; 1638: 731-741. DOI: 10.18287/1613-0073-2016-
     1638-731-741

Introduction
The study of the nature and patterns of occurrence of oscillations in dynamic sys-
tems are traditionally of interest [1]. Particular attention is paid to the study of
self-oscillatory regimes, which usually are unfavorable for the functioning of the
technical systems. Hence we should be able to predict and suppress such modes.
On the other hand, there is a number of modern technologies, such as chemical
industry and the energy sector, requiring the creation of the reactors substan-
tially non-equilibrium regimes, including oscillatory and, on the contrary, there
is the need to generate and control the self-oscillatory regimes. In this paper
we investigate several types of oscillations in the one model of electrochemical
reactor [2]: the oscillations of relaxation type, small oscillations (when the stable
limit cycle is small), and the critical oscillations (the canards [3, 4]).

Mathematical model of the Koper-Sluyters electrocatalytic
reaction
The Koper-Sluyters electrocatalytic reaction (KS–reaction) is a chemical reac-
tion corresponding to the following kinetic scheme (written by the symbols used
by the authors [2]):

        D/δ       ka       k
                         e
 Xbulk −−−→ Xsur ­ Xads −→ P + ne− .
                  kd

Here X is the single species, which diffuses towards the electrode where it suc-
cessively adsorbs and is electrochemically oxidized; D is the diffusion coefficient



Information Technology and Nanotechnology (ITNT-2016)                            731
Mathematical Modelling                                   Firstova N, Shchepakina E. Study...


of X; δ is the thickness of the Nernst diffusion layer; ka , ke , kd are the rate
constants for adsorption, desorption and electron transfer, respectively. The ox-
idation products P are assumed not to be adsorbed and to leave neighborhood
of the interface.
The mathematical model of the KS–reaction in dimensionless form is [2]

 du
    = −ka eγθ/2 u(1 − θ) + kd e−γθ/2 θ + 1 − u = f (u, θ),                               (1)
 dt

     dθ
 β      = ka eγθ/2 u(1 − θ) − kd e−γθ/2 θ − ke eα0 f E θ = g(u, θ),                      (2)
     dt
where u is the dimensionless interfacial concentration of X; θ is the dimension-
less amount of X that is adsorbed on the electrode surface; E is the electrode
potential; β is the coverage ratio of the adsorbate; α0 is the symmetry factor
for the electron transfer; and f = F/(RT ), where R, F and T have their usual
meaning [5]. The physical meaning of the parameter γ has always been a subject
of dispute. In most of the literature it is interpreted as an interaction param-
eter. Positive γ signifies attractive and negative γ signifies repulsive adsorbate
interactions.
Since the parameter β is small, the system (1), (2) is singularly perturbed. We
will investigate the dynamics of the solutions depending on the values of the
additional parameters of the system (1), (2). The study will be carried out
by using methods of the theory of singular perturbations [6, 7] and numerical
methods.
One of the major problems of chemical kinetics is a problem of management of
chemical processes. The goal of the control of the KS–reaction is the obtaining
a relatively high value of the reactant concentration in the framework of a safe
process. It has been found that this goal is realized during the critical mode.
The determination of the conditions of occurrence of a critical regime in the
considered chemical system is the main objective of this study.

Analysis of the slow curve
The equation of the slow curve [4,6] of (1), (2) is determined from the expression
g(u, θ) = 0 and has the form

        (kd e−γθ/2 + ke eα0 f E )θ
 u=                                .                                                     (3)
             ka eγθ/2 (1 − θ)

To study the stable and unstable parts we should find the jump points [4, 6] on
the slow curve (3).
The coordinates of the jump points are determined by the system:
 (
  g(u, θ) = 0,
   ∂g(u,θ)                                                                  (4)
     ∂θ    = 0.

For (1), (2) the system (4) has the form
(
  ka eγθ/2 u(1 − θ) − kd e−γθ/2 θ − ke eα0 f E θ = 0,
                                                                                         (5)
  ka u(1 − θ) γ2 eγθ/2 − ka ueγθ/2 − kd e−γθ/2 + kd e−γθ/2 θ γ2 − ke eα0 f E = 0.



Information Technology and Nanotechnology (ITNT-2016)                                   732
Mathematical Modelling                                  Firstova N, Shchepakina E. Study...


The system (5) is transcendental, and therefore, it is impossible to find the
solutions of the system in an analytical form, but with specific values of param-
eters, this system is solved numerically [8]. The jump points separate the stable
(attractive) and unstable (repelling) parts of the slow curve (3).
The shape of the slow curve depends on the value of the parameters. Depending
on the ratio of the parameters the three following cases are possible.
In the first case, 0 < γ < 4, the system (5) has no solutions. Therefore, the slow
manifold is either entirely stable or entirely unstable. Due to the fact that
  ∂g(u, θ)
           < 0,
      ∂θ
i.e.,
             γ                                          γ
 ka u(1 − θ) eγθ/2 − ka ueγθ/2 − kd e−γθ/2 + kd e−γθ/2 θ − ke eα0 f E < 0,
             2                                          2
the slow curve is stable in this case. Hence, the trajectories of the system (1),
(2) are attracted to the slow curve and then follow along it.
Let us consider the second case when γ ≈ 4. In this case, the system (5) has one
solution. The derivative ∂g(u, θ)/∂θ in the transition through this point does
not change its sign. So, we should find the inflection points. The coordinates of
the inflection points are determined by the system
(
 g(u, θ) = 0,
  ∂ 2 g(u,θ)                                                                            (6)
      ∂θ 2   = 0.
The system (6) in this case has a unique solution, and the slow curve is stable
in this case.
In the last case, when γ > 4, the system (5) has two solutions. The shapes of
the slow curve have the form as shown in Fig. 1.
Thus, for γ > 4 the jump points divide the slow curve into three parts, which are
zeroth approximations of the corresponding integral manifolds: near the stable
branches F1 and F3 there are the stable slow integral manifolds M1 and M3 ,
respectively; near the unstable branch F2 there is the unstable slow integral
manifold M2 .
A system’s trajectory, starting from an initial point in the basin of attraction of
the stable slow integral manifold M1 (or M3 ), will be attracted to it with the
velocity of the fast variable of order O( β1 ) as β → 0 and then follows along it
with the velocity of the slow variable, of the order O(1) as β → 0. The further
behavior of the trajectory will depend on the location of the critical point of the
system (1), (2).

Critical point
The critical points of the system (1), (2) are determined by the system
(
  g(u, θ) = 0,
  f (u, θ) = 0,
which for the system (1), (2) has the form:
(
  −ka eγθ/2 u(1 − θ) + kd e−γθ/2 θ + 1 − u = 0,
                                                                                        (7)
  ka eγθ/2 u(1 − θ) − kd e−γθ/2 θ − ke eαo f E θ = 0.



Information Technology and Nanotechnology (ITNT-2016)                                  733
Mathematical Modelling                                   Firstova N, Shchepakina E. Study...




 Fig. 1. The slow curve of the system (1), (2) for (top, left) 0 < γ < 4, (top, right)
                              γ ≈ 4, (bottom) γ > 4

From the system (7) we get:

 u∗ = 1 − ke eαo f E θ.                                                                  (8)

Substituting (8) into the second equation of the system (7), we obtain the equa-
tion that determines the value θ = θ∗ :

 ka eγθ/2 (1 − ke eαo f E θ)(1 − θ) − kd e−γθ/2 θ − ke eαo f E θ = 0.                    (9)

Thus, we obtain the critical point
  ¡                     ¢
A θ∗ , 1 − ke eαo f E θ∗ ,

where θ∗ is the solution of the equation (9).
The Jacobian matrix of the system (1), (2):
à                !
   ∂f (u,θ)   ∂f (u,θ)
     ∂u          ∂θ       ,
   ∂g(u,θ)    ∂g(u,θ)
     ∂u          ∂θ

has the characteristic equation

 λ2 + λξ1 + ξ2 = 0,

with the discriminant

 D = ξ12 − 4ξ2 ,

where
        ka γθ∗ /2                         γ            kd   ∗        γθ∗
 ξ1 =     e       (1 − ke eα0 f E θ∗ )(1 − (1 − θ∗ )) + e−γθ /2 (1 −     )
        β                                 2            β              2


Information Technology and Nanotechnology (ITNT-2016)                                   734
Mathematical Modelling                                  Firstova N, Shchepakina E. Study...

     ke α0 f E         ∗
 +     e       + ka eγθ /2 (1 − θ∗ ) + 1,
     β
        ka γθ∗ /2                         γ            kd   ∗        γθ∗
 ξ2 =     e       (1 − ke eα0 f E θ∗ )(1 − (1 − θ∗ )) + e−γθ /2 (1 −     )
        β                                 2            β              2
     ke α0 f E ka ke γθ∗ /2
 +     e      +     e       (1 − θ∗ )eα0 f E .
     β          β
The type of critical point and its coordinates depends on the value of the pa-
rameter γ. Let us consider the most interesting case when γ > 4. Without loss
of generality the parameters of the system are chosen to be ² = 0.2, γ = 8.99,
ka = 10, kd = 100, α0 = 0.05, f = 38, 7, E = 0.207564 unless other values are
specified in figure captions.
In [8] it has been shown that the critical point is a stable focus when it lies on
the stable part of the slow curve (see Fig. 2) and it is an unstable focus when
it lies on F2 . In the second case the relaxation oscillations are observed in the
system, see Fig. 3.




 Fig. 2. The slow curve (red line) and the trajectory (black line) of the system (1),
                                    (2); ke = 0.8




            Fig. 3. Relaxation oscillations in the system (1), (2); ke = 1.285


The transition between these two situations corresponds to the case when the
critical point coincides with the jump point, the stable equilibrium of the system
becomes unstable, and at the same instant the stable limit cycle is originated,
i.e., the Andronov–Hopf bifurcation occurs, see Fig.4. With further minor mod-
ifications of the control parameter, say ke , (other parameters are fixed), the
critical point moves on the unstable part of the slow curve, staying in small



Information Technology and Nanotechnology (ITNT-2016)                                  735
Mathematical Modelling                                  Firstova N, Shchepakina E. Study...


(of order O(β) as β → 0) neighbourhood of the jump point. As parameter ke
changes further this limit cycle grows, and at a value ke = ke∗ (so-called canard
point) it becomes the canard cycle (see Fig. 5) with the following it canard
explosion [3, 9–11]. Recall that the trajectories which at first move along the
stable slow integral manifold and then continue for a while along the unstable
slow integral manifold are called canards [4, 6].
From the first sight the threshold in the qualitative behaviour of the solutions of
the system corresponds to the Andronov–Hopf bifurcation point. However, when
the value of the control parameter is close to the Andronov–Hopf bifurcation
point, the size of the limit cycle is so small that the behavior of the system’
solution is practically indistinguishable from the slow mode. If, in the case of slow
regime, the trajectories approach the stable equilibrium, practically coinciding
with the jump point, in the later case they tend to a small limit cycle, nearly
coinciding with the same jump point. And only when the control parameter
attains the canard point, provided the equilibrium is on the unstable part of the
slow curve, but in the sufficiently small vicinity of the jump point, the qualitative
change in the system’s behavior can be observed. Namely, the growth of the limit
cycle occurs in such a way that it becomes possible to speak of the existence of
the canard trajectory. In other words, the appreciable change in size and/or in
form of the limit cycle is observed for small variations of the control parameter,
i.e. the canard explosion takes place. Thus, the canard point is the critical value
of the control parameter.




        Fig. 4. Andronov-Hopf bifurcation in the system (1), (2); ke = 0.92




Fig. 5. The slow curve (red line) and the canard (black line) of the system (1), (2);
                                   ke = 0.920529




Information Technology and Nanotechnology (ITNT-2016)                                  736
Mathematical Modelling                                  Firstova N, Shchepakina E. Study...


The Andronov-Hopf bifurcation in the model
The sufficient conditions for the Andronov–Hopf bifurcation occurring in the
system (1), (2) are [8]:

 ξ1 = 0,                                                                              (10)

 ξ12 − 4ξ2 < 0.                                                                       (11)

Taking into account (10) we can rewrite (11) as

 ξ2 > 0.                                                                              (12)

The expressions (10) and (12) in more detailed form are
                                   µ        ¶
 ka γθ∗ /2 ∗ ³ γ         ´ k
                            d   ∗       γθ∗    ke
   e      u 1 − (1 − θ∗ ) + e−γθ /2 1 −       + eα0 f E
 β             2           β             2     β

               ∗
 = −ka eγθ /2 (1 − θ∗ ) − 1,                                                          (13)
                                                            µ       ¶
 ka γθ∗ /2 ¡         α0 f E ∗
                              ¢³    γ      ∗
                                             ´ k
                                                 d −γθ ∗ /2     γθ∗
   e        1 − ke e       θ     1 − (1 − θ ) + e            1−
 β                                  2          β                 2

     ke α0 f E ka ke γθ∗ /2
 +     e      +     e       (1 − θ∗ )eα0 f E > 0.                                     (14)
     β          β
Substituting (13) in (14) yields
           ∗                            ∗
 ka ke eγθ /2 (1 − θ∗ )eα0 f E − β(ka eγθ /2 (1 − θ∗ ) + 1) > 0.                      (15)

Thus, the expressions (13) and (15) give us the sufficient condition for the
Andronov-Hopf bifurcation in the system under consideration. It should be
noted that for β → 0 the condition (15) is fulfilled for all values of parameters.
In the system (1), (2) the Andronov–Hopf bifurcation can occur either in the
neighborhood of jump point A1 or A2 . The Fig. 6 demonstrates the stable limit
cycle of the system (1), (2) arising via the Andronov–Hopf bifurcation in the
neighborhood of jump point A1 .




Fig. 6. (left) The slow curve (red line) and the limit cycle (black line) of the system
   (1), (2) and (right) plots of u(t) (blue line), θ(t) (yellow line) for ke = 2.40855




Information Technology and Nanotechnology (ITNT-2016)                                  737
Mathematical Modelling                                         Firstova N, Shchepakina E. Study...


Canards
The canards and the parameter value ke∗ allow asymptotic expansions in powers
of the small parameter β [4, 6, 11]:

 u = Φ(θ, β) = u0 (θ) + βu1 (θ) + β 2 u2 (θ) + . . . ,                                       (16)

 ke∗ = χ(β) = χ0 + βχ1 + β 2 χ2 + . . . .                                                    (17)
In order to find these asymptotic expansions for the canard and the canard point
we substitute the formal expansions (16) and (17) into the invariance equation [6]
 du
    g(u, θ) = βf (u, θ).
 dθ
which follows from the system (1), (2). As a result we obtain the following
equation:
 ³
   ka eγθ/2 (1 − θ)(u0 + βu1 + β 2 u2 + . . . ) − (χ0 + βχ1 + β 2 χ2 + . . . )eα0 f E θ
           ´
−kd eγθ/2 θ (u00 + βu01 + β 2 u02 + . . . ) = −βka eγθ/2 (1 − θ)(u0 + βu1 + β 2 u2 + . . . )
    ³                                            ´
+β kd eγθ/2 θ + 1 − u0 − βu1 − β 2 u2 + . . . .                                      (18)

On setting equal the coefficients of powers of β in the equation (18) we find
the functions u0 (θ), u1 (θ), . . . . To obtain the values χ0 , χ1 , . . . we require
the continuity of the functions ui (θ) (i = 0, 1, ...) at the jump point. This
requirement means that we glue the stable and the unstable integral manifolds
at the jump point and, as a result, construct the canard passing through this
point [4, 6, 9]. As a result we have:

            (kd e−γθ/2 + χ0 eα0 f E )θ
 u0 (θ) =                              ,                                                     (19)
                 ka eγθ/2 (1 − θ)


            −ka u0 (θ)(1 − θ)eγθ/2 + kd e−γθ/2 θ + 1 − u0 (θ) + χ1 eα0 f E θu00 (θ)
 u1 (θ) =                                                                           , (20)
                                      ka eγθ/2 u00 (θ)


        ka (1 − θ̄)eγ θ̄/2 − kd e−γ θ̄/2 θ̄
 χ0 =                                       ,                                                (21)
        (ka (1 − θ̄)eγ θ̄/2 − 1)eα0 f E θ̄


          ka u1 (θ̄)(1 − θ̄)eγ θ̄/2 + u1 (θ̄) + ka u1 (θ̄)u01 (θ̄)(1 − θ̄)eγ θ̄/2
 χ1 = −                                                                           ,          (22)
                                     eα0 f E θ̄u01 (θ̄)

where the value θ = θ̄ corresponding to the jump point can be calculate from
the system (4). The equations (19)–(22) define the first–order approximations
for the canard and the canard point of the system (1), (2) in a neighborhood of
the jump point (u(θ̄), θ̄). It should be noted that we can construct the canard
either in the neighborhood of jump point A1 (by gluing the stable slow integral
manifold M1 and the unstable one M2 , see Fig. 7) or in the neighborhood of
A2 (by gluing the stable slow integral manifold M3 and M2 , see Fig. 8). If it is



Information Technology and Nanotechnology (ITNT-2016)                                         738
Mathematical Modelling                                  Firstova N, Shchepakina E. Study...


necessary to glue stable and unstable slow invariant manifolds at the both jump
points simultaneously, we should use two control parameters and as a result we
obtain a canard cascade [12].




Fig. 7. (left) The slow curve (red line) and the canard (black line) of the system (1),
      (2) and (right) plots of u(t) (blue line), θ(t) (yellow line) for ke = 2.4055




Fig. 8. (left) The slow curve (red line) and the canard (black line) of the system (1),
    (2) and (right) plots of u(t) (blue line), θ(t) (yellow line) for ke = 0.9205283


Conclusion
In the paper the dynamical model of the electrochemical reactor has been inves-
tigated. The critical regime separating the basic types of the regimes, slow and
relaxation, was modelled with the help of the integral manifolds of variable sta-
bility. This approach was used in [13–23] for modelling of the critical phenomena
in chemical systems.
The bifurcation point at which the supercritical Andronov–Hopf bifurcation
takes place as well as the canard point of the control parameter at which the
system has the canard cycle have been determined analytically. It is shown that
the critical mode is modelled by the canard. The obtained results is of utmost
importance for several applications in chemical kinetics, as they can be used to
determine the dynamics of the process in the chemical system for given initial
conditions.

Acknowledgements
This work is supported in part by the Russian Foundation for Basic Research
(grant 14-01-97018-p) and the Ministry of Education and Science of the Rus-
sian Federation under the Competitiveness Enhancement Program of Samara
University (2013-2020).



Information Technology and Nanotechnology (ITNT-2016)                                  739
Mathematical Modelling                                    Firstova N, Shchepakina E. Study...


References
   1. Gol’dshtein VM, Sobolev VA, Yablonskii GS. Relaxation self–oscillations in chemical
      kinetics: a model, conditions for realization. Chem. Eng. Sci., 1986; 41: 2761–2766.
   2. Koper MTM, Sluyters JH. Instabilities and oscillations in simple models of
      electrocatalytic surface reactions. Journal of Electroanalytical Chemistry, 1994; 371(1):
      149–159.
   3. Benoit E, Calot JL, Diener F, Diener M. Chasse au canard [The duck shooting]. Collec-
      tanea Mathematica, 1981; 31–32: 37–119. [in French]
   4. Shchepakina E, Sobolev V. Black swans and canards in laser and combustion models.
      Chapter 8 (51 p.) in “Singular Perturbations and Hysteresis”, ed. by M. Mortell, R.
      O’Malley, A. Pokrovskii, V. Sobolev, Philadelphia: SIAM, 2005.
   5. Berthier F, Diard JP, Nugues S. On the nature of the spontaneous oscillations observed
      for the Koper–Sluyters electrocatalytic reaction. Journal of Electroanalytical Chemistry,
      1997; 436(1): 35–42.
   6. Shchepakina E, Sobolev V, Mortell MP. Singular Perturbations: Introduction to System
      Order Reduction Methods with Applications. In: Springer Lecture Notes in Mathematics,
      Cham: Springer International Publishing, 2014; 2114.
   7. Kononenko LI, Sobolev VA. Asymptotic expansion of slow integral manifolds. Sib.
      Math. J., 1994; 35: 1119-1132.
   8. Firstova N. Study of the critical phenomena in the model of electrochemical reactor.
      Vestnik of Samara State University, 2013; 110(9/2): 221–226. [in Russian]
   9. Shchepakina E, Korotkova O. Condition for canard explosion in a semiconductor optical
      amplifier. Journal of the Optical Society of America B: Optical Physics, 2011; 28(8):
      1988–1993. DOI: 10.1364/JOSAB.28.001988.
   10.Shchepakina E, Korotkova O. Canard explosion in chemical and optical systems. Dis-
      crete and Continuous Dynamical Systems - Series B, 2013; 18(2): 495–512. DOI:
      10.3934/dcdsb.2013.18.495.
   11.Pokrovskii A, Rachinskii D, Sobolev V, Zhezherun A. Topological degree in analysis of
      canard-type trajectories in 3-D systems. Applicable Analysis, 2011; 90(7): 1123–1139.
      DOI:10.1080/00036811.2010.511193.
   12.Sobolev V. Canard Cascades. Discrete and Continuous Dynamical Systems. Series B,
      2013; 18(3): 513–521.
   13.Babushok VI, Gol’dshtein VM, Sobolev VA. Critical conditions for thermal explosion
      with reactant consumption.Combust.Sci.andTech., 1990; 70: 81–89.
   14.Gorelov GN, Sobolev VA. Mathematical modeling of critical phenomena in thermal ex-
      plosion theory. Combust. and Flame, 1991; 87: 203– 207.
   15.Gorelov GN, Sobolev VA. Duck–trajectories in a thermal explosion problem. Appl.
      Math. Lett., 1992; 5(6): 3–6.
   16.Shchepakina E. Critical conditions of selfignition in dusty media. Journal of Advances in
      Chemical Physics, 2001; 20(7): 3-9.
   17.Shchepakina EA. Attracting/repelling integral surfaces in combustion problems. Mathe-
      matical Models and Computer Simulations, 2002; 14(3), 30–42. [in Russian]
   18.Shchepakina EA. Black swans and canards in self-ignition problem. Nonlinear Analysis:
      Real World Applications, 2003; 4: 45-50. DOI: 10.1016/S1468-1218(02)00012-3.
   19.Shchepakina E. Singular perturbations in the problem of modeling of safe combustion
      regimes. Mathematical Models and Computer Simulations, 2003; 15(8): 113–117. [in
      Russian]
   20.Shchepakina E. Canards and black swans in model of a 3-D autocatalator. Journal of
      Physics: Conference Series, 2005; 22: 194–207.
   21.Golodova ES, Shchepakina EA. Modeling of safe combustion at the maximum tempera-
      ture. Mathematical Models and Computer Simulations, 2009; 1(2): 322–334. DOI:
      10.1134/S207004820902015X.
   22.Sobolev VA, Tropkina EA. Asymptotic expansions of slow invariant manifolds and re-
      duction of chemical kinetics models. Comput. Mathematics and Math. Physics, 2012;
      52(1): 75–89.




Information Technology and Nanotechnology (ITNT-2016)                                    740
Mathematical Modelling                                  Firstova N, Shchepakina E. Study...




   23.Shchepakina EA. Critical phenomena in a model of fuel’s heating in a porous medium.
      CEUR Workshop Proceedings, 2015; 1490: 179-189. DOI: 10.18287/1613-0073-2015-
      1490-179-189.




Information Technology and Nanotechnology (ITNT-2016)                                  741