Asymptotic methods are used to explore the role of one of the homogeneous chemical reactions, water protolysis, in a recent numerical model for the electrochemical pickling of steel. The analysis helps to interpret the somewhat surprising results of the original numerical model: that the local current density profile at the pickled strip does not depend on whether the water protolysis reaction is included in the model or not. The analysis also gives a useful systematic approach for determining qualitative estimates for the width of reaction boundary layers at cell electrodes, which can be of help when designing a computational mesh for the numerical solution of electrochemical models that have such layers.
Electrochemical pickling is a process that is used in the steelmaking industry to remove a surface oxide layer that is formed during the manufacture of thin steel strips. In the process, an example of which is shown in figure 1, a steel strip is passed between pairs of anodic and cathodic electrodes in an electrochemically neutral electrolyte, usually sodium sulphate, Na2SO4. When a current is passed through the cell, the ‘pickling’ reaction, i.e. the removal of the oxide layer, thought to consist predominantly of chromium oxide (Cr2O3), occurs according to 1.1 as do other electrochemical reactions that result in the evolution of oxygen and hydrogen, 1.2 and 1.3 respectively. Further information on many aspects of the pickling of austenitic stainless steels can be found in the survey by Li & Celis (2003) and the thesis by Ipek (2006).
The realization that the industrial process is not particularly efficient (Henriet 1996) has led to attempts at mathematical modelling in order to optimize it. Following earlier work that focused on laboratory situations with a non-moving strip (Ipek et al. 2002, 2006), Ipek et al. (2007) derived a model that accounts for the electrochemical reactions described above, as well as for the translational motion of the steel strip. Three variations of the model were derived, which either included or omitted some or all of the leading chemical reactions that are thought to occur in the electrolyte bulk, and the model equations were solved numerically using finite-element methods. While this gave qualitatively reasonable results, the drawback of a purely numerical approach for this application was the long computation time. This is due both to the size of the system of partial differential equations that must be simultaneously solved for, and the need to resolve the concentration boundary layers adjacent to reacting surfaces; in addition, artificial streamwise diffusion is necessary in order to maintain numerical stability, particularly as the strip velocity is increased. Furthermore, it should be noted that these drawbacks are already prevalent for the simplest of the three variants, which indicates that further meaningful extensions of this canonical model for electrolytic pickling, incorporating two-phase flow (Ipek et al. 2008) or further bulk reactions (Ipek 2006), will always be hindered by computational limitations.
To alleviate some of these difficulties, Vynnycky & Ipek (2009a) derived an asymptotic model for the simplest of the three variants mentioned above, and were able to demonstrate that it gave the same current density as the numerical models, although at considerably reduced computational cost. The variant considered was the one for which all bulk reactions are neglected. It is therefore of interest to see whether, and how, an asymptotic structure can be derived for the case when bulk reactions are not neglected. Although these are believed to be of importance for modelling the process, a puzzling feature of the results from the full numerical model of Ipek et al. (2007) was that the inclusion of the bulk reactions did not significantly affect the current density profiles. Furthermore, although Vynnycky & Ipek (2009b) presented preliminary considerations on analysing the asymptotic structure for one of the two models with bulk reactions, the results are far from complete.
In this context, a primary concern is whether the inclusion of reaction terms affects the asymptotic structure given by Vynnycky & Ipek (2009a). On the one hand, while electrochemical problems involving convection, diffusion and migration give rise to boundary layers whose thicknesses can be easily classified, it is clear that the inclusion of reaction terms necessitates a case-by-case approach. It is notable that while authors commonly refer to a thin reaction layer adjacent to an electrode, a qualitative estimate for its thickness is often not given (Yin et al. 1991; Nylén et al. 2007; Van Damme et al. 2009); the only available yardstick appears to be a rough approximation calculated from a formula derived by Levich (1962) for a pseudo-first-order homogeneous reaction (Van den Bossche et al. 1995, 1996).
The purpose of this paper is to extend the asymptotic model of Vynnycky & Ipek (2009a) to include one of the bulk reactions, water protolysis, that was present in the earlier numerical model of Ipek et al. (2007). In §2, we recapitulate the model of Ipek et al. (2007). In §3, we non-dimensionalize the governing equations, and propose an asymptotic structure for the model; an unexpected feature is the existence of a transition layer within the concentration boundary layers. Also, we demonstrate analytically why the current density is not affected by the inclusion of the water protolysis reaction in the model. In §4, the analysis is shown to be consistent with the earlier numerical solution of the problem. Conclusions are drawn in §5.
2. Model formulation
In this paper, we will, for simplicity, refer to reduced models (1) and (2) in Ipek et al. (2007) as models I and II, respectively. First, we recap the relevant model equations, which are for the concentrations, ci, of five ionic species (i=Cr2O72−, H+, Na+, OH−, SO42−) the electric potential in the electrolyte, Φ(e), and the electric potential in the steel strip, Φ(b). The model geometry is as given in figure 2 and consists of a two-dimensional region of length L(b) and width Deb+d; anode and cathode electrodes at x=0, each of length L, at a distance Die apart and held at electric potentials ±U/2, respectively; and a steel strip with electrical conductivity κ(b) that moves vertically upwards with speed V (b), and is at a horizontal distance Deb from the electrodes. Because of symmetry, only one half of the electrode-steel strip section is modelled; thus, d denotes the half-width of the steel strip, so that the centreline of the strip lies at x=Deb+d. The remainder of the boundary at x=0 is assumed to be impermeable and insulated. Electrolyte is assumed to enter at y=0 and exit at y=L(b). Characteristic values for d,Deb, Die, L, L(b), U, V (b) and κ(b) are given in table 1, as are other model parameters that are introduced in the main text later; a full list of symbols is given in the supplementary material.
(a) Governing equations
Three mechanisms, diffusion, migration and convection, are assumed to contribute to the transport of ionic species in the electrolyte. Assuming the applicability of transport equations in dilute solutions, see e.g. Newman & Thomas-Alyea (2004), the molar flux, Ni, of the ionic species i can be expressed via the Nernst–Planck equation as, for i=Cr2O72−, H+, Na+, OH−, SO42−, 2.1 where u is the hydrodynamic velocity of the electrolyte, Di is the diffusion coefficient and zi the charge number for species i. The quantities F (=96485 C mol−1), R (=8.314 J mol s−1) and T are the Faraday constant, the universal gas constant and the absolute temperature, respectively. The differential material balance for species i is given, in steady state, by 2.2 where Ri represents the production of species i through chemical reactions in the electrolyte. In this paper, the only homogeneous chemical reaction that we will consider will be the water protolysis reaction, 2.3
where and are the relevant rate constants for the forward and backward reactions, respectively. The expressions for Ri in equation (2.2) are then given by 2.4 and 2.5 In addition, the solution is assumed to be electrically neutral, which is expressed by 2.6 Equations (2.2) and (2.6) then provide a consistent description of transport processes in the dilute electrolyte. In addition, the current density, i, is calculated from the flux of charged species, and is given by Faraday’s Law as 2.7 Substituting equation (2.1) into equation (2.7), we can identify the electrolyte conductivity, κ(e), as 2.8 Also, u is taken to be the simplest possible profile that is consistent with the fixed electrodes and the translational motion of the steel strip; consequently, we take u=(u,v), where 2.9
For the potential field in the steel strip, conservation of charge gives 2.10
(b) Boundary conditions
At reacting surfaces, Faraday’s Law is used to relate local ionic fluxes to the current density. At all of these interfaces, electrochemical reactions are assumed to proceed according to Tafel kinetics, which seems reasonable in view of the high current densities (of the order of several kA m−2) at which the process is typically operated. In the following, iO2, iH2 and iCr2O3 are the current densities (A m−2), i0,O2 and i0,H2 are the exchange current densities (A m−2) for reactions (1.2) and (1.3), respectively, αO2 and αH2 are the transfer coefficients and and are the equilibrium potentials. Values for αO2, αH2, i0,O2, i0,H2 were determined from experimental polarization curves; see table 2 and Ipek et al. (2005) for details. In addition, and as already discussed at length in earlier work (Ipek et al. 2005, Ipek 2006), experimental results were used to establish the fraction, Λ, of the total current that is attributable to reaction (1.1); as in Vynnycky & Ipek (2009a), we set Λ=0.1.
Taking into account the oxygen evolution reaction (1.2) gives 2.11 2.12 2.13 2.14 and 2.15 where 2.16 with where pHeq, with pH.
From the hydrogen evolution reaction (1.3), we have 2.17 2.18 2.19} 2.20 and 2.21 where 2.22 with where pHeq.
From the oxygen and hydrogen evolution reactions (1.2) and (1.3), as well as chromium oxide removal according to equation (1.1), we obtain 2.23 2.24 2.25 2.26 and 2.27 where 2.28 2.29 and 2.30 with Continuity of the normal current at the steel strip/electrolyte interface gives the boundary condition 2.31 where 2.32
At the inlet at y=0, the concentrations must be such that there is local equilibrium, which implies that 2.33 In practice, one would expect to be able to measure cNa2SO4 in the homogeneous electrolyte entering the domain. We obtain immediately that 2.34 moreover, we take The electroneutrality condition gives 2.35 which reduces to , giving 2.36 as in the case when there are no bulk reactions.
The respective values for , , , and are as given in table 3, and the boundary condition is written as 2.37
At the outlet boundary, we assume that the convective flux in the axial direction of the channel, for all ionic species, is dominant. Consequently, we set the sum of the diffusion and migration fluxes to zero, i.e. for , H+, Na+, OH−, , 2.38 Combining this with the electroneutrality condition, we obtain in simpler form 2.39
(vi) Insulated and symmetry boundaries
At the boundary of the electrolyte domain that is assumed insulated, i.e. x=0, L≤y≤L+Die, we have, for , H+, Na+, OH−, 2.40 Similarly, at the boundaries of the steel strip, which are either symmetry planes or insulated, we have 2.41
With a view to determining the correct asymptotic structure for this problem, we begin by non-dimensionalizing with 3.1 In what follows, the solution for the strip is exactly the same as that given in Vynnycky & Ipek (2009a), so we omit that here and focus solely on the electrolyte.
Substituting equation (3.1) into equation (2.2) gives 3.2 where 3.3 3.4 where and . Noting from table 1 that and hence that we observe that 3.5 this motivates investigating the asymptotic regime for which Pe≫1, Π≫1, ϵ≪1, λ1≪1. In particular, since Pe≫1, we can expect boundary layers, and we will need to consider these and the bulk separately. To do this, we proceed as in Vynnycky & Ipek (2009a) by considering asymptotic expansions in the bulk of the form 3.6 and 3.7 As for the case with no bulk reactions, 3.8 where, at leading order in ϵ, 3.9 and 3.10 in effect, the asymptotic expansions for these minority ions begin at O(ϵ). Consequently, equation (3.2) reduces to 3.11 In what follows, it will be convenient to denote the leading-order solution for ϕ(e) in the bulk as ϕbulk.
We consider the anode, the cathodic portion of the strip, the cathode and the anodic portion of the strip in turn; the complexity of the analysis increases in this order, with preliminary results for the anode being necessary for use later.
Consider first the material balance equations for H+ and OH− ions. In dimensionless form and neglecting the second derivatives in Y , we have 3.12 where is related to ϕ(e) by If there is now to be a balance at leading order between reaction and diffusion terms in equation (3.12), then we should have, for H+, 3.13 where [XH+] denotes the thickness of the proposed reaction layer for H+, and where Hence, 3.14 On this length scale, the magnitude of the convective term is so that the reactive–diffusive balance is consistent if For OH− ions, on the other hand, if the reaction terms are to be of significance, we must have 3.15 where [XOH−] denotes the thickness of the proposed reaction layer for OH−, and denotes the order of magnitude of COH− in this layer; we would then need Since OH− ions are not produced at the anode, we can expect at most COH−∼ϵ; hence, this inequality will be satisfied for the parameters given in equation (3.5). Nevertheless, it is instructive to see how the analysis differs for the three cases:
If we have, from equation (3.12) for i=H+ and OH− at leading order, 3.16 and 3.17 where and Hence, we can determine CH+ first and use this to find from equation (3.17). However, this expression for will clearly not satisfy boundary condition (2.14) at X=0, for which a nested sub-layer of thickness ϵ(λ1Pe)−1/2 is required. Setting we have 3.18 where in addition, we must have the matching condition 3.19 Notice that the profile for is entirely unaffected by this sub-layer.
If we have, at leading order, 3.20 and 3.21 Hence, for the most part, the reaction term still has no effect at leading order in the equation for , although for a different reason from that in (i): the equation for COH− holds the reaction at equilibrium. Once again, a sub-layer of thickness ϵ(λ1Pe)−1/2 is necessary for which satisfies equation (3.18). However, equation (3.20), with the reaction term removed, cannot be valid all the way down to since equation (3.21) ceases to be valid in the sub-layer; hence, at first sight, boundary condition (2.12) cannot be applied directly. Consequently, this means that, unlike in case (i), there will be a coupled problem to solve for and in the sub-layer. Thus, instead of equation (3.20), we will have 3.22 which at leading order gives just 3.23 if Hence, will be conserved through the sub-layer, meaning that equation (2.12) can in fact be applied to equation (3.20), with equation (3.21) being used to eliminate the reaction term. In consequence, the main difference between cases (i) and (ii) is that in equation (3.18) must be replaced by CH+, which is solved for using equation (3.23).
If the details are more or less the same as for (ii), although (3.20) and (3.21) are replaced by 3.24 and 3.25 where This time, a sub-layer of thickness ϵ(λ1Pe)−1/2 is required, in order to complete the details for
In summary, we have shown that the reaction term has no effect at all in any of the cases on the conservation equation for H+ at leading order; hence, it is the same as it was when the bulk reaction was omitted. A further consequence of this is that we will obtain, after some manipulation as in Vynnycky & Ipek (2009a), 3.26 where 3.27
(b) Boundary layer at cathodic strip
The material balance equations for H+ and OH− ions are, in dimensionless form and neglecting the second derivatives in Y , 3.28 where is related to ϕ(e) by If there is now to be a balance at leading order between reaction and diffusion terms, then we should have, for OH−, [XOH−]∼(λ1Pe)−1/2, where [XOH−] denotes the thickness of the proposed reaction layer for OH−, and Since the magnitude of the convective term is O(1), the reactive–diffusive balance is consistent only if λ1≫1; thus, one would not expect to see a reaction layer at all at this boundary for COH− for the parameters given in equation (3.5). Consequently, the reaction terms in equation (3.28) are negligible at leading order, and the situation closest resembles case (i) in §3a. Furthermore, equation (3.28) can satisfy the relevant boundary condition, equation (3.26), at X=1, since it is not affected by the nested sub-layer for the ion that participates in the bulk reaction, H+, which will nevertheless be necessary and which we discuss next.
For 1−X∼Pe−1/2, we have 3.29 where Since ϵ2≪λ1, the leading order balance gives just Naturally, this solution cannot satisfy boundary condition (2.24), and a nested sub-layer in is required. Setting we obtain 3.30 once again, the left-hand side can be neglected, since ϵ2≪λ1. As at the anode, the decoupling of the problem for H+ at this layer from the leading-order problem for OH− is evident.
The remaining details are straightforward. The cases analogous to (i)–(iii) in §3a for the anode are determined by the magnitude of the parameter λ1, as opposed to we omit the details here. However, here, as the anode, the conservation equation for the current-carrying ion, i.e. OH−, is at leading order the same as it was when the bulk reaction was omitted. Again, as in Vynnycky & Ipek (2009a), we obtain 3.31 where ϕ(b)bulk=ϕ(b)−ϕbulk(1,Y) and 3.32
The situation at the cathode is more complicated than that at the anode, because of the advection of H+ ions from the anode. Hence, in model II, in addition to having COH−∼O(1), as one would expect since OH− ions are being produced at the cathode, we also have CH+∼O(1). In view of the added complexity, we consider only the case that is relevant to the computations, so that ; thus, the concentration boundary layer will be of thickness Pe−1/3. Once again, setting we have 3.33 and 3.34 The numerical results to be presented later in §4 suggest that COH−∼O(ϵ2) when but CH+∼O(ϵ2), COH−∼O(1) near the electrode, and that there must be a transition layer between these two regions. For we set and obtain from equations (3.33) and (3.34) at leading order 3.35 and 3.36 For the region near the electrode, we note that COH− undergoes a rapid change in boundary condition from equation (2.40) to equation (2.20) at i.e. the leading edge of the cathode. Although the mathematical implications of this may be undocumented within the electrochemical literature, it is already well established in the viscous boundary-layer theory of fluid mechanics that a boundary layer develops a two-tier structure in order to accommodate a change in boundary condition (Smith 1974; Veldman 1975; Cebeci et al. 1979; Vynnycky 1998). The location of the interface between the two tiers is similar to that considered recently by Vynnycky & Borg (2010) for the forced convection of excess supporting electrolyte between parallel electrodes; this implies that there is sub-layer for which where Setting 3.37 equations (3.33) and (3.34) become, at leading order, 3.38 and 3.39 a further point concerns the lower tier, which will have a sub-layer of thickness ϵ(λ1Pe)−1/2 for although the details are similar to those already given for at the anode, and we therefore do not repeat them here.
From the above, we expect the transition layer to be where so we set where χ≪1 and will have to be determined. For the transition layer, we will have so that, on setting equations (3.33) and (3.34) become 3.40 and 3.41 where φc∼O(1) is a strictly positive function of ζ; more specifically, and must be strictly positive since the electrical potential increases from its value at the cathode as increases. From equations (3.40) and (3.41), it is clear that we require leading to 3.42 and 3.43 The boundary conditions are 3.44 3.45 3.46 and 3.47 where and are positive O(1) functions that are to determined; we show how to do so presently. From equations (3.42) and (3.43), we see that 3.48 where is a function of ζ; from equations (3.44)–(3.47), 3.49 and hence must be negative.
To demonstrate that this analysis is self-consistent, it would be beneficial to provide a solution to equations (3.42)–(3.47); however, there appears to be no analytical solution, and the equations must be solved numerically. Prior to doing this, we consider the behaviour of and as As try where A+,B+ and m+ are all functions of ζ; also, we must have m+<0. Then, 3.50 and 3.51 thus So, since clearly we can have m+<0 if we take the negative sign. Note that equation (3.50) the expression for m+, whereas equation (3.51) provides the ratio A+/B+. On the other hand, as we set where A−,B− and m− are all functions of ζ; this leads to as required, and As confirmation of the feasibility of this asymptotic structure, we provide a numerical solution to equations (3.42)–(3.47), which is given in figure 3. Note, however, that this figure is only illustrative, in the sense that we have taken arbitrary values for φc and —in this case, φc=1 and In practice, neither will be known a priori, but will have to be solved for. We return in §3e to the issue of how to determine φc. on the other hand, is found by solving equation (3.39) for COH− in the lower tier and equation (3.35) for CH+ in the upper tier, subject to two interface conditions at The first comes from equation (3.49), although now expressed in upper and lower tier variables, rather than transition layer variables: 3.52 The second comes from the continuity of current, which can be reduced to 3.53 this makes use of equation (2.7) and the fact that N and are both continuous across the transition layer.
Thus far, we have demonstrated a self-consistent asymptotic structure for the boundary layer at the cathode. However, it is evident that the reaction terms are not negligible for this boundary layer, as they were shown to be at the anode and at the cathodic portion of the steel strip, and it is not therefore immediately obvious whether manipulations similar to those that lead to equations (3.26) and (3.31) can once again be used to yield a boundary condition in terms of ϕbulk only. Starting with the governing equations at leading order for this layer, i.e. 3.54 and 3.55 we obtain 3.56 which is the same as that derived by Vynnycky & Ipek (2009a) when bulk reactions were not accounted for; the reason for this coincidence is that and Integrating once with respect to gives 3.57 where 3.58 with 3.59 Hence, for i=H+,Na+,SO42−, we arrive at 3.60 where 3.61 and the boundary conditions in equations (2.18), (2.19) and (2.21) can be expressed as 3.62 Thence, analogous to equations (3.26) and (3.31), we obtain, as desired, 3.63
(d) Boundary layer at anodic strip
Setting we have, analogous to equations (3.33) and (3.34), 3.64 and 3.65 The numerical results to be presented later in §4 suggest that COH−∼O(ϵ2) near the strip, but CH+∼O(ϵ2), COH−∼O(1) when , and that there must be a transition layer between these two regions. For we set and obtain from equations (3.33) and (3.34) at leading order 3.66 and 3.67 For the region near the strip, we note that CH+ undergoes a rapid change in boundary condition from equation (2.40) to equation (2.20) at where if δ≪1; by analogy with the cathode, this is where the local current starts to increase rapidly, although its location will not be known a priori in general. Once again, a two-tier structure is required, although this time where . Setting 3.68 equations (3.33) and (3.34) become, at leading order, 3.69 and 3.70
From the above, we expect the transition layer to be where so we set where χ≪1 and will have to be determined. For the transition layer, we will have so that, on setting equations (3.33) and (3.34) become 3.71 and 3.72 where φb,a∼O(1) is a strictly negative function of ζ; more specifically, and must be strictly negative since the electric potential decreases from its value at the steel strip as increases. From equations (3.71) and (3.72), it is clear that we require the details thereafter are exactly the same as in equations (3.42)–(3.49), although with the roles of and reversed, and so we do not repeat them here. Once again, the thickness of the transition layer is found to be ϵ(λ1Pe)−1/2, as for the cathode, although it has been arrived at differently.
Lastly, we consider whether a development analogous to that in equations (3.54)–(3.63) is possible. Starting with 3.73 and 3.74 we obtain 3.75 which can then be integrated to give 3.76 where 3.77 and 3.78 with 3.79 once again, this is the same as when bulk reactions were not accounted for. Finally, we arrive at, for 3.80 where 3.81 At we will have simply 3.82 Matching to the bulk then gives 3.83
Figure 4 summarizes the various boundary layers, transition layers and sub-layers that the analysis has uncovered for model II.
As regards a systematic strategy for solving the asymptotic problem, this turns out to be the same as that given by Vynnycky & Ipek (2009a), when water protolysis was not included: ϕbulk and ϕ(b) can be solved for ahead of solving for the concentrations in the boundary layers. In particular, note also that once ϕbulk becomes known, it is possible to obtain the values of φc and φb,a required to completely determine the transition layers in §3c,d.
4. Comparison of asymptotic and numerical results
First, as a reference for the results of model II, we show in figures 5 and 6 the profiles for CH+ and COH−, respectively, at Y =0.25 and 0.75, as computed by Ipek et al. (2007) for model I. Focusing first on the curves for Y =0.25, it is evident that the concentration boundary layer near X=0 is much thicker than that at X=1. Also, CH+ is much higher than the bulk value near X=0, whereas COH− is much higher than its bulk value near X=1; this is to be expected in view of the electrochemical reactions that occur at the anode and the cathodic portion of the strip, respectively. Similar trends are observed for the curves at Y =0.75, although with the roles of H+ and OH− reversed. However, the additional features are that CH+ and COH− are high near X=0 and X=1, respectively, even though H+ and OH− ions are not being produced in the electrode reactions at the cathode and the anodic portion of the strip, respectively. This is because H+ and OH− ions are advected up from the anode and cathodic strip, respectively.
Figures 7 and 8 show the profiles for CH+ and COH−, respectively, at Y =0.25 and 0.75, as computed for model II in Ipek et al. (2007). For these computations, so we would expect the thickness of the boundary layer at Y =0.25 for CH+ near X=0, and for COH− near X=1, in model II to be of the same order of magnitude as that for model I. In both cases, the concentration of the ion that does not participate in the electrochemical reaction at that boundary is considerably smaller than the bulk value; this supports the scalings used in §3a,b. On the other hand, the curves for model II for Y =0.75 display the rapid transition in CH+ and COH− within each boundary layer; once again, these figures support the scalings used in §3c,d. It is, however, unlikely that the numerical results can resolve the details within the transition layer, nor the details of the sub-layers for the minority non-reacting ions: the analysis indicates that these should be of dimensionless thickness 10−8, which is considerably smaller than any length scale that Ipek et al. (2007) had envisaged that they would need to resolve. However, the analysis also shows why there was no actual need to resolve these length scales, as regards determining the current density correctly: the leading-order problem for the electric potential, and hence the current density, decouples from that for the concentrations.
This paper has used asymptotic methods to analyse the role of the water protolysis bulk reaction in a recent numerical model for the electrochemical pickling of steel (Ipek et al. 2007). Starting from an asymptotic model for the problem in the absence of bulk reactions (Vynnycky & Ipek 2009a), significant analytical progress was once again possible. In particular, our asymptotic results were able to explain why the current density profiles obtained in the numerical model for the cases with and without the water protolysis reaction were almost identical. Furthermore, we have provided a systematic approach for the mathematical treatment of reaction layers that may be of use in more general electrochemical systems.
In the context of the modelling of electrochemical pickling, the next step will to be to include the disassociation of sulphuric acid, a further bulk reaction that was initially thought to be of importance in the earlier numerical model (Ipek et al. 2007); however, the inclusion of that reaction was found not to affect the results significantly either, and further asymptotics should be able to explain the reason for this also.
The first author acknowledges the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland Mathematics Initiative grant 06/MI/005.
- Received January 1, 2011.
- Accepted February 22, 2011.
- This journal is © 2011 The Royal Society