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