- Email: [email protected]
A mathematical model for current oscillations at the active-passive transition in metal electrodissolution
31
J. Ekctroanal. Chem., 347 (1993) 31-48 Elsevier Sequoia S.A., Lausanne
JEC 02439
A mathematical model for current oscillations at the active-passive transition in metal electrodissolution M.T.M. Koper and J.H. Sluyters Department of Electrochemistry, University of Utrecht, Padualaan 8, 3584 CH Utrecht (Netherlands) (Received 26 May 1992; in revised form 21 August 1992)
Abstract The classical model of Franck and FitzHugh for current oscillations at the active-passive transition in metal electrodissolution is reconsidered without making the assumption of a discontinuous change in the reaction kinetics at the Flade potential. A simple analytical calculation shows that this model cannot support sustained oscillatory behavior. A new model is presented which incorporates the role of the uncompensated ohmic cell resistance and the dissolving metal ions. The bifurcation structure of this model is studied numerically in some detail and regions of bistability, oscillations and complex (chaotic) behavior are identified.
INTRODUCTION
Recent years have witnessed a revived interest in the study of electrochemical instabilities. Electrochemical systems appear to be fruitful model objects in the search for and study of multiple steady states, oscillations and deterministic chaos. Among the wide variety of electrochemical processes capable of supporting such typical non-equilibrium behavior, the oscillations in the electrodissolution of various metals have attracted considerable attention [l-51. The electrodissolution of iron in sulphuric acid [6-91, gold in hydrochloric acid [81, copper in phosphoric acid [ll-131 or acidic chloride [141, nickel in sulphuric acid [151, cobalt in phosphoric acid [16] and zinc in aqueous hydroxides [17] have long been known to exhibit sustained oscillatory behavior. In many cases, these phenomena appear as current oscillations that take place at potentials close to the active-passive transition. Despite the large number of papers on oscillatory behavior of electrodissoluting metals, the physicochemical interpretation of these oscillations has remained poor. Most authors refer to the classical work of Franck and FitzHugh (FFI-I) [61 who developed, more than 30 years ago, a simple two-variable model which describes 0022-0728/93/$06.00
0 1993 - Elsevier Sequoia S.A. All rights reserved
32
the oscillations as a continuous cycling of the electrode between an active and a passive state. Although the model was intended primarily to explain the oscillations in the electrodissolution of iron in sulphuric acid, it is nowadays considered to be a simple general model for oscillations at the active-passive transition in a number of other electrodissolving metals as well. The key feature of the FFH model is the supposedly discontinuous change in reaction kinetics at the Flade potential, the potential which relates the active-passive transition to the local pH [3]. In a number of recent publications [18-221, several authors have carried out stability and bifurcation analyses of the FFH model and slightly modified versions, and discontinuous kinetics seems well established as a necessary prerequisite for obtaining oscillations in electrodissolving systems. However, its physical origin remains obscure, and in our opinion the necessity for discontinuous kinetics has been accepted somewhat too easily without really convincing arguments other than its ability to describe oscillations. However, Wang et al. [21] criticized the discontinuity in the FFH model and eliminated it by smoothing out its Heaviside function. This approach is an improvement in the sense that the model acquires nicer mathematical properties, but its physical meaning still remains obscure. In this paper, a reconsideration of the original FFH model is given on the basis of the above mentioned objection. Our study starts with a simple linear stability analysis of the chemical equations associated with the chemical reaction scheme assumed by FFH, a mechanism known as the dissolution-precipitation (DP) model in the passivation literature [23], by simply treating the model’s separate steps as elementary. This model allows a unique stationary state that appears to be stable under all circumstances, and does not predict sustained oscillations. The remainder of the paper will be devoted to the formulation of a new model for current oscillations at the active-passive transition in metal electrodissolution. Instead of assuming discontinuous kinetics, we incorporate into the model an uncompensated ohmic potential drop in the electrolyte solution between the working and reference electrodes. It is well known that accounting for ohmic drop is important in the interpretation of the steady-state and oscillatory behavior of iron electrodissolution in sulphuric acid [24-281, and its crucial role in destabilizing a negative faradaic impedance has been emphasized in several theoretical studies [29-311. The new model is presented in two versions. The first version simply couples the DP mechanism, as worked out by FFH, with the ohmic drop. Oscillatory behavior is then indeed obtained, but the model is dealt with only briefly as it still suffers from a number of flaws, leading to nonphysical behavior. The second version of the model should overcome these flaws by implementing the role of the dissolving metal ions and the formation of a porous salt layer on the electrode. This three-variable model gives oscillation profiles in good agreement with experiment, and for certain values of the parameters we show that the model predicts complex and chaotic oscillations. We study this model in parameter space in order to elucidate the influence of the applied potential, the bulk acidity and the mass
33
transport on the qualitative aspects of the dynamic behavior. The paper ends with a discussion which summarizes our results and conclusions and indicates some of the remaining shortcomings of our model. AI1 models presented in this paper are studied in their scaled or dimensionless form. The numerical values that we choose for the different parameters are arbitrary to the extent that they should give behavior that qualitatively mimics the experiments. We have not attempted to relate the predictions of the model quantitatively to experimental results. While the models are general and not intended to model any system in particular, they also remain low-dimensional and thus simplified descriptions of the real systems. In our opinion such a quantitative approach will in most cases be virtually out of the question because of the complexity of real systems. DISSOLUTION-PRECIPITATION
MODEL
A popular and simple model for the active-passive transition in metal electrodissolution, which also lies at the basis of the FFH model, is a reaction scheme known as the dissolution-precipitation mechanism [23]. It consists of two parallel reaction steps: k, Me-Me”++me(1) Me+nH,OG
2
[MeO,],,,
+ 2nH++ 2n e-
(2)
The faradaic current-voltage relationships of these reactions are depicted graphically in Fig. 1. In the following, these steps will be treated as elementary.
‘F
I
Me-Me
m++me-
MeOn+2nH++2ne.
Fig. 1. Typical current-voltage relationships for reactions (1) and (2) of the DP mechanism.
The equilibrium potential of reaction (2), which depends linearly on the pH, is known as the Flade potential [3]. For reasons of mathematical simplicity, we assume m = 1 and 2n = 1 in the remainder of the paper. This may be quantitatively incorrect, but this does not matter in our qualitative approach. It also means that we assume reactions (2) to be first order in [H+l. It is shown below that this has no consequences for the number of steady states nor for their stability properties. It is also assumed that the negative counter-ion present in the electrolyte solution is the monovalent species A-. Franck and FitzHugh [6] make a number of simplifications in the description of the mass transport. (a) The concentration of Me*+ is so small that its contribution to the migration current is negligible. (b) H+ carries all the current in the solution, which is equivalent to assuming that its transference number is unity at every location in the electrolyte solution. Given assumption (a) this is not as inaccurate as might seem at first sight. If electroneutrality prevails at every location x in the solution (i.e. the electrical double layer has negligible thickness), then [H’lx = [A-], so that the transference
number for H+ becomes
b+[H+lx
tH+(X) =h,+[H+],
A n+
= &++A,-
+ A,-[A-],
=
1
if An+>> A,-, which is not unreasonable because H+ is usually much more mobile than A-. (c) Deviations from a diffusion layer with thickness 6 and with a linear concentration profiles can be neglected. Assumptions (a)-(c) allow us to derive a simple equation for the time evolution of the surface concentration c,(t) of an arbitrary species j in the following way. From Fig. 2 it follows that the shaded area ;[c,(t
+dt)
-c&)]6
(3)
should equal
(Jctif - Jmigr + JF,j)dt =
Dkb”,k- cowl s
(4)
where the J are particle fluxes, J,j is the Faraday flux associated with particle j, Z is total current, tj is the transference number, A is the total electrode surface area and F is Faraday’s constant. Equating (3) and (4) and taking the limit dt -P 0, we obtain 2tjz 2 dc, 20 -= (5) F(C~UL co) FSA + gJ,,j dt It should be noted that this equation is strictly valid only if the transference number rj #f(x). This is true for a binary electrolyte under electroneutrality, a
electrode
surface
-
Jdif
I Cbulk
I I I I
61
-x
Fig. 2. The change in the linear concentration profile near the electrode in the time interval dt and the various fluxes in its determination.
property contained in assumptions (a) and (b). Further, it should be noted that eqn. (5) disregards the (possible) existence of a diffusion potential across the diffusion layer. The total faradaic current IF = FA(.Z,,C,,+ .ZF,C2j) associated with reactions (1) and (2) will be ZF=AF[kl(l-t’)+k2(1-8)-k_,e[Z-Z+]]
(6)
where 8 is the fractional coverage of the electrode with the passive oxide layer MeO,. Since FFH assume that the electrode potential is truly constant, even under oscillatory conditions, there is no charging current and I, equals the total current I: I, = z
(7)
From eqns. Cl), (2) and W(7) and assumptions (a)-(c), we can write the following two evolution equations for the H+ concentration [H+], at the electrode-electrolyte interface and the passive layer coverage 8: d[H+]a 2D - dt = +H+lmc pg
= k,(l
- [H+
Id - ;&,(I) + JF,(2)) + g,,(2)
- 0) - k_#[H+],
where D is the diffusion coefficient of H+ and ~3 is the surface concentration of MeO, at maximum coverage (recall that tH+= 1). The passive layer is thus assumed to be a two-dimensional lattice gas of negligible thickness following
36
langmuirian adsorption laws, i.e. without particle interactions. eqn. (8), we obtain d[H+],,
Inserting eqn. (6) in
20
- [H+],) - &I6’) (10) dt At this stage we note a mistake in FFH’s original derivation. They assumed a priori that reaction (2) does not contribute to the material balance for [H+l because this contribution would be small (the last term in eqn. (8) is supposed to be zero), although they did account for the current associated with this reaction being carried by the H+ by means of migration. However, the two terms cancel exactly in the [H+], balance because of assumption (b) in the model, as is clear in the above derivation. Therefore the balance in the original FFH model is incorrect. If we denote [H+lbulk by h, and [H+] by h, and we rescale time according to 2Dt/A --) t, we obtain -
=
@H+ltm
dh dt = d(h, -h)
- d”*(l-
(1)
8)k,
de dt = E[ k,( 1 - e) - k_*eh] where d = A/S’, k; A”*/D
E = A’/*/2/3
and the rate constants are resealed according to
+ ki.
We shall refer to this model as model I. Model I also lies at the basis of the FFH model. At this stage, FFH assume that the passivation kinetics can be treated discontinuously and they obtain two sets of two coupled differential equations each, one valid for potentials lower than the Flade potential and one valid for higher potentials. In this way limit-cycle oscillations were obtained which resemble those of the experiments rather well. It should be noted that the introduction by FFH of an overpotential, which is the difference between the electrode potential and the (pH-dependent) Flade potential, is nothing more than a reformulation of the problem and is in fact completely redundant. We prefer to work with the equation for [H+], because it is simpler. Detailed discussions of the FFH model can be found in the literature [4,5,18-221. The stability analysis of model I can be carried out without the discontinuity assumption. The stationary states h, and OS, of model I can be calculated explicitly: 1
%s = 1+
= -B+
h ss
(k_,/k,)h,
(11)
[B*+dhh(k-2/k2)11’* 2( k-,/k,)
where B = 1 - hb(k_*/k2) + k1(k_2/k2)/d’/2 but it is chemically irrelevant.)
(12) (There is another solution for h,
37
200
x
0’
”
(b)
(a)
Fig. 3. Two typical current-voltage characteristics for model I with ey = 0, ey = 20, d = 1, h, = 1, /~!~/k! = 0.1, I = q(l- 0&,, q = 100 and c any value: (a) /cp = 0.001; (b) kf = 0.01.
The stability of these stationary states can be investigated using linear stability analysis [5]. By superimposing small perturbations such that h = h,, + 6h and 8 = 0, + 619,we obtain two linear homogeneous differential equations -d
d”2k
-Ek,e,
1
-e( k, + k_,h,)
(13)
whose Jacobian matrix has the following eigenvalues 5: t2+tT+A=0
(14)
where T = d + c( k, + k_,h,,) A = de( k, + k_,h,)
+ cd1’2kIk2es
Because both T and A are always positive, the eigenvalues are strictly negative and consequently all perturbations die out; the stationary state is always asymptotically stable. The decay is oscillatory when 44 > T2, but no sustained oscillations are possible in model I. This result is in contrast with the results from the original FFH model and shows that the oscillations predicted by that model must be caused by the assumption that the passivation kinetics are discontinuous. Figure 3 shows two typical current-voltage characteristics which qualitatively resemble typical passivation characteristics. In Fig. 3(b) we see a mass-transferlimited dissolution due to the limited rate at which the H+ ions diffuse towards the electrode-electrolyte interface (the limiting current is proportional to dh,). For the rate constants, we have assumed Tafel laws according to k, = k; exp[ (e - e;)/2] k, = k; exp[ (e - e”,)/2] k_,=k”_,
exp[-(e-e;)/21
(1%
38
where e = FE/RT and E is the electrode potential. In Fig. 3 we took e’, = 0.0 for convenience. It should be noted that a reaction order of [H+] different from unity, say p, does not alter the number of steady states nor their stability properties. The expression for h, would simply be eqn. (112 to the power l/p. This is a specific property of model I and a direct result of assumptions (a)-(c), and it will not necessarily hold for the more elaborate models we shall consider below. DISSOLUTION-PRECIPITATION
MODEL WITH OHMIC DROP
Electrodissolution is usually accompanied by high electric currents and thus by appreciable ohmic drops in the electrolyte solution. It has been argued by several authors [16,24-281 that ohmic drop plays an important role in the interpretation of experimentally obtained current-voltage relationships. In many cases, a hysteresis which can be ascribed entirely to the potential drop between the working and reference electrode is observed [27]. Moreover, it is known that a negative faradaic impedance (dZ,/d E)-‘, in terms of which passivation can be understood (see Fig. 3), needs an ohmic series resistance to give rise effectively to instabilities such as hysteresis (or bistability, which is the same) or spontaneous current oscillations [24,29-311. Let R, be the ohmic resistance between working and reference electrode. If the potential on the electrode is identified with zero, the potential just outside the electrical double layer with E and the (controlled) potential of the reference electrode with V, eqn. (7) has to be modified according to I=-
V-E R,
=I,+z,=ACg
+z,
Using the dimensionless equations are obtained: de -= dt dh
potential
+q[k,(l-e)+k,(l+k_,Bh] 1 u-e
x=~(h,-h)-d1~2--(-)+d1’2[k2(1-8)-k_28h]
de
e = FE/RT, the following scaled differential
dt = +,(1-
(11)
e) - k,ehl
with LJ= FV/RT, r = there is no additional double-layer capacity referred to as “model quasi-steady state for
2DC,R,
and q = F2A’12/2 RTC,. The model assumes that ohmic drop associated with the passive layer and that the per unit surface area C, is a constant. This model will be II”. A useful version of model II is that corresponding to a 8 (i.e. dfZ/dt = O), which can then be equated to eqn. (11).
39
time
0
10
0
time
I”
(b)
(0)
Fig. 4. Typical current oscillation profiles from model II, with q = 100, ky = 0.01, ki = 0.01, k!!, = 0.001, r = 0.1, u = 30, d = 1, h, = 1, ef = 0, et = 20 and u = 30: (a) E = 200, (b) E = 01.
Numerical integration shows that it is indeed possible to obtain oscillatory solutions for model II (a Runge-Kutta routine with automatic step-size adaptor was used to solve the equations [32]). Figure 4 shows two examples, one in the quasi-steady-state assumption for 8. Although these waveforms do not resemble the experimental ones very closely, a more serious objection to model II is that h can become negative during these oscillations. This is also manifested in the oscillation profile in Fig. 4(b), in which a high current plateau exists at a higher current than the diffusion-limited current plateau in the steady-state curves of Fig. 3. This is caused by a negative H+ concentration. This artifact is a result of the unrealistic assumptions (a) and (b) discussed in the foregoing section. Clearly, a more realistic description becomes desirable. Therefore we end our study of model II at this stage and turn to a more complete description of the physical processes occurring during metal electrodissolution in the next section. IMPROVED
DISSOLUTION-PRECIPITATION
MODEL WITH OHMIC DROP
In this section models I and II are improved by replacing assumptions (a> and (b) with more realistic mass-transport descriptions. In this model all three ionic species (H+, Me+ and A-) migrate and carry their part of the current flowing through the cell, each in proportion to their transference number. Electroneutrality still prevails at every location x in the electrolyte solution:
WL + [Me% = [A-lx
(17)
which eliminates [A-] as a variable. The local transference j is given by
ti(x)= [H+],(h,++
[jldj AA-) + [Me’]J
hr,,le++ AA-)
number tj(x) of species
(18)
40
When Me+ is introduced as a variable in the model description, under conditions where the electrode is active, the Me+ concentration at the electrode-electrolyte interface may reach extremely high values. In reality, however, the solubility product K,, of Me+A- is exceeded, and a porous salt film precipitates onto the electrode. This has been observed, for instance, in the electrodissolution of iron in acidic sulphate electrolytes [33,34]. Then the rate at which the salt layer can dissolve into the neighboring solution determines the rate of the metal electrodissolution. Therefore formation of a salt layer may cause a limiting current in the current-voltage characteristic because the salt film dissolution is essentially independent of potential [35]. The exact kinetic equations for the formation of a salt film are complex [36]; presumably these kinetics are not essential for understanding oscillations because the salt film has no passivating properties. It is usually assumed that the passivating oxide layer forms underneath the salt film, closest to the electrode material. In order to account for the limited metal electrodissolution rate due to salt film formation, but at the same time to avoid complicated equations which may distract from the essence of the unstable behavior, we take the following opportunistic approach. We assume that a porous salt film has precipitated onto the electrode. To a first very crude approximation we assume that its conducting properties are no different from those of the electrolyte solution. This means there is no additional ohmic drop due to the salt layer. A metal ion that has come into being through reaction (1) travels through the salt film and reaches the salt film-electrolyte interface. Its probability of dissolving into the electrolyte is unity if the metal ion concentration in the neighboring solution is zero and zero when this concentration is beyond some critical value [Me+],,,. For intermediate situations, we assume a linear relationship to hold. This can be translated into a corrected rate constant k::
(19) where 0 I [Me+], I [Me+],,,. The critical concentration solubility product K,, of the MeA salt:
[Me+],,,
is related to the
[Me+lcrdA-l~= K,, Using eqn. (17), this is equivalent to
[Me+],,i, = -i[H’]o
+ (i[H+]i
+K,I)1’2
(21)
At this stage let us note that eqn. (5) can no longer be used to derive differential equations for [H+], and [Me+], since in this new model both t,+ and tMe+ are functions of the distance x from the electrode. There are two alternatives for deriving the new differential equations: (i) one can introduce an averaged transference number fj = (l/b)/ttj(x)dx and insert this in eqn. (5), or (ii) instead of assuming a linear concentration profile one can assume the concentration to be
41
uniform in a region of size 6 adjacent to the electrode and to attain its bulk value outside this region. The latter case would simply give eqn. (5) with the right-hand side divided by a factor of 2. In scaled terms, both approaches yield the same result except that the expressions for the various parameters differ by a factor of 2. Therefore, denoting [Me+], by u, and assuming for simplicity the quasi-steadystate assumption for 8, the following three scaled differential equations can be obtained: de v-e - = -qk,*(l -e> dt r dh dt = d( h, - h) - thd1’2 du dt-
- --LYdu + d”‘( 1 - 0)k:
(III) - t,d1’2
where a=_=-=-
D Me+
A Me+
A,-
D H+
A,+
A,+
h t, =
h( 1 +(Y) + 2~~4 (YU
t”=
h(l+cr)
+2au
u, = - $h + ( $h2 + K,$” and 8 follows from eqn. (11) with h for h,,. It should be noted that we disregard the salt film thickness in order to keep the model mathematically as simple as possible. In Figs. 5-10 we summarize some of the steady-state and dynamic properties of model III. In Fig. 5 the steady-state current-voltage curve is given for low ohmic resistance. The curve is stable and smooth without abrupt transitions or hysteresis phenomena. We have also simulated some linear sweep voltammograms for higher ohmic resistance. In the simulation, we took v = v” + ct. Results are depicted in Fig. 6. For fast sweep rates c, there is an abrupt transition at the active-passive transition; at lower sweep rates, this transition is accompanied by oscillations. The simulated sweep is in qualitative agreement with several experimental sweeps [3,16,17,231.
42
Fig. 5. Steady-state current-voltage relationship for model III at low ohmic resistance K,, = 16 and a = 0.2. Other parameter values are as in Fig. 4.
r = 0.02,
Oscillatory time series associated with a typical limit cycle solution of model III are shown in Fig. 7. The relaxation oscillation for the Z-f curve resembles that of the experiments. It should be noted from the 0-t curve that the electrode indeed effectively switches between an entirely active and an entirely passive state. For certain adjustments of the various parameters, the model predicts complex or even chaotic oscillations. Figure 8(a) shows a typical “mixed-mode” oscillation, which is made up of small and large excursions. As a function of one parameter, mixed-mode oscillations are known to follow complicated sequences but we shall not dwell upon them here (but see ref. 5). A chaotic solution of model III is illustrated in Fig. 8(b) with its Z-t behavior.
300
3oc
I
0 1"
5.
(a)
"
__
Yb
10
"
35
(b)
Fig. 6. Linear sweep voltammograms for r = 0.1: (a) sweep rate c = 1; (b) c = 0.04. Other parameter values are as in Fig. 5.
43
h
0
tme
0
50
time
5r
,L 0
50
(b
(a)
I
0
ti me
50
time
II
(c)
Fig. 7. Typical oscillatory time series for limit cycle solution of model III with u = 30, r = 0.1, h, = 1.2 and other parameter values as in Fig. 5; (a) Z-r behavior; (b) h-t behavior; (c) u-t behavior; (d) 13-t behavior.
Figures 9 and 10 illustrate the bifurcation behavior of model III in more detail. In Fig. 9 we show two typical parameter diagrams which map the stability of the steady-state solution(s) of model III onto the d-v and h,-u planes respectively. Within the curves labelled SN (saddle-node bifurcations), model III supports three simultaneous steady-state solutions. One of hem is always unstable. At the Hopf t
200
I
I/ time (b)
Fig. 8. Complex oscillations in model III; (a) mixed-mode oscillation for u = 30.52, r = 0.1 and d = 1; (b) chaotic oscillation for u = 24, r = 0.032505 and d = 0.5. Other parameter values are as in Fig. 5.
44
12
2 Hopf
0 “b
4
I
1
4
04
25
20
30
A”
35
0 20
25
_
v
30
35
(b)
(a) _
Fig. 9. Linear stability diagrams of model III which indicate regions of unique stable steady-state, oscillations and multiple steady states, and associated bifurcations (r = 0.1): (a) d-u diagram for h, = 1; (b) hb-u diagram for d = 1. Other parameter values are as in Fig. 5.
0
0.5
-
‘.o *
1.5
2.0
Ic)
Fig. 10. Bifurcation diagrams for model III: full curves, stable steady states or minimum and maximum of stable oscillation; broken curves, unstable steady states; Hopf, Hopf bifurcation; SN, saddle-node bifurcation; shaded region MMO, region where mixed-mode oscillations are obtained. (a) Z-h, diagram, u = 30, d = 1; (bl Z-u diagram, hi, = 1, d = 1; (cl Z-d diagram, u = 30, h, = 1.
45
bifurcation, a stable steady state becomes unstable and the system starts oscillating spontaneously [5]. Three representative sections of Fig. 9 are shown in Figs. 10(a)-10(c), where we have also indicated the occurrence of mixed-mode oscillations. These figures show the influence of the various parameters on the oscillation amplitude and also whether the transition from steady-state to oscillatory behavior is abrupt or smooth. For instance, with increasing bulk acidity, h,, the amplitude and the period of oscillation decreases, in agreement with the experimental results for iron in H,SO, [7]. (In short, the bifurcation diagrams of Fig. 10 were computed as follows. First, numerical integration provided an estimation of a steady-state expected to have a large basin of attraction. This result was used as an initial estimate for the Newton-Raphson method of iteratively solving the steady-state equations [32], and as a starting point for the one-parameter continuation of the steady-state. During the one-parameter continuation, the linear eigenvalues of the steady state were monitored and checked for passage through a saddle-node bifurcation (one zero eigenvalue) or a Hopf bifurcation (complex pair of eigenvalues with zero real part). Oscillation extrema were obtained from numerical integration. Figure 9 was drawn from an assembly of such one-parameter continuations. Usually, the electrode potential e was used as a parameter since this avoids instabilities in the NewtonRaphson algorithm close to saddle-node bifurcations.) DISCUSSION AND CONCLUSIONS
In this paper we have described a simple and general mathematical model to explain the current oscillations that are often observed at the active-passive transition in various metal electrodissolutions. The oscillations in the model arise from the coupling of a negative faradaic impedance due to the passivation kinetics with an uncompensated ohmic series resistance attributable to the electrolyte solution. It has also been shown that the oscillations in the.classical FFH model must be due to their discontinuous treatment of the passivation kinetics. When treated without assuming this discontinuity, this model can describe only stable steady-state behavior and never oscillations. Some authors have stated that discontinuous kinetics are a necessary prerequisite for obtaining oscillations. Apart from the objection that it is unclear how such discontinuous kinetics would arise physically, we believe that the present model shows that it is not necessary to incorporate discontinuous kinetics in order to obtain oscillations. Once again we emphasize that our model suffers from a number of oversimplifications and therefore it should not be judged on its quantitative properties. Although the experimental behavior is usually of low dimension (limit cycles or low-dimensional complex behavior), this only means that a low-dimensional qualitative model will suffice; a quantitative approach needs a far more detailed and laborious description of all the physicochemical processes in the system. In particular, the DP mechanism is surely not a reliable candidate for a quantitative
46
treatment of electrodissolution and passivation kinetics, but only provides the simplest possible interpretation of the negative faradaic impedance associated with the passivation phenomenon. For instance, the present model cannot explain the bistability and hysteresis that Epelboin et al. [37] observed in the current-voltage relationship for iron in H,SO, at low pH after correction for the ohmic drop. They were able to reproduce this observation qualitatively by assuming a more complex dissolution scheme which involves the formation of a non-passivating and a passivating intermediate. Such hysteresis can also be reproduced when an attracting mean field interaction between the adsorbing particles is introduced [38], i.e. a kind of Frumkin adsorption. When coupled to a second differential equation, Frumkin adsorption can give rise to limit cycle solutions (oscillations) without the explicit introduction of the electrode potential or mass transport (see, for example, the work by Talbot and Oriani [39], Brynn-Hilbert and Murphy [40] and Kado and Kunitomi [41]). Such an approach seems to us to be of a highly speculative nature since it often lacks experimental confirmation and it disregards many non-trivial processes appropriate to electrochemical systems. It should be noted that our model also differs in some important features from a recent model proposed by Haim et al. [42] who described periodic phenomena during anodic nickel dissolution. Their model applies to the transpassive region instead of the active-passive transition, it does not account for mass transport, which plays an important role in our model, and finally it describes potential oscillations under galvanostatic conditions whereas our model describes current oscillations under fixed potential conditions. Another point is that our treatment of the passive layer as a two-dimensional lattice gas that does not need a spatial description is certainly not always justified. The formation and dissolution of the passive layer is likely to give rise to islands of activation/passivation resulting in non-uniform current and concentration distributions. This would introduce a spatial component in the mathematical description of the system. Also, our simplification of the diffusion problem in terms of a linear diffusion layer misses an essential point. Under time-dependent conditions, such as spontaneous oscillations, the linear diffusion-layer approximation is no longer valid. There will be a kind of delay or relaxation associated with the diffusion process in trying to restore the linear or steady-state concentration profile, and a simple way in which this can be accounted for qualitatively is through the introduction of a second diffusion layer [43]. We observed that, with the introduction of this diffusion-layer relaxation, mixed-mode oscillations can appear where simple oscillations were obtained in the original model. This means that this effect can be very important in the interpretation of complex and chaotic oscillations in electrochemical systems. In fact, it has been shown for the indium thiocyanate electrochemical oscillator that, with the incorporation of diffusion-layer relaxation, excellent agreement with the experimentally observed complex bifurcation sequences is obtained [44]. Finally, throughout this paper we have assumed that the uncompensated ohmic cell resistance is entirely attributable to the electrolyte solution and could thus be
41
considered a constant. However, both the passive layer and the salt layer will contribute to the total ohmic drop in a non-negligible fashion, and it is not very likely that this contribution will remain unchanged during the oscillations. It also remains to be seen in what respect the complex and chaotic oscillations predicted by the model relate to the experimental chaos. Although they show some unusual similarities with the experimental results for a number of systems, this property of the model should not be overestimated. Many metal electrodissolutions display bifurcation scenarios that may still be traced to unique properties of these systems themselves. Nevertheless, despite these objections and shortcomings, the present model may prove to be a helpful tool in the interpretation and understanding of instabilities and (a)periodic oscillations at the active-passive transition in many metal electrodissolutions. REFERENCES 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
26 27 28 29
ES. Hedges and J.E. Myers, The Problem of Physico-chemical Periodic&y, Arnold, London, 1926. U.F. Franck, Z. Elektrochem., 62 (1958) 649. K-J. Vetter, Elektrochemische Kinetik, Springer, Berlin, 1961. J. Wojtowicz in J.O’M. Bockris and B.E. Conway (Eds.), Modern Aspects of Electrochemistry, Vol. 9, Plenum Press, New York, 1973, p. 47. SK. Scott, Chemical Chaos, Clarendon Press, Gxford, 1991. U.F. Franck and R. FitsHugh, Z. Elektrochem., 65 (1961) 156. J.J. Podesta, R.C.V. Piatti and A.J. Arvia, J. Electrochem. Sot., 126 (1979) 1363. G. Beni and S. Hackwood, J. Appl. Electrochem., 14 (1984) 623. C.B. Diem and J.L. Hudson, AIChE J., 33 (19871218. J.J. Podesta, R.C.V. Piatti and A.J. At-via, Electrcchim. Acta, 24 (1979) 633. F.N. Albahadily and M. Schell, J. Chem. Phys., 88 (1988) 4312. F.N. Albahadily, J. Ringland and M. Schell, J. Chem. Phys., 90 (1989) 813. M. Schell and F.N. Albahadily, J. Chem. Phys., 90 (1989) 822. M.R. Bassett and J.L. Hudson, J. Phys. Chem., 93 (1989) 2731. 0. Lev, A. WoIffberg, M. Sheintuch and L.M. Pismen, Chem. Eng. Sci., 43 (1988) 1339. D. Sazou and M. Pagitsas, J. Electroanal. Chem., 323 (19921 247. M.C.H. McKubre and D.D. Macdonald, J. Electrochem. Sot., 128 (1981) 524. J.B. Talbot, R.A. Oriani and M.J. Dicarlo, J. Electrochem. Sot., 132 (1985) 1545. SM. Sharland, Corros. Sci., 27 (1987) 289. A.J. Pearlstein and J.A. Johnson, J. Electrochem. Sot., 136 (1989) 1290. Y. Wang, J.L. Hudson and N.I. Jaeger, J. Electrochem. Sot., 137 (1990) 485. M. Pagitsas and D. Sazou, Electrochim. Acta, 36 (1991) 1301. N. Sato and G. Okamoto in JG’M. Bockris, B.E. Conway, E. Yeager and R.E. White (Eds.), Comprehensive Treatise of Electrochemistry, Vol. 4, Plenum Press, New York, 1984, p. 193. I. Epelboin, C. Gabrielli, M. Keddam. J.-C. Lestrade and H. Takenouti, J. Electrochem. Sot., 119 (19721 1632. I. Epelboin, C. Gabrielli and M. Keddam in J.O’M. Bockris, B.E. Conway, E. Yeager and R.E. White (Eds.), Comprehensive Treatise of Electrochemistry, Vol. 9, Plenum Press, New York, 1984, p. 151. C. Gabrielli, M. Keddam, C. Rakotomavo and H. Takenouti, J. Electrochem. Sot., 127 (1980) 2648. C.G. Law Jr. and J. Newman, J. Electrochem. Sot., 126 (1979) 2150. P. Russell and J. Newman, J. Electrochem. Sot., 134 (1987) 1051. H. Degn, Trans. Faraday Sot., 64 (1968) 1348.
48 30 R. de Levie, J. Electroanal. Chem., 25 (1970) 257. 31 M.T.M. Koper, Electrochim. Acta, 37 (1992) 1771. 32 W.H. Press, B.P. PIannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes, The Art of Scientific Computing, Cambridge University Press, Cambridge, 1989. 33 R. A&ire, D. Emsberger and T.R. Beck, J. Electrochem. Sot., 125 (1978) 1383. 34 M.E. Orazem and M.G. Miller, J. Electrochem. Sot., 134 (1987) 393. 35 L. Kiss, Kinetics of Electrochemical Metal Dissolution, Elsevier, Amsterdam, 1988. 36 P. Russell and J. Newman, J. Electrochem. Sot., 133 (1986) 59. 37 I. Epelboin, C. Gabrielli, M. Keddam and H. Takenouti, Electrochim. Acta, 20 (1975) 913. 38 G.L. Griffen, J. Electrochem. Sot., 131 (1984) 18. 39 J.B. Talbot and R.A. Oriani, Electrochim. Acta, 30 (1985) 1277. 40 D. Btynn Hilbert and S.V. Murphy, J. Electtochem. Sot., 138 (1991) L30. 41 T. Kado and N. Kunitomi, J. Electrochem. Sot., 138 (1991) 3312. 42 D. Haim, 0. Lev, L.M. Pismen and M. Sheintuch, J. Phys. Chem., 96 (1992) 2676. 43 M.T.M. Koper and P. Gaspard, J. Chem. Phys., 96 (1992) 7797. 44 M.T.M. Koper, P. Gaspard and J.H. Sluyters, J. Chem. Phys., in press.
J. Ekctroanal. Chem., 347 (1993) 31-48 Elsevier Sequoia S.A., Lausanne
JEC 02439
A mathematical model for current oscillations at the active-passive transition in metal electrodissolution M.T.M. Koper and J.H. Sluyters Department of Electrochemistry, University of Utrecht, Padualaan 8, 3584 CH Utrecht (Netherlands) (Received 26 May 1992; in revised form 21 August 1992)
Abstract The classical model of Franck and FitzHugh for current oscillations at the active-passive transition in metal electrodissolution is reconsidered without making the assumption of a discontinuous change in the reaction kinetics at the Flade potential. A simple analytical calculation shows that this model cannot support sustained oscillatory behavior. A new model is presented which incorporates the role of the uncompensated ohmic cell resistance and the dissolving metal ions. The bifurcation structure of this model is studied numerically in some detail and regions of bistability, oscillations and complex (chaotic) behavior are identified.
INTRODUCTION
Recent years have witnessed a revived interest in the study of electrochemical instabilities. Electrochemical systems appear to be fruitful model objects in the search for and study of multiple steady states, oscillations and deterministic chaos. Among the wide variety of electrochemical processes capable of supporting such typical non-equilibrium behavior, the oscillations in the electrodissolution of various metals have attracted considerable attention [l-51. The electrodissolution of iron in sulphuric acid [6-91, gold in hydrochloric acid [81, copper in phosphoric acid [ll-131 or acidic chloride [141, nickel in sulphuric acid [151, cobalt in phosphoric acid [16] and zinc in aqueous hydroxides [17] have long been known to exhibit sustained oscillatory behavior. In many cases, these phenomena appear as current oscillations that take place at potentials close to the active-passive transition. Despite the large number of papers on oscillatory behavior of electrodissoluting metals, the physicochemical interpretation of these oscillations has remained poor. Most authors refer to the classical work of Franck and FitzHugh (FFI-I) [61 who developed, more than 30 years ago, a simple two-variable model which describes 0022-0728/93/$06.00
0 1993 - Elsevier Sequoia S.A. All rights reserved
32
the oscillations as a continuous cycling of the electrode between an active and a passive state. Although the model was intended primarily to explain the oscillations in the electrodissolution of iron in sulphuric acid, it is nowadays considered to be a simple general model for oscillations at the active-passive transition in a number of other electrodissolving metals as well. The key feature of the FFH model is the supposedly discontinuous change in reaction kinetics at the Flade potential, the potential which relates the active-passive transition to the local pH [3]. In a number of recent publications [18-221, several authors have carried out stability and bifurcation analyses of the FFH model and slightly modified versions, and discontinuous kinetics seems well established as a necessary prerequisite for obtaining oscillations in electrodissolving systems. However, its physical origin remains obscure, and in our opinion the necessity for discontinuous kinetics has been accepted somewhat too easily without really convincing arguments other than its ability to describe oscillations. However, Wang et al. [21] criticized the discontinuity in the FFH model and eliminated it by smoothing out its Heaviside function. This approach is an improvement in the sense that the model acquires nicer mathematical properties, but its physical meaning still remains obscure. In this paper, a reconsideration of the original FFH model is given on the basis of the above mentioned objection. Our study starts with a simple linear stability analysis of the chemical equations associated with the chemical reaction scheme assumed by FFH, a mechanism known as the dissolution-precipitation (DP) model in the passivation literature [23], by simply treating the model’s separate steps as elementary. This model allows a unique stationary state that appears to be stable under all circumstances, and does not predict sustained oscillations. The remainder of the paper will be devoted to the formulation of a new model for current oscillations at the active-passive transition in metal electrodissolution. Instead of assuming discontinuous kinetics, we incorporate into the model an uncompensated ohmic potential drop in the electrolyte solution between the working and reference electrodes. It is well known that accounting for ohmic drop is important in the interpretation of the steady-state and oscillatory behavior of iron electrodissolution in sulphuric acid [24-281, and its crucial role in destabilizing a negative faradaic impedance has been emphasized in several theoretical studies [29-311. The new model is presented in two versions. The first version simply couples the DP mechanism, as worked out by FFH, with the ohmic drop. Oscillatory behavior is then indeed obtained, but the model is dealt with only briefly as it still suffers from a number of flaws, leading to nonphysical behavior. The second version of the model should overcome these flaws by implementing the role of the dissolving metal ions and the formation of a porous salt layer on the electrode. This three-variable model gives oscillation profiles in good agreement with experiment, and for certain values of the parameters we show that the model predicts complex and chaotic oscillations. We study this model in parameter space in order to elucidate the influence of the applied potential, the bulk acidity and the mass
33
transport on the qualitative aspects of the dynamic behavior. The paper ends with a discussion which summarizes our results and conclusions and indicates some of the remaining shortcomings of our model. AI1 models presented in this paper are studied in their scaled or dimensionless form. The numerical values that we choose for the different parameters are arbitrary to the extent that they should give behavior that qualitatively mimics the experiments. We have not attempted to relate the predictions of the model quantitatively to experimental results. While the models are general and not intended to model any system in particular, they also remain low-dimensional and thus simplified descriptions of the real systems. In our opinion such a quantitative approach will in most cases be virtually out of the question because of the complexity of real systems. DISSOLUTION-PRECIPITATION
MODEL
A popular and simple model for the active-passive transition in metal electrodissolution, which also lies at the basis of the FFH model, is a reaction scheme known as the dissolution-precipitation mechanism [23]. It consists of two parallel reaction steps: k, Me-Me”++me(1) Me+nH,OG
2
[MeO,],,,
+ 2nH++ 2n e-
(2)
The faradaic current-voltage relationships of these reactions are depicted graphically in Fig. 1. In the following, these steps will be treated as elementary.
‘F
I
Me-Me
m++me-
MeOn+2nH++2ne.
Fig. 1. Typical current-voltage relationships for reactions (1) and (2) of the DP mechanism.
The equilibrium potential of reaction (2), which depends linearly on the pH, is known as the Flade potential [3]. For reasons of mathematical simplicity, we assume m = 1 and 2n = 1 in the remainder of the paper. This may be quantitatively incorrect, but this does not matter in our qualitative approach. It also means that we assume reactions (2) to be first order in [H+l. It is shown below that this has no consequences for the number of steady states nor for their stability properties. It is also assumed that the negative counter-ion present in the electrolyte solution is the monovalent species A-. Franck and FitzHugh [6] make a number of simplifications in the description of the mass transport. (a) The concentration of Me*+ is so small that its contribution to the migration current is negligible. (b) H+ carries all the current in the solution, which is equivalent to assuming that its transference number is unity at every location in the electrolyte solution. Given assumption (a) this is not as inaccurate as might seem at first sight. If electroneutrality prevails at every location x in the solution (i.e. the electrical double layer has negligible thickness), then [H’lx = [A-], so that the transference
number for H+ becomes
b+[H+lx
tH+(X) =h,+[H+],
A n+
= &++A,-
+ A,-[A-],
=
1
if An+>> A,-, which is not unreasonable because H+ is usually much more mobile than A-. (c) Deviations from a diffusion layer with thickness 6 and with a linear concentration profiles can be neglected. Assumptions (a)-(c) allow us to derive a simple equation for the time evolution of the surface concentration c,(t) of an arbitrary species j in the following way. From Fig. 2 it follows that the shaded area ;[c,(t
+dt)
-c&)]6
(3)
should equal
(Jctif - Jmigr + JF,j)dt =
Dkb”,k- cowl s
(4)
where the J are particle fluxes, J,j is the Faraday flux associated with particle j, Z is total current, tj is the transference number, A is the total electrode surface area and F is Faraday’s constant. Equating (3) and (4) and taking the limit dt -P 0, we obtain 2tjz 2 dc, 20 -= (5) F(C~UL co) FSA + gJ,,j dt It should be noted that this equation is strictly valid only if the transference number rj #f(x). This is true for a binary electrolyte under electroneutrality, a
electrode
surface
-
Jdif
I Cbulk
I I I I
61
-x
Fig. 2. The change in the linear concentration profile near the electrode in the time interval dt and the various fluxes in its determination.
property contained in assumptions (a) and (b). Further, it should be noted that eqn. (5) disregards the (possible) existence of a diffusion potential across the diffusion layer. The total faradaic current IF = FA(.Z,,C,,+ .ZF,C2j) associated with reactions (1) and (2) will be ZF=AF[kl(l-t’)+k2(1-8)-k_,e[Z-Z+]]
(6)
where 8 is the fractional coverage of the electrode with the passive oxide layer MeO,. Since FFH assume that the electrode potential is truly constant, even under oscillatory conditions, there is no charging current and I, equals the total current I: I, = z
(7)
From eqns. Cl), (2) and W(7) and assumptions (a)-(c), we can write the following two evolution equations for the H+ concentration [H+], at the electrode-electrolyte interface and the passive layer coverage 8: d[H+]a 2D - dt = +H+lmc pg
= k,(l
- [H+
Id - ;&,(I) + JF,(2)) + g,,(2)
- 0) - k_#[H+],
where D is the diffusion coefficient of H+ and ~3 is the surface concentration of MeO, at maximum coverage (recall that tH+= 1). The passive layer is thus assumed to be a two-dimensional lattice gas of negligible thickness following
36
langmuirian adsorption laws, i.e. without particle interactions. eqn. (8), we obtain d[H+],,
Inserting eqn. (6) in
20
- [H+],) - &I6’) (10) dt At this stage we note a mistake in FFH’s original derivation. They assumed a priori that reaction (2) does not contribute to the material balance for [H+l because this contribution would be small (the last term in eqn. (8) is supposed to be zero), although they did account for the current associated with this reaction being carried by the H+ by means of migration. However, the two terms cancel exactly in the [H+], balance because of assumption (b) in the model, as is clear in the above derivation. Therefore the balance in the original FFH model is incorrect. If we denote [H+lbulk by h, and [H+] by h, and we rescale time according to 2Dt/A --) t, we obtain -
=
@H+ltm
dh dt = d(h, -h)
- d”*(l-
(1)
8)k,
de dt = E[ k,( 1 - e) - k_*eh] where d = A/S’, k; A”*/D
E = A’/*/2/3
and the rate constants are resealed according to
+ ki.
We shall refer to this model as model I. Model I also lies at the basis of the FFH model. At this stage, FFH assume that the passivation kinetics can be treated discontinuously and they obtain two sets of two coupled differential equations each, one valid for potentials lower than the Flade potential and one valid for higher potentials. In this way limit-cycle oscillations were obtained which resemble those of the experiments rather well. It should be noted that the introduction by FFH of an overpotential, which is the difference between the electrode potential and the (pH-dependent) Flade potential, is nothing more than a reformulation of the problem and is in fact completely redundant. We prefer to work with the equation for [H+], because it is simpler. Detailed discussions of the FFH model can be found in the literature [4,5,18-221. The stability analysis of model I can be carried out without the discontinuity assumption. The stationary states h, and OS, of model I can be calculated explicitly: 1
%s = 1+
= -B+
h ss
(k_,/k,)h,
(11)
[B*+dhh(k-2/k2)11’* 2( k-,/k,)
where B = 1 - hb(k_*/k2) + k1(k_2/k2)/d’/2 but it is chemically irrelevant.)
(12) (There is another solution for h,
37
200
x
0’
”
(b)
(a)
Fig. 3. Two typical current-voltage characteristics for model I with ey = 0, ey = 20, d = 1, h, = 1, /~!~/k! = 0.1, I = q(l- 0&,, q = 100 and c any value: (a) /cp = 0.001; (b) kf = 0.01.
The stability of these stationary states can be investigated using linear stability analysis [5]. By superimposing small perturbations such that h = h,, + 6h and 8 = 0, + 619,we obtain two linear homogeneous differential equations -d
d”2k
-Ek,e,
1
-e( k, + k_,h,)
(13)
whose Jacobian matrix has the following eigenvalues 5: t2+tT+A=0
(14)
where T = d + c( k, + k_,h,,) A = de( k, + k_,h,)
+ cd1’2kIk2es
Because both T and A are always positive, the eigenvalues are strictly negative and consequently all perturbations die out; the stationary state is always asymptotically stable. The decay is oscillatory when 44 > T2, but no sustained oscillations are possible in model I. This result is in contrast with the results from the original FFH model and shows that the oscillations predicted by that model must be caused by the assumption that the passivation kinetics are discontinuous. Figure 3 shows two typical current-voltage characteristics which qualitatively resemble typical passivation characteristics. In Fig. 3(b) we see a mass-transferlimited dissolution due to the limited rate at which the H+ ions diffuse towards the electrode-electrolyte interface (the limiting current is proportional to dh,). For the rate constants, we have assumed Tafel laws according to k, = k; exp[ (e - e;)/2] k, = k; exp[ (e - e”,)/2] k_,=k”_,
exp[-(e-e;)/21
(1%
38
where e = FE/RT and E is the electrode potential. In Fig. 3 we took e’, = 0.0 for convenience. It should be noted that a reaction order of [H+] different from unity, say p, does not alter the number of steady states nor their stability properties. The expression for h, would simply be eqn. (112 to the power l/p. This is a specific property of model I and a direct result of assumptions (a)-(c), and it will not necessarily hold for the more elaborate models we shall consider below. DISSOLUTION-PRECIPITATION
MODEL WITH OHMIC DROP
Electrodissolution is usually accompanied by high electric currents and thus by appreciable ohmic drops in the electrolyte solution. It has been argued by several authors [16,24-281 that ohmic drop plays an important role in the interpretation of experimentally obtained current-voltage relationships. In many cases, a hysteresis which can be ascribed entirely to the potential drop between the working and reference electrode is observed [27]. Moreover, it is known that a negative faradaic impedance (dZ,/d E)-‘, in terms of which passivation can be understood (see Fig. 3), needs an ohmic series resistance to give rise effectively to instabilities such as hysteresis (or bistability, which is the same) or spontaneous current oscillations [24,29-311. Let R, be the ohmic resistance between working and reference electrode. If the potential on the electrode is identified with zero, the potential just outside the electrical double layer with E and the (controlled) potential of the reference electrode with V, eqn. (7) has to be modified according to I=-
V-E R,
=I,+z,=ACg
+z,
Using the dimensionless equations are obtained: de -= dt dh
potential
+q[k,(l-e)+k,(l+k_,Bh] 1 u-e
x=~(h,-h)-d1~2--(-)+d1’2[k2(1-8)-k_28h]
de
e = FE/RT, the following scaled differential
dt = +,(1-
(11)
e) - k,ehl
with LJ= FV/RT, r = there is no additional double-layer capacity referred to as “model quasi-steady state for
2DC,R,
and q = F2A’12/2 RTC,. The model assumes that ohmic drop associated with the passive layer and that the per unit surface area C, is a constant. This model will be II”. A useful version of model II is that corresponding to a 8 (i.e. dfZ/dt = O), which can then be equated to eqn. (11).
39
time
0
10
0
time
I”
(b)
(0)
Fig. 4. Typical current oscillation profiles from model II, with q = 100, ky = 0.01, ki = 0.01, k!!, = 0.001, r = 0.1, u = 30, d = 1, h, = 1, ef = 0, et = 20 and u = 30: (a) E = 200, (b) E = 01.
Numerical integration shows that it is indeed possible to obtain oscillatory solutions for model II (a Runge-Kutta routine with automatic step-size adaptor was used to solve the equations [32]). Figure 4 shows two examples, one in the quasi-steady-state assumption for 8. Although these waveforms do not resemble the experimental ones very closely, a more serious objection to model II is that h can become negative during these oscillations. This is also manifested in the oscillation profile in Fig. 4(b), in which a high current plateau exists at a higher current than the diffusion-limited current plateau in the steady-state curves of Fig. 3. This is caused by a negative H+ concentration. This artifact is a result of the unrealistic assumptions (a) and (b) discussed in the foregoing section. Clearly, a more realistic description becomes desirable. Therefore we end our study of model II at this stage and turn to a more complete description of the physical processes occurring during metal electrodissolution in the next section. IMPROVED
DISSOLUTION-PRECIPITATION
MODEL WITH OHMIC DROP
In this section models I and II are improved by replacing assumptions (a> and (b) with more realistic mass-transport descriptions. In this model all three ionic species (H+, Me+ and A-) migrate and carry their part of the current flowing through the cell, each in proportion to their transference number. Electroneutrality still prevails at every location x in the electrolyte solution:
WL + [Me% = [A-lx
(17)
which eliminates [A-] as a variable. The local transference j is given by
ti(x)= [H+],(h,++
[jldj AA-) + [Me’]J
hr,,le++ AA-)
number tj(x) of species
(18)
40
When Me+ is introduced as a variable in the model description, under conditions where the electrode is active, the Me+ concentration at the electrode-electrolyte interface may reach extremely high values. In reality, however, the solubility product K,, of Me+A- is exceeded, and a porous salt film precipitates onto the electrode. This has been observed, for instance, in the electrodissolution of iron in acidic sulphate electrolytes [33,34]. Then the rate at which the salt layer can dissolve into the neighboring solution determines the rate of the metal electrodissolution. Therefore formation of a salt layer may cause a limiting current in the current-voltage characteristic because the salt film dissolution is essentially independent of potential [35]. The exact kinetic equations for the formation of a salt film are complex [36]; presumably these kinetics are not essential for understanding oscillations because the salt film has no passivating properties. It is usually assumed that the passivating oxide layer forms underneath the salt film, closest to the electrode material. In order to account for the limited metal electrodissolution rate due to salt film formation, but at the same time to avoid complicated equations which may distract from the essence of the unstable behavior, we take the following opportunistic approach. We assume that a porous salt film has precipitated onto the electrode. To a first very crude approximation we assume that its conducting properties are no different from those of the electrolyte solution. This means there is no additional ohmic drop due to the salt layer. A metal ion that has come into being through reaction (1) travels through the salt film and reaches the salt film-electrolyte interface. Its probability of dissolving into the electrolyte is unity if the metal ion concentration in the neighboring solution is zero and zero when this concentration is beyond some critical value [Me+],,,. For intermediate situations, we assume a linear relationship to hold. This can be translated into a corrected rate constant k::
(19) where 0 I [Me+], I [Me+],,,. The critical concentration solubility product K,, of the MeA salt:
[Me+],,,
is related to the
[Me+lcrdA-l~= K,, Using eqn. (17), this is equivalent to
[Me+],,i, = -i[H’]o
+ (i[H+]i
+K,I)1’2
(21)
At this stage let us note that eqn. (5) can no longer be used to derive differential equations for [H+], and [Me+], since in this new model both t,+ and tMe+ are functions of the distance x from the electrode. There are two alternatives for deriving the new differential equations: (i) one can introduce an averaged transference number fj = (l/b)/ttj(x)dx and insert this in eqn. (5), or (ii) instead of assuming a linear concentration profile one can assume the concentration to be
41
uniform in a region of size 6 adjacent to the electrode and to attain its bulk value outside this region. The latter case would simply give eqn. (5) with the right-hand side divided by a factor of 2. In scaled terms, both approaches yield the same result except that the expressions for the various parameters differ by a factor of 2. Therefore, denoting [Me+], by u, and assuming for simplicity the quasi-steadystate assumption for 8, the following three scaled differential equations can be obtained: de v-e - = -qk,*(l -e> dt r dh dt = d( h, - h) - thd1’2 du dt-
- --LYdu + d”‘( 1 - 0)k:
(III) - t,d1’2
where a=_=-=-
D Me+
A Me+
A,-
D H+
A,+
A,+
h t, =
h( 1 +(Y) + 2~~4 (YU
t”=
h(l+cr)
+2au
u, = - $h + ( $h2 + K,$” and 8 follows from eqn. (11) with h for h,,. It should be noted that we disregard the salt film thickness in order to keep the model mathematically as simple as possible. In Figs. 5-10 we summarize some of the steady-state and dynamic properties of model III. In Fig. 5 the steady-state current-voltage curve is given for low ohmic resistance. The curve is stable and smooth without abrupt transitions or hysteresis phenomena. We have also simulated some linear sweep voltammograms for higher ohmic resistance. In the simulation, we took v = v” + ct. Results are depicted in Fig. 6. For fast sweep rates c, there is an abrupt transition at the active-passive transition; at lower sweep rates, this transition is accompanied by oscillations. The simulated sweep is in qualitative agreement with several experimental sweeps [3,16,17,231.
42
Fig. 5. Steady-state current-voltage relationship for model III at low ohmic resistance K,, = 16 and a = 0.2. Other parameter values are as in Fig. 4.
r = 0.02,
Oscillatory time series associated with a typical limit cycle solution of model III are shown in Fig. 7. The relaxation oscillation for the Z-f curve resembles that of the experiments. It should be noted from the 0-t curve that the electrode indeed effectively switches between an entirely active and an entirely passive state. For certain adjustments of the various parameters, the model predicts complex or even chaotic oscillations. Figure 8(a) shows a typical “mixed-mode” oscillation, which is made up of small and large excursions. As a function of one parameter, mixed-mode oscillations are known to follow complicated sequences but we shall not dwell upon them here (but see ref. 5). A chaotic solution of model III is illustrated in Fig. 8(b) with its Z-t behavior.
300
3oc
I
0 1"
5.
(a)
"
__
Yb
10
"
35
(b)
Fig. 6. Linear sweep voltammograms for r = 0.1: (a) sweep rate c = 1; (b) c = 0.04. Other parameter values are as in Fig. 5.
43
h
0
tme
0
50
time
5r
,L 0
50
(b
(a)
I
0
ti me
50
time
II
(c)
Fig. 7. Typical oscillatory time series for limit cycle solution of model III with u = 30, r = 0.1, h, = 1.2 and other parameter values as in Fig. 5; (a) Z-r behavior; (b) h-t behavior; (c) u-t behavior; (d) 13-t behavior.
Figures 9 and 10 illustrate the bifurcation behavior of model III in more detail. In Fig. 9 we show two typical parameter diagrams which map the stability of the steady-state solution(s) of model III onto the d-v and h,-u planes respectively. Within the curves labelled SN (saddle-node bifurcations), model III supports three simultaneous steady-state solutions. One of hem is always unstable. At the Hopf t
200
I
I/ time (b)
Fig. 8. Complex oscillations in model III; (a) mixed-mode oscillation for u = 30.52, r = 0.1 and d = 1; (b) chaotic oscillation for u = 24, r = 0.032505 and d = 0.5. Other parameter values are as in Fig. 5.
44
12
2 Hopf
0 “b
4
I
1
4
04
25
20
30
A”
35
0 20
25
_
v
30
35
(b)
(a) _
Fig. 9. Linear stability diagrams of model III which indicate regions of unique stable steady-state, oscillations and multiple steady states, and associated bifurcations (r = 0.1): (a) d-u diagram for h, = 1; (b) hb-u diagram for d = 1. Other parameter values are as in Fig. 5.
0
0.5
-
‘.o *
1.5
2.0
Ic)
Fig. 10. Bifurcation diagrams for model III: full curves, stable steady states or minimum and maximum of stable oscillation; broken curves, unstable steady states; Hopf, Hopf bifurcation; SN, saddle-node bifurcation; shaded region MMO, region where mixed-mode oscillations are obtained. (a) Z-h, diagram, u = 30, d = 1; (bl Z-u diagram, hi, = 1, d = 1; (cl Z-d diagram, u = 30, h, = 1.
45
bifurcation, a stable steady state becomes unstable and the system starts oscillating spontaneously [5]. Three representative sections of Fig. 9 are shown in Figs. 10(a)-10(c), where we have also indicated the occurrence of mixed-mode oscillations. These figures show the influence of the various parameters on the oscillation amplitude and also whether the transition from steady-state to oscillatory behavior is abrupt or smooth. For instance, with increasing bulk acidity, h,, the amplitude and the period of oscillation decreases, in agreement with the experimental results for iron in H,SO, [7]. (In short, the bifurcation diagrams of Fig. 10 were computed as follows. First, numerical integration provided an estimation of a steady-state expected to have a large basin of attraction. This result was used as an initial estimate for the Newton-Raphson method of iteratively solving the steady-state equations [32], and as a starting point for the one-parameter continuation of the steady-state. During the one-parameter continuation, the linear eigenvalues of the steady state were monitored and checked for passage through a saddle-node bifurcation (one zero eigenvalue) or a Hopf bifurcation (complex pair of eigenvalues with zero real part). Oscillation extrema were obtained from numerical integration. Figure 9 was drawn from an assembly of such one-parameter continuations. Usually, the electrode potential e was used as a parameter since this avoids instabilities in the NewtonRaphson algorithm close to saddle-node bifurcations.) DISCUSSION AND CONCLUSIONS
In this paper we have described a simple and general mathematical model to explain the current oscillations that are often observed at the active-passive transition in various metal electrodissolutions. The oscillations in the model arise from the coupling of a negative faradaic impedance due to the passivation kinetics with an uncompensated ohmic series resistance attributable to the electrolyte solution. It has also been shown that the oscillations in the.classical FFH model must be due to their discontinuous treatment of the passivation kinetics. When treated without assuming this discontinuity, this model can describe only stable steady-state behavior and never oscillations. Some authors have stated that discontinuous kinetics are a necessary prerequisite for obtaining oscillations. Apart from the objection that it is unclear how such discontinuous kinetics would arise physically, we believe that the present model shows that it is not necessary to incorporate discontinuous kinetics in order to obtain oscillations. Once again we emphasize that our model suffers from a number of oversimplifications and therefore it should not be judged on its quantitative properties. Although the experimental behavior is usually of low dimension (limit cycles or low-dimensional complex behavior), this only means that a low-dimensional qualitative model will suffice; a quantitative approach needs a far more detailed and laborious description of all the physicochemical processes in the system. In particular, the DP mechanism is surely not a reliable candidate for a quantitative
46
treatment of electrodissolution and passivation kinetics, but only provides the simplest possible interpretation of the negative faradaic impedance associated with the passivation phenomenon. For instance, the present model cannot explain the bistability and hysteresis that Epelboin et al. [37] observed in the current-voltage relationship for iron in H,SO, at low pH after correction for the ohmic drop. They were able to reproduce this observation qualitatively by assuming a more complex dissolution scheme which involves the formation of a non-passivating and a passivating intermediate. Such hysteresis can also be reproduced when an attracting mean field interaction between the adsorbing particles is introduced [38], i.e. a kind of Frumkin adsorption. When coupled to a second differential equation, Frumkin adsorption can give rise to limit cycle solutions (oscillations) without the explicit introduction of the electrode potential or mass transport (see, for example, the work by Talbot and Oriani [39], Brynn-Hilbert and Murphy [40] and Kado and Kunitomi [41]). Such an approach seems to us to be of a highly speculative nature since it often lacks experimental confirmation and it disregards many non-trivial processes appropriate to electrochemical systems. It should be noted that our model also differs in some important features from a recent model proposed by Haim et al. [42] who described periodic phenomena during anodic nickel dissolution. Their model applies to the transpassive region instead of the active-passive transition, it does not account for mass transport, which plays an important role in our model, and finally it describes potential oscillations under galvanostatic conditions whereas our model describes current oscillations under fixed potential conditions. Another point is that our treatment of the passive layer as a two-dimensional lattice gas that does not need a spatial description is certainly not always justified. The formation and dissolution of the passive layer is likely to give rise to islands of activation/passivation resulting in non-uniform current and concentration distributions. This would introduce a spatial component in the mathematical description of the system. Also, our simplification of the diffusion problem in terms of a linear diffusion layer misses an essential point. Under time-dependent conditions, such as spontaneous oscillations, the linear diffusion-layer approximation is no longer valid. There will be a kind of delay or relaxation associated with the diffusion process in trying to restore the linear or steady-state concentration profile, and a simple way in which this can be accounted for qualitatively is through the introduction of a second diffusion layer [43]. We observed that, with the introduction of this diffusion-layer relaxation, mixed-mode oscillations can appear where simple oscillations were obtained in the original model. This means that this effect can be very important in the interpretation of complex and chaotic oscillations in electrochemical systems. In fact, it has been shown for the indium thiocyanate electrochemical oscillator that, with the incorporation of diffusion-layer relaxation, excellent agreement with the experimentally observed complex bifurcation sequences is obtained [44]. Finally, throughout this paper we have assumed that the uncompensated ohmic cell resistance is entirely attributable to the electrolyte solution and could thus be
41
considered a constant. However, both the passive layer and the salt layer will contribute to the total ohmic drop in a non-negligible fashion, and it is not very likely that this contribution will remain unchanged during the oscillations. It also remains to be seen in what respect the complex and chaotic oscillations predicted by the model relate to the experimental chaos. Although they show some unusual similarities with the experimental results for a number of systems, this property of the model should not be overestimated. Many metal electrodissolutions display bifurcation scenarios that may still be traced to unique properties of these systems themselves. Nevertheless, despite these objections and shortcomings, the present model may prove to be a helpful tool in the interpretation and understanding of instabilities and (a)periodic oscillations at the active-passive transition in many metal electrodissolutions. REFERENCES 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
26 27 28 29
ES. Hedges and J.E. Myers, The Problem of Physico-chemical Periodic&y, Arnold, London, 1926. U.F. Franck, Z. Elektrochem., 62 (1958) 649. K-J. Vetter, Elektrochemische Kinetik, Springer, Berlin, 1961. J. Wojtowicz in J.O’M. Bockris and B.E. Conway (Eds.), Modern Aspects of Electrochemistry, Vol. 9, Plenum Press, New York, 1973, p. 47. SK. Scott, Chemical Chaos, Clarendon Press, Gxford, 1991. U.F. Franck and R. FitsHugh, Z. Elektrochem., 65 (1961) 156. J.J. Podesta, R.C.V. Piatti and A.J. Arvia, J. Electrochem. Sot., 126 (1979) 1363. G. Beni and S. Hackwood, J. Appl. Electrochem., 14 (1984) 623. C.B. Diem and J.L. Hudson, AIChE J., 33 (19871218. J.J. Podesta, R.C.V. Piatti and A.J. At-via, Electrcchim. Acta, 24 (1979) 633. F.N. Albahadily and M. Schell, J. Chem. Phys., 88 (1988) 4312. F.N. Albahadily, J. Ringland and M. Schell, J. Chem. Phys., 90 (1989) 813. M. Schell and F.N. Albahadily, J. Chem. Phys., 90 (1989) 822. M.R. Bassett and J.L. Hudson, J. Phys. Chem., 93 (1989) 2731. 0. Lev, A. WoIffberg, M. Sheintuch and L.M. Pismen, Chem. Eng. Sci., 43 (1988) 1339. D. Sazou and M. Pagitsas, J. Electroanal. Chem., 323 (19921 247. M.C.H. McKubre and D.D. Macdonald, J. Electrochem. Sot., 128 (1981) 524. J.B. Talbot, R.A. Oriani and M.J. Dicarlo, J. Electrochem. Sot., 132 (1985) 1545. SM. Sharland, Corros. Sci., 27 (1987) 289. A.J. Pearlstein and J.A. Johnson, J. Electrochem. Sot., 136 (1989) 1290. Y. Wang, J.L. Hudson and N.I. Jaeger, J. Electrochem. Sot., 137 (1990) 485. M. Pagitsas and D. Sazou, Electrochim. Acta, 36 (1991) 1301. N. Sato and G. Okamoto in JG’M. Bockris, B.E. Conway, E. Yeager and R.E. White (Eds.), Comprehensive Treatise of Electrochemistry, Vol. 4, Plenum Press, New York, 1984, p. 193. I. Epelboin, C. Gabrielli, M. Keddam. J.-C. Lestrade and H. Takenouti, J. Electrochem. Sot., 119 (19721 1632. I. Epelboin, C. Gabrielli and M. Keddam in J.O’M. Bockris, B.E. Conway, E. Yeager and R.E. White (Eds.), Comprehensive Treatise of Electrochemistry, Vol. 9, Plenum Press, New York, 1984, p. 151. C. Gabrielli, M. Keddam, C. Rakotomavo and H. Takenouti, J. Electrochem. Sot., 127 (1980) 2648. C.G. Law Jr. and J. Newman, J. Electrochem. Sot., 126 (1979) 2150. P. Russell and J. Newman, J. Electrochem. Sot., 134 (1987) 1051. H. Degn, Trans. Faraday Sot., 64 (1968) 1348.
48 30 R. de Levie, J. Electroanal. Chem., 25 (1970) 257. 31 M.T.M. Koper, Electrochim. Acta, 37 (1992) 1771. 32 W.H. Press, B.P. PIannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes, The Art of Scientific Computing, Cambridge University Press, Cambridge, 1989. 33 R. A&ire, D. Emsberger and T.R. Beck, J. Electrochem. Sot., 125 (1978) 1383. 34 M.E. Orazem and M.G. Miller, J. Electrochem. Sot., 134 (1987) 393. 35 L. Kiss, Kinetics of Electrochemical Metal Dissolution, Elsevier, Amsterdam, 1988. 36 P. Russell and J. Newman, J. Electrochem. Sot., 133 (1986) 59. 37 I. Epelboin, C. Gabrielli, M. Keddam and H. Takenouti, Electrochim. Acta, 20 (1975) 913. 38 G.L. Griffen, J. Electrochem. Sot., 131 (1984) 18. 39 J.B. Talbot and R.A. Oriani, Electrochim. Acta, 30 (1985) 1277. 40 D. Btynn Hilbert and S.V. Murphy, J. Electtochem. Sot., 138 (1991) L30. 41 T. Kado and N. Kunitomi, J. Electrochem. Sot., 138 (1991) 3312. 42 D. Haim, 0. Lev, L.M. Pismen and M. Sheintuch, J. Phys. Chem., 96 (1992) 2676. 43 M.T.M. Koper and P. Gaspard, J. Chem. Phys., 96 (1992) 7797. 44 M.T.M. Koper, P. Gaspard and J.H. Sluyters, J. Chem. Phys., in press.