## Abstract

In this paper, we address the phenomenon of temporal, self-sustained oscillations which have been observed under quite general conditions in solid oxide fuel cells. Our objective is to uncover the fundamental mechanisms giving rise to the observed oscillations. To this end, we develop a model based on the fundamental chemical kinetics and transfer processes which take place within the fuel cell. This leads to a three-dimensional dynamical system, which, under typical operating conditions, is rationally reducible to a planar dynamical system. The structural dynamics of the planar dynamical system are studied in detail. Self-sustained oscillations are shown to arise through Hopf bifurcations in this planar dynamical system, and the key parameter ranges for the occurrence of such oscillations are identified.

## 1. Introduction

Solid oxide fuel cells (SOFCs) are electrochemical devices which generate electricity by continuously converting hydrogen to water, without the need for combustion. The device can be constructed simply from five components, namely the anode, cathode, electrolyte and two interconnect wires, making it ideal for generating electricity when compared with conventional, fossil fuel-burning power technologies, which tend to have many moving parts. In addition, an SOFC produces high-grade waste heat during regular operation which is a desirable quality in many applications. For larger power applications, numerous SOFCs may be stacked together in series using interconnect materials to achieve the desired output. Such interconnecting materials are often also designed to accommodate other functions such as air and fuel channelling.

We model the operation of an SOFC by the two half cell reactions which occur at opposite electrodes, separated by a solid electrolyte which is usually made from yttria-stabilized zirconia. The reactions are, respectively,
1.1
and
1.2
The electrolyte ideally acts as a barrier to anything other than the ions produced at the cathode via the half cell reaction (1.2), at temperatures approximately between 600 and 1000^{°}*C*. Once the O^{2−} ions have migrated across, the electrolyte they react, via (1.1), at what is known as the triple phase boundary (TPB) within the anode, where the electrode, electrocatalyst and reactants/products all meet. The electrons liberated by this reaction generate current by flowing around an external circuit back to the cathode, where they are free to participate in reaction (1.2) again. The overall lumped reaction for an SOFC is given as,
1.3
It has been widely observed and reported that electrochemical systems may exhibit oscillatory dynamics, and SOFCs have been shown to produce both current oscillations in potentiostatic mode [1], and voltage oscillations in galvanostatic mode [2]. These electrical fluctuations are undesirable in applications, which require a steady power supply. Therefore, to avoid such behaviour, it is essential to identify and understand the fundamental mechanisms giving rise to oscillatory dynamics, and to identify the key parameters and parameter regimes, associated with this behaviour. It is this question that we address in detail in this paper. We develop a rational model based on the fundamental chemical kinetics and transfer processes operating in the SOFC. This leads to a three-dimensional dynamical system. The dimensionless form of this dynamical system allows for the identification of the key dimensionless parameters which govern the behaviour in the SOFC. At typical operating conditions, we are able to rationally reduce the three-dimensional dynamical system to a planar dynamical system. We then perform a detailed bifurcation analysis of this planar dynamical system, which reveals the occurrence of equilibrium saddle–node, Hopf, homoclinic and periodic saddle–node bifurcations. The Hopf, homoclinic and periodic saddle–node bifurcations give rise to self-sustained periodic dynamics. We identify the key parameters associated with these bifurcations, and the regions in parameter space where self-sustained oscillatory dynamics occur. Although a broad range of experimental data on SOFC oscillatory dynamics is, at present, unavailable, the predictions established in this paper appear to be broadly in line with, for example, results given in [3–5], which include bifurcations, multiple steady states and oscillatory phenomena. This is encouraging in underpinning the model we have developed here as capturing the fundamentals of oscillatory dynamics in SOFCs. Specifically, the predictions of this model have enabled us to inform, design and establish a structured, experimental programme in the EPSRC-supported Doctoral Training Centre for Hydrogen Fuel Cells and their Applications, which allows us to control the model parameters using various operating conditions, in order to induce oscillations from an SOFC. A significant point is that the model developed here is based on simple and fundamental physical and chemical mechanisms, enabling each of the dimensionless parameters in the model to be readily determined and interpreted, in any particular experimental realization. This enables direct and clear comparison of model predictions and experimental data. The early stage results from this programme are encouraging in supporting the predictions of the model and are enabling a deeper understanding of the mechanisms and key parameters associated with SOFC oscillatory dynamics to be established.

The paper is structured as follows. In §2, we introduce and develop the model, together with its dimensionless form and its rational reduction. The significant dimensionless parameters in the model are identified. Section 3 is devoted to a detailed bifurcation study of the planar dynamical system modelling the SOFC dynamics. Particular attention is paid to the parameter regimes in which the dynamics give rise to self-sustained oscillations. In §4, we give complimentary numerical integrations of the dynamical system for the distinct parameter regimes determined in §3. In §5, we draw conclusions and make suggestions for further developments. A brief of symbols and notations is provided in appendix A.

## 2. The model

Here, we develop a model describing the dynamics of an SOFC. The model is based on fundamental chemical kinetics coupled with Fickian mass transfer through a porous medium. The fuel used will be methane, and the model will account for the principle electrochemical and chemical reactions which take place at the anode. The oxidation of hydrogen, which contributes to the electrical current, occurs at the TPB within the anode, and the other chemical reactions are promoted by the anode catalyst at the surface. The reactions, which take place within an SOFC, can be widely varied and include the electrochemical oxidation of carbon monoxide, which contributes current as well as direct methane oxidation. However, the direct oxidation of methane will be neglected in this model, because, for most methane-fuelled SOFCs, the behaviour is dominated by the reaction scheme set out in (2.1)–(2.4) [6]. The fundamental reaction scheme adopted is as follows:
2.1
2.2
2.3
2.4
Here, *k*_{i}, (*i*=1,2,3,4) are the reaction rate constants for each respective reaction step. As the operating temperature range of an SOFC is approximately 600–1000^{°}*C* [6,7], all the reactants will be in the gas phase. This means that the water involved in the reactions will be in the form of steam, and there will be no water management issues to resolve which are commonly associated with the proton exchange membrane fuel cells (PEM) [8]. This is due to the negligible effect that hydration has on an SOFC electrolyte’s conductivity, when compared with that of PEM, which often requires careful management in both the gas and liquid phase.

The reactions involving *H*_{2}O are chemical reactions representing the endothermic steam reforming of methane (2.1), and the slightly exothermic water–gas shift (WGS) reaction (2.2). Reactions (2.3) and (2.4) determine the electrochemical oxidation of hydrogen and carbon monoxide, respectively. In the electrochemical reactions, electrons are liberated, then the current is collected and distributed along an external circuit.

Generally, in an SOFC, a variety of hydrocarbon fuels can be used although the shorter the hydrocarbon chain, the less carbon coking will occur. The general steam reforming reaction for hydrocarbons is 2.5 In the absence of steam, the direct oxidation of methane may occur [6], via the reaction, 2.6 or more generally, 2.7 This reaction leads to the build-up of solid carbon deposits on the anode, which then reduces the reaction surface area, and hence degrades the performance. In this work, we have included steam in the fuel stream, and adopted (2.1) as the steam reforming reaction. Future work may include the reaction given by (2.6).

A further simplification can be made to the reaction scheme (2.1)–(2.4), because the reaction given by (2.2) proceeds much faster than the other reactions. In general, for reactions (2.1)–(2.4) in SOFCs, it has been confirmed [6,9,10] that
2.8
which allows us to reduce the reaction scheme (2.1)–(2.4) to
2.9
and
2.10
with the composite reaction (2.9) being governed by the slowest component rate of reaction. This reduction is supported by Ho *et al.* [11], who note that the reaction step (2.2) reaches equilibrium because it is kinetically fast and almost all of the CO is consumed in this reaction. Any remaining CO may participate in the reaction given by (2.4) which contributes electrical current; however, the CO oxidation rate is around two to three times slower than that of hydrogen oxidation. Hence, the dominant current contribution is from hydrogen oxidation alone. This is confirmed by Yakabe *et al.* [12], who found that the WGS reaction (2.2) was fast enough to significantly reduce the concentration polarization downstream of the fuel inlet.

Specifically, when calculating electrical current, the accuracy of the model will be increased by considering the electrochemical oxidation of CO in addition to the hydrogen oxidation. However, in the present model, we aim primarily to capture the qualitative characteristics of the most predominant reaction kinetics of an SOFC in order to investigate the principle dynamic behaviour, and so the exclusion of reaction (2.4) is justified on the grounds outlined above. Henceforth, the model will be based on the reaction scheme (2.9) and (2.10).

We now introduce the chemical concentrations of the reactant species as
2.11
The configuration within which the reaction scheme (2.9) and (2.10) is embedded is illustrated in figure 1. There are two principle zones. In zone B, the species CH_{4}, H_{2} and H_{2}O (steam) are present and advected through the system from the inlet to the outlet and, in particular, we note that oxygen is available in abundance from the oxidant stream through the porous cathode, so that, throughout, we may regard the oxygen concentration *x* as fixed at the oxidant stream concentration. The inlet flow contains only CH_{4} and H_{2}O, at inlet concentrations *a*_{0} and *c*_{0}, respectively. We restrict attention to the situation when the flow rate *q* is sufficiently rapid, relative to the Fickian transfer rates between zone B and zone A, so that zone B provides a steady transfer pool of CH_{4} and H_{2}O, at the inlet concentration, for Fickian transfer into zone A. This is a reasonable assumption as at higher velocities the depletion zone of a reactant stream flowing over a catalytic surface becomes very small in comparison with the width of the channel [13]. The main reactions, given by the reduced scheme (2.9) and (2.10), take place in zone A. Within zone A, O_{2} is present at the pool concentration *x*, whereas CH_{4}, H_{2}O and H_{2} are transferred across the boundary between zone A and zone B via Fickian transfer. Zone A represents the porous anode where the reactions take place. Zone B is the fuel stream channel. We need to only consider the reaction rate equations in region *A* because this is where the reactions occur. The associated reaction rate and transfer balances for *a*, *b* and *c* lead, via (2.9) and (2.10), to the three coupled, nonlinear ordinary differential equations, namely
2.12
2.13
2.14
Here, *D*_{a}, *D*_{b} and *D*_{c} are respectively, the effective diffusion coefficients for the Fickian transfer of methane, hydrogen and water (as steam) to the reaction site, and *h* is the scale thickness of the anode. The coefficient *A*/*V* is the surface area of transfer-to-volume ratio, whereas *a*_{0} and *c*_{0} represent the inlet concentrations of species *a* and *c*, respectively. The first term of each rate equation represents the transfer of the particular species into the reaction zone A, whereas the remaining terms represent the production/consumption of the species as given by the reduced reaction scheme (2.9) and (2.10).

We now establish typical scales for each variable in equations (2.12)–(2.14). Let *a*∼*a*_{s}, *b*∼*b*_{s}, *c*∼*c*_{s}, *t*∼*T* with *a*_{s},*b*_{s},*c*_{s} and *T* being typical scales for *a*,*b*,*c* and *t*, respectively. A structured balance of terms in equations (2.12)–(2.14) then requires,
2.15
This results in the following dimensional scales
2.16
We non-dimensionalize the equations (2.12)–(2.14) with respect to the scales we have determined in (2.16). We write *a*=*a*_{s}*a*^{′}, *b*=*b*_{s}*b*^{′}, *c*=*c*_{s}*c*^{′}, *t*=*Tt*^{′}. Equations (2.12)–(2.14) then become
2.17
2.18
2.19
where primes are dropped for convenience.

In equations (2.17)–(2.19), we have introduced the following dimensionless parameters.
The parameter *ε* measures the ratio of reaction rates for the reaction scheme (2.9) and (2.10). For *ε*≫1, reaction (2.10) is the faster, whereas for *ε*≪1, reaction (2.9) is the faster. The parameter measures the ratio of the diffusivity of methane to the diffusivity of water into the reaction zone A. Similarly, the parameter measures the ratio of the diffusivity of hydrogen to the diffusivity of water into the reaction zone A. The parameters and represent the ratios of inlet concentration to the scaled concentration of methane and water, respectively.

Experimental studies in SOFCs confirm that, in general, 0<*ε*≪1 which corresponds to reaction (2.9) being significantly faster than reaction (2.10). Specifically, this is confirmed by the following two results. Park *et al.* [10] found through simulation that there is accumulation of water within the cell owing to the faster consumption of water by the reforming reactions than being consumed by the electrochemical oxidation of hydrogen. This is supported by You *et al.* [14], who found experimentally that there was hydrogen in the anode exit stream indicating that the reforming reactions were producing more hydrogen than was able to be oxidized. These observations support the condition that 0<*ε*≪1.

In what follows, we consider 0<*ε*≪1 with . This is supported by typical values in SOFC experiments [15–17], which lead to order of magnitudes,
Under these conditions, equations (2.12)–(2.14) can be rationally reduced. At leading order in *ε*, equations (2.17)–(2.19) become,
2.20
2.21
2.22
From equation (2.22), it follows immediately that
2.23
It is straightforward to establish that, when 0<*ε*≪1, the manifold (2.23) is attracting in the invariant quadrant *a*,*b*,*c*≥0 for the dynamical system (2.19). On substituting for *c* from (2.23) into (2.20) and (2.21), we obtain
2.24
and
2.25

Thus, we have arrived at a two-dimensional dynamical system for (*a*(*t*),*b*(*t*)) in *t*≥0 which models the temporal dynamics of the concentrations *a*(*t*) and *b*(*t*) in the anode of an SOFC, with the concentration *c*(*t*) then following via (2.23). In this study, we analyse the situation when which corresponds to a low steam concentration in the inlet flow. The leading order form of equations (2.24) and (2.25) is then,
2.26
and
2.27
which determine the temporal dynamics of the concentrations (*a*(*t*),*b*(*t*)). We now consider in detail the structure of the planar dynamical system (2.26)–(2.27), which we refer to henceforth as [D–S].

## 3. The phase portrait of [D–S]

Here, we establish some general properties of [D–S] in the (*a*,*b*) phase plane which eventually enable us to construct the global phase portrait of [D–S]. We examine the phase portrait of [D–S] in the quadrant , where
3.1
It is convenient to write [D–S] as
where ** x**=(

*a*,

*b*) and, 3.2 and 3.3

### (a) Invariant and attracting sets

We first observe that in ,
3.4
where ** i** and

**are unit vectors pointing in the positive**

*j**a*- and

*b*-directions, respectively. We conclude from (3.4) that is a positively invariant region for [D–S]. Specifically, we observe that the non-negative

*a*-axis is a phase path for [D–S], with an equilibrium point at .

Next, consider that phase path passing through at *t*=0, and denote the phase path by (*a*,*b*)=(*a*^{i}(*t*),*b*^{i}(*t*)) in *t*≥0. As is a positively invariant region for [D–S], then we can conclude immediately that
3.5
for all *t*≥0. It then follows from (2.26) that,
3.6
An integration of (3.6) then gives
3.7
from which we obtain, with (3.5),
3.8
It also follows from (2.26) and (2.27) that, with *w*=4*a*^{i}+*b*^{i},
3.9
where
3.10
An integration of (3.9) establishes that
3.11
from which it follows via (3.5) that
3.12

and so,
3.13
Thus, we have established *a priori* bounds for (*a*^{i}(*t*),*b*^{i}(*t*)) in (3.8) and (3.13), and so we can conclude that the phase path (*a*^{i}(*t*),*b*^{i}(*t*)) exists globally in *t*≥0, and, moreover, remains bounded in for all *t*≥0. We can immediately conclude, via the Poincaré–Bendixson theorem [18] that the phase path (*a*^{i}(*t*),*b*^{i}(*t*)) must approach one of the following as , namely

an equilibrium point of [D–S] in .

a periodic orbit of [D–S] in and

a homoclinic or heteroclinic orbit of [D–S] in .

Moreover, it follows from (3.7), (3.11) and (i)–(iii) that the *ω*-limit set of the phase path (*a*^{i}(*t*),*b*^{i}(*t*)) must be contained in , and so the phase path (*a*^{i}(*t*),*b*^{i}(*t*)) must

enter the region at some

*t**≥0 and remain in for all*t*≥*t**,

or,

approach the region (having at least one limit point on the boundary of

*R*) as

where *R* is the interior of a quadrilateral, given by,
3.14
and is illustrated in figure 2. It also follows from (i) and (ii), together with (I) and (II), that
3.15
where *E* is the set of equilibrium points of [D–S] in , *P* is the set of periodic orbits of [D–S] in and H is the set of homoclinic and heteroclinic orbits of [D–S] in .

### (b) Equilibrium points and equilibrium point bifurcations

We now identify all equilibrium points of [D–S] in , and their associated bifurcation structure. Let (*a*,*b*)=(*α*,*β*) be an equilibrium point of [D–S] in . It immediately follows from (3.14) and (3.15), that, , and so,
3.16
Substitution into [D–S] then requires that
3.17
and
3.18
It follows immediately that
3.19
is an equilibrium point. All other equilibrium points must satisfy
3.20
and
3.21
with . It follows from (3.21) and (3.20) that *α*>0 and *β*>0. Then, via (3.21)
3.22
By substituting from (3.22) into (3.20) and rearranging, we obtain
3.23
where
3.24
and
3.25
We must now examine the roots of (3.23) for . It is straightforward to establish that (3.23) has no roots when . Thus, when , the [D–S] has just one equilibrium point in , being ** e_{0}** as identified in (3.19). However, for , the situation is different. Specifically, regarding as

*fixed*, with as an

*unfolding parameter*and as a

*bifurcation parameter*, then there is a value , for each , such that (3.23) has no roots for , has a single root (of multiplicity two) for , and has two simple roots for . We conclude that, in addition to the equilibrium point

**given in (3.19), [D–S] has no equilibrium points in when . However, an equilibrium point saddle–node bifurcation occurs when , creating two additional equilibrium points in when . We denote these two additional equilibrium points by, 3.26 with**

*e*_{0}*α*

_{+}>

*α*

_{−}. It follows from §3 that

**,**

*e*_{+}**∈**

*e*_{−}*R*. Moreover, in the unfolding plane, the curve has the parametric representation, 3.27

with . In fact, the parameter *α* in (3.27) represents the *a*-coordinate of the bifurcating equilibrium point. We are now in a position to sketch the equilibrium point bifurcation diagrams for [D–S] on the and the bifurcation planes, and the equilibrium point bifurcation curves in the unfolding plane. Qualitative sketches of the bifurcation diagrams are given in figures 3 and 4 (which are obtained from (3.19) and the locus of (3.23), together with (3.22)), together with the locus of the equilibrium points (with increasing ) in on the (*a*,*b*) phase plane in figure 5. The qualitative structure of the unfolding plane is sketched in figure 6 (via the locus of (3.27)). We see that, for any , , when , then the only equilibrium point of [D–S] is that representing the unreacting state. However, with , then [D–S] has two additional equilibrium points representing fully reacting states. We now discuss the local dynamic stability of the equilibrium points ** e_{0}**,

**and**

*e*_{+}**.**

*e*_{−}### (c) Local stability of equilibrium points

Here, we examine the local stability of the equilibrium points of [D–S] via the linearization theorem. The eigenvalues of the associated linearization of [D–S] at each respective equilibrium point ** e_{0}**,

**and**

*e*_{+}**are the roots of the characteristic equation, 3.28 where 3.29 and 3.30**

*e*_{−}It follows from (3.28)–(3.30), after a little algebra, that,

is a hyperbolic stable node for all and , .

is a hyperbolic saddle point for all and , .

However, an analysis of (3.28)–(3.29) establishes that there exists a value , which depends upon , with, for all , and for which,

is a hyperbolic stable node for all and .

is a hyperbolic unstable node or spiral for all and .

] is a hyperbolic stable node or spiral for and .

Here, for all , and the eigenvalues of the associated linearization of [D–S] at ** e_{−}** when are a purely imaginary pair. It is also straightforward, via (3.23) and (3.28)–(3.30), to obtain the following parametrization for the curve in the unfolding plane, namely
3.31
for parameter . In fact, the parameter

*α*in (3.31) represents the

*a*-coordinate of the equilibrium point

**when and .**

*e*_{−}We next observe that the curves and intersect tangentially at the point,
3.32
in the unfolding plane, with,
3.33
When the eigenvalues of the associated linearization of [D–S] at the coincident equilibrium points ** e_{−}**=

**are both zero. This point in the unfolding plane is the single unfolding point, with fixed. In particular, it has a parametrization with , which is readily obtained from (3.27) and (3.31) as, 3.34 for , with the parameter**

*e*_{+}*α*being the

*a*-coordinate of the coincident equilibrium points

**and**

*e*_{+}**at the unfolding point , The curves (referred to as SN) and (referred to as H) divide the unfolding plane into three disjoint regions**

*e*_{−}*A*,

*B*and

*C*, and this is illustrated in figure 7. We have now characterized the local stability properties of the equilibrium points

**,**

*e*_{0}

*e*_{+}and

**. We can next use this information to locate local bifurcations at the equilibrium points**

*e*_{−}**,**

*e*_{0}**and**

*e*_{+}**.**

*e*_{−}### (d) Local bifurcations

We now examine the possibility of local bifurcations occurring at equilibrium points following the conclusions of §3*c*. No local bifurcations occur at the equilibrium point ** e_{0}**, because this equilibrium point remains a hyperbolic stable node throughout

*A*∪

*B*∪

*C*. The equilibrium points

**and**

*e*_{+}**are created into at a saddle–node bifurcation in the (**

*e*_{−}*a*,

*b*) phase plane at (

*a*,

*b*)=

**=**

*e*_{+}**when . No further local bifurcations occur at**

*e*_{−}**for all because**

*e*_{+}**remains as a hyperbolic saddle point. However, for each , a further local bifurcation occurs at**

*e*_{+}**, when , and this is a Hopf bifurcation creating a unique limit cycle into either of (stable) or (unstable). Further details of the nature of this Hopf bifurcation are postponed until §3**

*e*_{−}*g*.

### (e) Heteroclinic connections

A straightforward centre manifold calculation, when and (*a*,*b*)∼** e_{±}**, reveals that when the equilibrium points

**and**

*e*_{+}**have a unique heteroclinic connection, which we label as , as shown in figure 8. This heteroclinic connection can subsequently be broken only in by either a local bifurcation at**

*e*_{−}**and/or**

*e*_{+}**, or a global bifurcation away from**

*e*_{−}**and**

*e*_{+}**.**

*e*_{−}### (f) Phase portrait at infinity

Here, we examine the structure of the phase portrait of [D–S] in the neighbourhood of the arc at infinity in . A straightforward application of the Poincaré projection establishes that in , the circular arc at infinity forms a phase path, containing exactly two equilibrium points at infinity at, and , where is an unstable node and is a saddle point. The arc at infinity in forms a heteroclinic connection, connecting to . The phase portrait at infinity in is sketched in figure 9. Observe that no phase paths are attracted to the arc at infinity in , in accordance with the conclusions of §3*a*.

### (g) Periodic orbits and bifurcation to periodic orbits

We consider the periodic orbits of [D–S] in , and their associated bifurcation structure. First, we observe that as is a positively invariant region for [D–S], then any periodic orbit of [D–S] which has non-trivial intersection with , must lie wholly within . Thus, no periodic orbit in can surround the equilibrium point ** e_{0}**. However, it follows via index theory, that any periodic orbit which may exist in must surround at least one finite equilibrium point in . We conclude, via §3

*b*, that there are no periodic orbits in for . Similarly, index theory rules out the existence of any periodic orbits in for (the index of

**at SN is zero). Furthermore, any periodic orbit which may exist in for must surround the equilibrium point**

*e*_{±}**, and no other equilibrium point, which again follows from index theory.**

*e*_{−}As identified in §3*b*, a Hopf bifurcation occurs for each , at the equilibrium point, ** e_{−}** when,
3.35
A simple (but lengthy) calculation, supported by numerical integration, determines the criticality, or degeneracy, of this Hopf bifurcation. There are three cases. In particular, there is a value , depending on , with , and for which the Hopf bifurcation at

**, when is degenerate. The three cases are**

*e*_{−}*Case* (*a*): . In this case, the Hopf bifurcation at ** e_{−}**, when , is supercritical, creating a (unique) unstable, limit cycle in for each , bifurcating out of the equilibrium point

**.**

*e*_{−}*Case* (*b*): . In this case, the Hopf bifurcation at ** e_{−}**, when is degenerate, and supercritical, creating a (unique) unstable limit cycle in for each , bifurcating out of the equilibrium point

**.**

*e*_{−}*Case* (*c*): . In this case, the Hopf bifurcation at ** e_{−}**, when , is subcritical, creating a (unique) stable limit cycle in for each , bifurcating out of the equilibrium point

**.**

*e*_{−}In each case, we will refer to this periodic orbit as *P*^{−}, with period *T*^{−} and amplitude (relative to the equilibrium point ** e_{−}**)

*A*

^{−}, which depend upon and (with fixed ). In particular, it is straightforward to establish that 3.36 as , while 3.37 as for all , with, 3.38 Now, numerical evidence and consistency on a circle around the unfolding point in the unfolding plane establishes that for , a homoclinic bifurcation occurs when, , via the formation of a homoclinic orbit on the saddle point equilibrium point

**. This homoclinic bifurcation creates a unique unstable limit cycle in for , and unique stable limit cycle in for , with the homoclinic bifurcation changing from subcritical for , to supercritical for . The limit cycle is created at finite amplitude from the homoclinic orbit on**

*e*_{+}**at . Numerical evidence establishes that (with depending on ). On the unfolding plane, we now have the situation which is illustrated in figure 10. The change in criticality of the Hopf bifurcation at induces a periodic saddle–node bifurcation for each , when, , with for all . The periodic saddle–node bifurcation is absorbed at the homoclinic bifurcation at its change in criticality, when . The periodic saddle–node bifurcation gives rise to both an unstable (outer) and stable (inner) limit cycle surrounding the equilibrium point**

*e*_{+}**, being created at finite amplitude as increases through , with a single bistable limit cycle appearing at , from which they are created in . The inner limit cycle is that which we have labelled as**

*e*_{−}*P*

^{−}. We now label the outer limit cycle as

*P*

^{+}.

We remark that the limit cycle created by the homoclinic bifurcation is created at finite amplitude and infinite period, whereas the two limit cycles created by the periodic saddle–node bifurcation are created at finite amplitude and finite period.

Numerical evidence then confirms that, for each , the full unfolding plane is as illustrated in figure 11, where the plane is now subdivided into the disjoint regions *A*′–*F*′. We are now in a position to use figure 11 to sketch bifurcation diagrams for periodic orbits at each , with the amplitude *A* represented as on the periodic orbit. These are illustrated in figure 12.

In figure 12, represents the amplitude of the periodic orbit at the homoclinic bifurcation and *A*_{p} represents the amplitude of the periodic orbit at the periodic saddle–node bifurcation. It should also be noted that in the (*a*,*b*) phase plane, each periodic orbit, *P*^{−} and *P*^{+}, must lie in the region , as a consequence of §3*a*, and in particular, (3.15).

The equilibrium point bifurcation diagram in figure 3 can now be combined with the bifurcation diagrams for periodic orbits in figure 12 to construct the full bifurcation diagrams illustrated in figure 13. Again, we remark that, in the (*a*,*b*) phase plane, all equilibrium points, periodic orbits and homoclinic orbits identified in the bifurcation diagrams of figure 13 lie in the region , as a consequence of §3*a* and (3.15).

### (h) Attractors for phase paths in

We can now determine the attracting set for each phase path in . With , in what follows, we will denote the phase path passing through (*a*_{0},*b*_{0}) at *t*=0, when *t*≥0, as (via §3*a*). We denote the *ω*-limit set of *P*_{0} by . In addition, at the relevant parameter value in the unfolding plane, we denote the stable limit cycle as , the unstable limit cycle as , and the bistable limit cycle as . We denote the homoclinic orbit on the saddle point ** e_{+}** as , and the stable manifold of the saddle–point

**as . Similarly, we write, , , , . We now examine**

*e*_{+}*ω*

_{0}for each

*P*

_{0}with , for each of the regions

*A*′–

*F*′ on the unfolding plane. In addition, we introduce the notation ∂

*A*′

*B*′ to be the boundary between regions

*A*′ and

*B*′ in the unfolding plane, with similar notation for all other boundaries in the unfolding plane. We have 3.39 3.40 3.41 3.42 3.43 3.44 3.45 3.46 3.47 with being the interior of , and being the interior of . 3.48 3.49 3.50 3.51 with being the interior of . 3.52 3.53 3.54 3.55 3.56 3.57 3.58 3.59 3.60 with being the stable manifold of the saddle–node

*ε*

_{±}. 3.61 3.62 3.63 3.64 3.65 with being the interior of the homoclinic orbit on

**. 3.66 3.67 3.68 3.69 3.70 3.71 3.72 3.73**

*e*_{+}All cases are now complete. We can now immediately identify from the above cases the stable attractors in , and this is exhibited in table 1. It is appropriate at this point to interpret the nature of the stable attractors as power generating states for the SOFC. The equilibrium point *ε*_{0} represents a dormant, non-power generating steady state, whereas the equilibrium point *ε*_{−} represents a constant, non-zero, power generating steady state. The stable periodic orbit represents a time-periodic, non-zero, power generating state.

### (i) Global phase portraits

Here, we sketch the global phase portraits for the regions *A*′−*F*′ in the unfolding plane via the information established in §3. The phase portraits on the boundaries between the regions *A*′–*F*′ may also be readily constructed, however, they are omitted for brevity. The qualitative sketches of the global phase portraits are given in figure 14, where stable limit cycles are shown with a single dash, while unstable ones have a double dash. In addition, the stable manifold of the equilibrium point ** e_{+}** is shown with a triple dash.

## 4. Numerical integration

Here, we exhibit sample numerical integrations to [D–S]. We take throughout, and select from illustrative regions in the unfolding plane. [D–S] is integrated forward into *t*>0, from initial conditions , using the MATLAB *ode*23*s* solver. The graphs of *a*(*t*) and *b*(*t*) against *t*≥0 are shown in figure 15. The initial conditions are chosen to illustrate the stable attracting sets in table 1 and are shown in table 2.

In figure 15*a,* all phase paths approach the stable equilibrium point as , in accordance with the phase portrait sketched in figure 14*a* for the region *A*′ in the unfolding plane.

Figure 15*b* shows the phase paths which approach the stable equilibrium points ** e_{0}** and

**=(**

*e*_{−}*α*

_{−},

*β*

_{−}), representative of region

*F*′ in the unfolding plane, whose phase portrait was sketched in figure 14

*b*. Phase paths starting within the interior of spiral onto the equilibrium point

**, whereas those outside approach the equilibrium point**

*e*_{−}**.**

*e*_{0}In figure 15*c*, which represents the region *E*′ of the unfolding plane, with phase portrait sketched in figure 14*c*, we see the same attracting set as in figure 15*b*. Unique to the region *E*′, however, the set of initial conditions which approach the equilibrium point ** e_{−}** is bounded by the unstable limit cycle , as opposed to being bounded by the stable manifold , of the equilibrium point

**for region**

*e*_{+}*F*′. This is demonstrated in figure 15

*c*by taking the initial conditions for phase path 1 just inside the unstable limit cycle, and the initial conditions for phase path 2 just outside the unstable limit cycle (table 2).

Figure 15*d* which represents region *B*′ in the unfolding plane, with phase portrait sketched in figure 14*d*, has phase paths which approach the stable equilibrium point ** e_{0}**, although we see the appearance now of the unstable equilibrium point

**and the saddle point**

*e*_{−}**. Phase paths can be seen spiralling away from the equilibrium point**

*e*_{+}**towards the stable equilibrium point**

*e*_{−}**.**

*e*_{0}Figure 15*e*–*g* shows various phase paths for region *C*′ in the unfolding plane, whereas the phase paths shown in figure 15*h*–*j* represent region *D*′. In figure 15*e*, showing phase path 1 in region *C*′, and (15*h*), showing phase path 1 in region *D*′, starting close to the equilibrium point ** e_{−}** the phase paths can be seen winding onto the stable limit cycle . Figure 15

*f*,

*g*, shows phase paths 2 and 3 in region

*C*′, can be seen approaching different attractors in the phase plane, namely the stable limit cycle in figure 15

*f*and the equilibrium point

**in figure 15**

*e*_{−}*g*. The initial conditions are taken to be just inside, and just outside of the unstable limit cycle . Phase path 2 shown in figure 15

*i*, representing region

*D*′ in the unfolding plane, approaches the stable limit cycle but starts much further out than phase path 2 in figure 15

*f*. The set of points which approach in region

*D*′ are bounded only by the stable manifold of the equilibrium point

**, and not by the unstable limit cycle , as in region**

*e*_{+}*C*′. Hence, initial conditions may be taken much further from the equilibrium point

**and still approach . Figure 15**

*e*_{−}*j*shows phase path 3 from region

*D*′ in the unfolding plane, which starts outside of the stable manifold of the equilibrium point

**and approaches the equilibrium point**

*e*_{+}**.**

*e*_{0}## 5. Discussion

### (a) Model summary

A model has been derived from fundamental chemical kinetics and transfer processes, for an internal reforming SOFC which uses methane as the primary fuel. The system was considered to be governed by the reaction scheme set out in (2.1)–(2.4), in accordance with [6]. The dynamical system (2.12)–(2.14), which describes the temporal dynamics of the chemical species within the anode, was derived based on this reaction scheme, in addition to the transfer processes within the cell. Typical scales were chosen and the dynamical system (2.12)–(2.14) was non-dimensionalized with respect to these scales, leading to a reduced planar dynamical system [D–S], represented by equations (2.26)–(2.27). This system describes the temporal dynamics of the concentrations (*a*(*t*),*b*(*t*)), which are the non-dimensional concentrations of methane and hydrogen within the SOFC, respectively.

### (b) Analysis

Having obtained the rationally reduced planar dynamical system [D–S], it has been established that the quadrant is a positively invariant region, and the existence of attracting sets for [D–S] within this region have been determined, with the *ω*-limit set for phase paths in being given by (3.15). Following this, the equilibrium points of [D–S] were identified and were analysed for stability. It was found that there exists a stable equilibrium point *e*_{0}, irrespective of the parameter regime, but that two more equilibrium points emerge via a saddle–node bifurcation as the parameter passes through the point . One of these (*e*_{+}) is a saddle point, whereas the other (*e*_{−}) is either stable or unstable depending on the parameter regime. The equilibrium point (*e*_{−}) changes stability via a Hopf bifurcation as the parameter passes through the point . The nature of the limit cycles associated with the Hopf bifurcations, which appear in the phase plane, are determined by the regions of the unfolding plane. Up to two limit cycles may coexist in the phase plane, surrounding the equilibrium point *e*_{−}, however, at most, one of these limit cycles is stable. Sketches of the phase portraits in the various regions of the unfolding plane may be seen in figure 14. Sample numerical integrations were given in §4 to illustrate the stable attracting sets of [D–S].

### (c) Conclusion

We conclude from the model presented in this paper, that the temporal dynamics of an internal reforming SOFC fuelled with methane are capable of displaying a wide range of behaviours. These include various bifurcations, multiple steady states and oscillatory phenomena, which is consistent with the literature [3,4,5]. However, the approach here is unique in identifying the key parameters which determine the bifurcations, and onset of nonlinear oscillations, as the diffusivities of the chemical species to the reaction zone, and the composition of the fuel stream. Practical operating conditions for an SOFC may be directly linked to the key parameters in the model, and thus oscillations may be induced by controlling these parameters in such a way that they lie within the correct regions of the unfolding plane.

### (d) Further work

As the data on oscillations in SOFCs which is available in the literature, are not specific to an internal reforming SOFC fuelled with methane, experiments have been planned to evaluate the model based on the findings in this paper. The results of these experiments will be published in a subsequent paper. In addition to this, the model may be extended in a number of ways in order to make it more general, beginning with taking the parameter to be of O(1). We may then include the temporal dynamics of the fuel stream in the model, thermal effects and also account for the decomposition/direct oxidation of methane. The model presented here indicates that it is of great importance to have accurate values for both diffusivities of chemical species, as well as the reaction rate constants under the conditions which are to be tested.

## Funding statement.

J.D.S. acknowledges financial support from the EPSRC supported Doctoral Training Centre in Hydrogen, Fuel Cells and their Applications at the University of Birmingham.

## Appendix A. Primary symbols and notations

symbols | description |
---|---|

a | concentration of CH_{4} |

b | concentration of H_{2} |

c | concentration of H_{2}O |

a_{s} | scaled concentration of CH_{4} |

b_{s} | scaled concentration of H_{2} |

c_{s} | scaled concentration of H_{2}O |

T | scaled time constant |

ε | ratio of reaction rate constants (dimensionless) |

dimensionless diffusivity of CH_{4} | |

dimensionless diffusivity of H_{2} | |

initial concentration of CH_{4} (dimensionless) | |

initial concentration of H_{2}O (dimensionless) | |

x | pool concentration of O_{2} |

(e_{0},e_{+},e_{−}) | equilibrium points |

quadrant of the phase plane defined by Q={(a,b)=a>0 and b>0} | |

stable limit cycle | |

unstable limit cycle |

- Received August 15, 2013.
- Accepted December 19, 2013.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.