Study of Static and Dynamic Behavior of a Membrane Reactor for Hydrogen Production

Abstract

The paper investigates the stability and bifurcation phenomena that can occur in membrane reactors for the production of hydrogen by ammonia decomposition. A simplified mixed model of the membrane reactor is studied and two expressions of hydrogen permeation are investigated. The effect of the model design and operating parameters on the existence of steady state multiplicity is discussed. In this regard, it is shown that the adsorption-inhibition effect caused by the competitive adsorption of ammonia can lead to the occurrence of multiple steady states in the model. The steady state multiplicity exists for a wide range of feed ammonia concentration and reactor residence time. The effect of the adsorption constant, the membrane surface area and its permeability on the steady state multiplicity is delineated. The analysis also shows that no Hopf bifurcation can occur in the studied model.

1. Introduction

Research in the use of renewable energy is being actively pursued in order to mitigate global climate challenges and preserve sustainability. In this regard, hydrogen is a promising clean energy that can be utilized in many applications in the chemical industry [1] and in the production of such systems as membrane fuel cells that can be used in electronic devices and transportation [2,3]. However, the use of hydrogen either in gaseous or liquid form is known to be hampered by difficulties in storage and transportation [4]. In order to overcome these logistic problems, active research is being carried out in the use of other energy carriers that can provide on site hydrogen production. Methanol and ammonia are two candidates for this purpose due to their high energy densities and limited storage problems. However, methanol reforming is known to produce greenhouse gases in the form of carbon oxides. Ammonia, on the other hand, is carbon free and can produce via catalytic or thermal cracking only hydrogen and nitrogen. Moreover, ammonia has a well-established storage and distribution infrastructure making it an economically feasible large scale energy carrier [4,5,6].

There has been active research in recent years in the use of catalytic membrane reactors for cracking ammonia to produce pure hydrogen [7,8,9,10,11,12,13,14,15]. de Nooijer et al. [10], for instance, studied the thermal stability and performance of Pd-Ag films used in membrane reactors for hydrogen production via the reforming of different feedstocks. Petriev et al. [11] carried out experimental modification of the surface of Pd-Ag films and showed that hydrogen permeability rate depends not only on the surface area but also on the structure of the membrane modifying layer. Lytkina et al. [12] carried out a comparative study of methanol steam reforming in conventional and membrane reactors. It was found that the use of a membrane reactor allowed the production of ultra pure hydrogen (i.e., CO-free) that can be directly fed to a polymer electrolyte membrane fuel cell. Barba et al. [13] investigated the merits of a Reformer and Membrane Module configuration which allowed to alternate between steam reforming (SR) reactions and hydrogen separation in the membrane module in order to overcome the thermodynamic limitations of the conventional (SR). Cechetto et al. [14] carried out an experimental decomposition of ammonia in a Pd-based membrane reactor over a Ru-based catalyst and compared its results with those of a conventional packed bed reactor. Membrane reactors were found to enjoy clear advantages over traditional reactors since the reaction and separation tasks are carried out in one single unit which results in a compact device that is more economical over conventional reactors. The use of membrane reactors offers an additional advantage in the case of ammonia decomposition to hydrogen. Compared to steam reforming of hydrocarbons, the ammonia decomposition achieve almost complete equilibrium conversion at temperatures close to 400 °C. However, at low temperatures the rate of hydrogen production is inhibited by hydrogen itself [14]. Using a membrane reactor could overcome both thermodynamic and kinetic limitations because the removal of hydrogen as it is produced would shift the equilibrium towards higher conversions. Operating at lower temperatures represents also an advantage in terms of energy requirements of the process.

While remarkable achievements have been accomplished to improve catalysts for ammonia decomposition, theoretical studies on using ammonia as feedstock for the production of hydrogen are limited [15,16]. Murmura et al. [17] carried out recently a numerical analysis of a model of membrane reactor for the production of pure hydrogen from ammonia. The authors concluded that a membrane reactor allowed the operation at low temperature, the selective removal of hydrogen and the enhancement of the reaction rate.

One area that deserves theoretical investigation is the stability behavior of such reactors. Unlike stirred tank reactors [18] or bioreactors [19], the theoretical analysis of stability of catalytic membrane reactors in general and the production of hydrogen in particular did not receive much attention in the literature. Reactors, including membrane ones, can present operational problems that can manifest themselves in form of steady state input, output multiplicities and undesired oscillations. Input multiplicities occur when different values of an input variable produce the same value of the output variable. Output multiplicities occur when the same value of an input variable produces different values of the output variable. The hysteresis phenomenon is the most common form of output multiplicity and is associated with the existence of a region of unstable behavior. Besides steady state multiplicity (input and output multiplicities), a nonlinear model can exhibit periodic or non-periodic oscillations [18,19]. Solutions that form a closed curve for some model parameters’ values are called limit cycle (or periodic orbit). This is associated with bifurcation from equilibria to limit cycles. This bifurcation is commonly called a Hopf bifurcation [18]. All such nonlinear phenomena (steady state multiplicity or oscillations) are generally detrimental to the operation of the reactor. An early detection of such operating regions in reactors would be useful since this would allow the removal or at least the reduction of these operational problems in the early stage of process design and would allow the selection of operating parameters that avoid such behavior. Recently, Murmura et al. [20] investigated the bistability in membrane reactors for the ammonia cracking to hydrogen. The authors studied both mixed reactor model and fixed bed reactor model and concluded that bistability is possible in such models.

Inspired by these results, the aim of the present work is to study in more detail the stability of a model of membrane reactor for the production of hydrogen from ammonia for two different expressions of hydrogen permeation rate, and to examine the effect of operating and design parameters on the existence of such behavior. The rest of the paper is structured as follows: We present the model in Section 2. In Section 3 we analyze the boundedness of the model solutions. In Section 4 we investigate the stability of the model and in Section 5 we present numerical simulations for the different expressions of hydrogen permeation rate.

2. Process Model

Ammonia NH${}_3$ catalytic decomposition is assumed to follow the reaction

$$\begin{array}{c}\hfill 2N{H}_3⟺3H_2+N_2.\end{array}$$

(1)

We consider an isothermal membrane reactor with pure ammonia feed. The membrane reactor is represented by the following set of ordinary differential equations resulting from mass balances on ammonia (denoted A) and hydrogen (denoted H).

$$\begin{array}{ccc}\hfill \frac{d{C}_A}{dt}& =& -2r+\frac{{\mathit{C}}_{\mathit{A}0}-{\mathit{C}}_{\mathit{A}}}{\theta },\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{ d{C}_H}{dt}& =& -\frac{{\mathit{C}}_{\mathit{H}}}{\theta }+3r-\frac{J{\mathit{a}}_{\mathit{m}}}{F\theta },\hfill \end{array}$$

(3)

t is time, $C_A$ and $C_H$ are the concentrations of ammonia and hydrogen respectively, F the volumetric flow rate, $\theta$ the reactor residence time (volume over volumetric flow rate), r the rate of reaction, $a_m$ the membrane area and J the hydrogen flux through the selective membrane. For the reaction rate $\left(r\right)$, studies have shown that at high temperatures, above 500 °C, the reaction rate depends only on the partial pressure of ammonia [21], whereas at temperatures lower than 500 °C and high hydrogen partial pressures, the reaction rate is inhibited by hydrogen and can be described by the following relation [21,22]:

$$\begin{array}{c}\hfill r=k{\left(\frac{{{\mathit{P}}_{\mathit{A}}}^2}{{{\mathit{P}}_{\mathit{H}}}^3}\right)}^{0.25}.\end{array}$$

(4)

where $P_A$ and $P_H$ are partial pressures of ammonia and hydrogen respectively and k is the reaction rate constant.

As to hydrogen permeation, assumed to be through dense palladium membranes, it is commonly described by Sieverts’ law [23]:

$$\begin{array}{c}\hfill J=\varphi (\sqrt{p_H^r}-\sqrt{p_H^p}),\end{array}$$

(5)

where $p_H^r$ and $p_H^p$ are partial pressures of hydrogen in the retentate and permeate sides, respectively, and $\varphi$ is the permeance of the membrane whose value depends on the membrane characteristics and on the operating temperature. Sieverts Law is widely applied in analyzing the steady state hydrogen permeation through Pd-based membranes. Alraeesi and Gardner [24], for instance, have assessed the validity of this law and modelled the effects of pressure and temperature on the accuracy of this power law.

Neglecting hydrogen partial pressure in the permeate side, and assuming an ideal gas law, J becomes:

$$\begin{array}{c}\hfill J=\varphi \sqrt{RT}\sqrt{C_H},\end{array}$$

(6)

where T is the temperature and R is the ideal gas constant. Hydrogen permeation is also known to be inhibited by the adsorption of ammonia on the membrane. A number of expressions were proposed in literature [25,26] to account for this inhibition. We will consider in this paper the following expression

$$\begin{array}{c}\hfill J=\varphi \frac{\sqrt{RT}\sqrt{C_H}}{1+K_{ad}{C}_A},\end{array}$$

(7)

where $K_{ad}$ is the adsorption constant of ammonia.

3. Model Analysis

Before we analyze the stability of the model, we first need to see whether the proposed model is mathematically correct and that it can represent real situations. This means that we need to check whether the solutions $(C_A,C_H)$ of the model remain bounded. We consider in the analysis that follows a general dependence of rate $\left(r\right)$ and flux $\left(J\right)$ on $C_A$ and $C_H$. To facilitate our analysis we reformulate the model by defining

$$\begin{array}{c}\hfill (C_A,C_H):=(x_1,x_2)=x\in \Omega ,\end{array}$$

(8)

and

$$\begin{array}{c}\hfill f(x_1,x_2)=(f_1(x_1,x_2),f_2(x_1,x_2)),\end{array}$$

(9)

where

$$\begin{array}{c}\hfill f_1(x_1,x_2)=-2r(x_1,x_2)-\frac{(x_1-C_{A0})}{\theta },\end{array}$$

(10)

and

$$\begin{array}{c}\hfill f_2(x_1,x_2)=-\frac{a_m}{F\theta }J(x_1,x_2)+3r(x_1,x_2)-\frac{x_2}{\theta }.\end{array}$$

(11)

$\Omega$ is the domain of interest for the model and is defined as follow:

$$\begin{array}{c}\hfill \Omega :=\{x\in (0,\infty )\times (0,\infty \left)\right\}.\end{array}$$

(12)

Then the system of differential equations can be written as

$$\dot{x}=f\left(x\right),x\left(0\right)=x_0>0.$$

(13)

Theorem 1.

The nonlinear autonomous model (13) has a bounded unique solution over Ω if $x\left(0\right)\in \Omega$.

Proof. 

Here both functions $f_1$ and $f_2$ are continuous and Lipchitz over $\Omega$ which guarantee the uniqueness of solutions. f is also a differentiable function over $\Omega$, and we have:

$$\begin{array}{c}\hfill |f_1|\le c_1+c_2\parallel x\parallel \le c_1+c_2\parallel x\parallel ,\end{array}$$

(14)

and

$$\begin{array}{c}\hfill |f_2|\le c_3|x_2|\le c_3\parallel x\parallel .\end{array}$$

(15)

Therefore

$$\begin{array}{c}\hfill {\parallel f\parallel }_{\infty }\le c_1+c\parallel x\parallel ,\end{array}$$

(16)

where $c_i$, $(i=1,3)$ depend on the parameters and $c=max\{c_2,c_3\}$. The above inequality shows a linear growth of the nonlinearity f which guarantees the uniqueness of solutions of (13) which is bounded [27,28]. □

4. Local Stability

In this section the stability of the steady state solutions are studied. We consider a general dependence of flux J on $C_A$ and $C_H$ but we consider the rate equation as defined by Equation (4), and written as function of concentrations,

$$\begin{array}{c}\hfill r(C_A,C_H)=k{\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{{\mathit{k}}_1{\mathit{C}}_{\mathit{H}}}^3}\right)}^{0.25}\mathrm{with}{k}_1=RT.\end{array}$$

(17)

The Jacobian matrix of the system is

$$\begin{array}{cc}\hfill \mathrm{J}& =\left[\begin{array}{cc}{j}_{11}& j_{12}\\ j_{21}& j_{22}\end{array}\right],\hfill \end{array}$$

(18)

where,

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}{j}_{11}=& -{\theta }^{-1}-0.50\frac{2k{\mathit{C}}_{\mathit{A}}}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}{\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}\right)}^{-0.75},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill j_{12}=& 0.75\frac{2k{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^4}{\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}\right)}^{-0.75},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} j_{21}=& 0.50\frac{3k{\mathit{C}}_{\mathit{A}}}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}{\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}\right)}^{-0.75}-\frac{\left(\frac{\partial }{\partial {\mathit{C}}_{\mathit{A}}}J\left({\mathit{C}}_{\mathit{A}},{\mathit{C}}_{\mathit{H}}\right)\right){\mathit{a}}_{\mathit{m}}}{F\theta },\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} j_{22}=& -{\theta }^{-1}-0.75\frac{3k{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^4}{\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}\right)}^{-0.75}-\frac{\left(\frac{\partial }{\partial {\mathit{C}}_{\mathit{H}}}J\left({\mathit{C}}_{\mathit{A}},{\mathit{C}}_{\mathit{H}}\right)\right){\mathit{a}}_{\mathit{m}}}{F\theta }.\hfill \end{array}$$

(22)

The trace and determinant of the matrix are:

$$\begin{array}{cc}\hfill Trace\left(J\right)=& -\left(T_1+T_2\right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}{T}_1& =\frac{2+\frac{\left(\frac{\partial }{\partial {\mathit{C}}_{\mathit{H}}}J\left({\mathit{C}}_{\mathit{A}},{\mathit{C}}_{\mathit{H}}\right)\right){\mathit{a}}_{\mathit{m}}}{F}}{\theta },\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{T}_2& =\frac{\left(0.5\times 2{\mathit{C}}_{\mathit{H}}+0.75\times 3{\mathit{C}}_{\mathit{A}}\right)k{\mathit{C}}_{\mathit{A}}}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^4}{\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}\right)}^{-3/4},\hfill \end{array}$$

(25)

and

$$\begin{array}{cc}\hfill Det\left(J\right)=& \left[{\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}\right)}^{-3/4}\right]A,\hfill \\ \hfill A=& 0.75\frac{2k{{\mathit{C}}_{\mathit{A}}}^2{\mathit{a}}_{\mathit{m}}\frac{\partial }{\partial {\mathit{C}}_{\mathit{A}}}J\left({\mathit{C}}_{\mathit{A}},{\mathit{C}}_{\mathit{H}}\right)}{\theta {\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^{4}F}\hfill \\ \hfill & +0.5\frac{\left(\frac{\partial }{\partial {\mathit{C}}_{\mathit{H}}}J\left({\mathit{C}}_{\mathit{A}},{\mathit{C}}_{\mathit{H}}\right)\right){\mathit{a}}_{\mathit{m}}}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3{\theta }^{2}F}\left(2.0{\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}\right)}^{3/4}{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3+1.0\times 2k{\mathit{C}}_{\mathit{A}}\theta \right)\hfill \\ \hfill & +\frac{1}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^4{\theta }^2}\left({\left(\frac{{{\mathit{C}}_{\mathit{A}}}^2}{{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^3}\right)}^{3/4}{\mathit{k}}_1{{\mathit{C}}_{\mathit{H}}}^4+0.75\times 3k{{\mathit{C}}_{\mathit{A}}}^2\theta +0.5\times 2k{\mathit{C}}_{\mathit{A}}\theta {\mathit{C}}_{\mathit{H}}\right).\hfill \end{array}$$

(26)

Case 1

For first case when J is given by

$$\begin{array}{cc}\hfill J=& \varphi \sqrt{\mathit{RT}}\sqrt{{\mathit{C}}_{\mathit{H}}},\hfill \end{array}$$

(27)

we have,

$$\begin{array}{cc}\hfill \frac{dJ}{d{C}_H}=& \frac{\varphi \sqrt{\mathit{RT}}}{2\sqrt{{\mathit{C}}_{\mathit{H}}}},\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \frac{dJ}{d{C}_A}=& 0.\hfill \end{array}$$

(29)

In this case the trace is negative and the determinant is positive. Thus the steady state solution is always stable.

Case 2

For the second case when J is given by

$$\begin{array}{cc}\hfill J=& \frac{\varphi \sqrt{\mathit{RT}}\sqrt{{\mathit{C}}_{\mathit{H}}}}{1+{\mathit{K}}_{\mathit{ad}}{\mathit{C}}_{\mathit{A}}},\hfill \end{array}$$

(30)

we have,

$$\begin{array}{cc}\hfill \frac{dJ}{d{C}_H}=& \frac{\varphi \sqrt{\mathit{RT}}}{2\sqrt{{\mathit{C}}_{\mathit{H}}}\left({\mathit{K}}_{\mathit{ad}}{\mathit{C}}_{\mathit{A}}+1\right)},\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \frac{dJ}{d{C}_A}=& -\frac{\varphi \sqrt{\mathit{RT}}\sqrt{{\mathit{C}}_{\mathit{H}}}{\mathit{K}}_{\mathit{ad}}}{{\left({\mathit{K}}_{\mathit{ad}}{\mathit{C}}_{\mathit{A}}+1\right)}^2}.\hfill \end{array}$$

(32)

In this case the trace is always negative but the determinant can be either positive or negative. Thus the equilibria can be either stable or unstable. However, a Hopf point can not exist in the model for both expressions of J since the trace is always strictly negative. We conclude therefore that adsorption inhibition is necessary for the occurrence of static multiplicity but oscillatory behavior is ruled out for the chosen expressions of permeation.

5. Numerical Simulations

The nominal values of the model parameters are shown in

Table 1. Nominal Values of Model parameters.

Parameter	Value	Source
C${}_{A0}$	89.36 (mol/m${}^3$)
k	50 (kgH${}_2$·Pa/m${}^3$·s)	[19]
P	5 $\times {10}^5$ (Pa)	[8,29]
T	400 °C	[8,29]
$\varphi$	2.5 × 10${}^{-5}$ (kgH${}_2$/s·m${}^2$·Pa${}^{0.5}$)	[19]
a${}_m$	25 × 10${}^{-4}$ (m${}^2$)	[8]
$\theta$	1.46 (s)	[8]

. There are two design parameters for the model: The permeance ($\varphi$) whose value depends on the membrane characteristics and on the operating temperature, and the membrane area $a_m$. The values of the design parameters were taken from experimentally validated studies [8,19]. As to the operating parameters, the pressure was selected to be $5\times {10}^5$ Pa in accordance with the work in [8,29].

For the operating temperature of 400 °C, it was selected in accordance with considerations mentioned earlier in model description, that operating at low temperatures lower than 500 °C will shift the equilibrium through the selective removal of hydrogen as well as favoring the reaction kinetics by removing the inhibiting component.

The nominal value of feed concentration shown in Table 1 was calculated from the given pressure and temperature conditions i.e., $C_{A0}=\frac{P}{RT}$. Since in practice the ammonia decomposition reaction is generally carried out at pressures between 10${}^5$ Pa and 5 × 10${}^5$ Pa and in the 300–500 °C temperature range, we considered in the simulations a range of feed concentration from 15 mol/m${}^3$ to 105 mol/m${}^3$. The simulations and the plots were extended sometimes above this range in order to show the range outside the multiplicity region. The nominal value of residence time (i.e., ratio of volume to volumetric flow rate) was taken to be 1.46 s which corresponded to a flow rate of 5.0 × 10${}^{-6}$ m${}^3$/s given that the volume of the reactor adopted from [8] is 7.3 × 10${}^{-6}$ m${}^3$. The residence time was changed by varying the feed flow rate in line with the experimental values in [8]. The maximum value of residence time was 42.9 s and corresponded to a small value of flow rate 1.7 × 10${}^{-7}$ m${}^3$/s. The smallest value of residence time was, on the other hand, 0.55 s and corresponded to a flow rate of 1.33 × 10${}^{-5}$ m${}^3$/s. While large values of residence times are not practical from a productivity point of view, the graphs nevertheless were extended to make sure that the model does not predict any other nonlinear phenomena at these values.

The numerical analysis was carried out using elementary continuation techniques with the help of the software AUTO [30]. The AUTO continuation package can trace out the entire steady state branches, locate the static limit points and continue these points in two parameters. Compared to dynamic simulation (i.e., time traces), continuation methods have the advantage of locating both stable and unstable solutions.

Figure 1 shows the continuity diagram for the first expression of the permeation rate (Equation (6)). As it can be seen the diagram features a single stable branch as the concentration of ammonia decreases with the residence time, while the concentration of hydrogen exhibits a maximum of C${}_H$ = 41.96 mol/m${}^3$ at $\theta =13.3$ s before decreasing. Figure 2 shows the continuity diagram when the feed concentration of ammonia is chosen as the bifurcation parameter. Again one stable branch characterizes the behavior of the model. As expected, the concentration of ammonia increases with its feed concentration while the same can be said about the hydrogen concentration which increases almost linearly. Figure 3 shows the bifurcation diagram for the second expression of the permeation rate (Equation (7)). A multiplicity in the form of hysteresis can be seen in the digram. The region of instability extends from $\theta =155.2$ s to $\theta =187.9$ s. The hydrogen concentration can be seen to go through a maximum of C${}_H$ = 117.42 mol/m${}^3$ at $\theta =98.10$ s that occurs outside the hysteresis region. Within the multiplicity region, abrupt changes in the feed conditions can push the operation of the membrane reactor from high conversion to low conversion conditions. The same multiplicity can be seen to occur when the ammonia concentration C${}_{A0}$ is varied (Figure 4). As C${}_{A0}$ increases, both ammonia and hydrogen concentrations increase but a region of hysteresis can be seen to exist between C${}_{A0}$ = 89.4 mol/m${}^3$ (P = 5 × 10${}^5$ Pa) and C${}_{A0}$ = 108.1 mol/m${}^3$ (P = 6.05 × 10${}^5$ Pa). This region can limit the performance of the reactor as the safe operation of the reactor dictates (barring the use of any feedback control) the operation away from the hysteresis region, that is to operate at the lower residence time (i.e., larger flow rates) since operating at larger residence times (lower flow rates) will decrease the productivity.

Next we examine the effect of some operating and design parameters of the membrane reactor on the existence of such multiplicity. Figure 5 shows the loci of the static limit points in the parameter space $(\theta ,C_{A0})$. Each curve represents the locus of the static limit point of Figure 3 and Figure 4. It can be seen that steady state multiplicity does not exist when $C_{A0}$ and $\theta$ are below the cusp point (A, $\theta =50.1$ s and C${}_{A0}$ = 31.2 mol/m${}^3$). However, if either $\theta$ or C${}_{A0}$ is increased, the range of steady state multiplicity increases. The effect of the adsorption constant on the multiplicity is shown in Figure 6a,b. It can be seen that for the nominal value of feed concentration $C_{A0}$, the instability exists only when the adsorption constant is below the cusp (k${}_{ad}$ = 192.2 m${}^3$/mol) that occurs at $\theta$ = 83.7 s. This means that, notwithstanding issues affecting productivity, the region of multiplicity can be avoided by operating the reactor at residence times lower than the cusp point. The range of multiplicity (in term of residence time) increases as the value of $k_{ad}$ decreases. On the other hand, for a given residence time, the region of multiplicity can be seen to exist for any value of feed concentration and can be seen to increase (in term of feed concentration) with the increase in the adsorption constant. Figure 6a,b showed therefore a conflicting effect of the adsorption constant. The selection of values of operating parameters (residence time $\theta$) and feed concentration (C${}_{A0}$) is to be made judiciously in order avoid or at least minimize the instability introduced by the adsorption constant.

The effect of the surface area of the membrane is shown in Figure 7a,b. The multiplicity in term of residence time occurs only when the surface area is larger than a critical value, a${}_m$ = 0.00167 m${}^2$, and the region increases as the surface area is increased. In term of C${}_{A0}$, Figure 6b shows that the multiplicity also exists only for values of $a_m$ larger than certain value, and increases in range as the surface area is increased. Finally, Figure 8 shows the effect of membrane permeance. The diagram is similar to Figure 7 and similar conclusions can be drawn concerning the range of hysteresis in term of $\theta$ and C${}_{A0}$. An important question to ask is whether the range of parameters in which the multiplicity exists can be found in practice. We have seen that when the feed concentration was the bifurcation parameter, the multiplicity region seemed to occur within practical values such as the ones reported in [8]. The values, on the other hand, of residence time at which instability occurred are larger than what would be met in practical operations. However, a note should be made about the ammonia adsorption constant. It is the only parameter for which reliable data are still not available [20], therefore the diagrams could change as more accurate values are obtained. Figure 6a has effectively shown the important effect of the adsorption constant. In any case, experimental investigation of these instability phenomena is still needed to confirm our results as well as those reported in similar previous studies [20,31].

6. Conclusions

The paper investigated the occurrence of steady state multiplicity in a membrane reactor for the cracking of ammonia into hydrogen. Using a simple mixed model for the reactor the analysis was carried out for two expressions of hydrogen permeation rates. An adsorption-inhibition form was found to be necessary for the occurrence of static multiplicity but no oscillatory behavior i.e., Hopf bifurcations was found in the studied expressions of hydrogen flux.

Instability in the form of hysteresis was found to exist for a wide range of operating conditions (residence time and feed ammonia concentration). From a practical point of view, the region of hysteresis is detrimental to the operation of the reactor and should be avoided. In this regard, for given residence time (feed flow rate) lowering feed ammonia concentration decreases the range of instability. On the other hand, for given feed ammonia concentration, operating at smaller residence time (high flow rate) decreases this multiplicity region. Larger values of membrane areas or its permeance were found to accentuate the region of multiplicity.

For given value of feed ammonia concentration, static multiplicity is expected only when the value of adsorption constant is smaller than a critical value. On the contrary, for given flow rate the range of static multiplicity, in term of feed ammonia concentration is expected to accentuate when the value of adsorption constant is large. A careful choice should be made for the residence time (i.e., flow rate) and feed concentration (i.e., pressure) in order to avoid the instability region while maintaining a reasonably high level of reactor productivity.

It could be well argued that the mixed model studied in this paper is too simple to represent the actual operation of the membrane reactor which should be modeled as a distributed parameter system (represented by partial differential equations). A rigorous modeling could include a heterogeneous fixed bed membrane reactor such as the one developed in [15] for the production of ultra-pure hydrogen where different reactor configurations were compared. Gomez-Garci et al. [32] also proposed valuable systematic design guidelines for the design and analysis of the membrane reactor for ammonia decomposition. Recently, Murmura et al. [33] also provided a review and critical analysis of the modeling of fixed bed membrane reactors for hydrogen from steam reforming reactions. Nevertheless, the merit of this study is that it is one of few studies, along with [20], that have tackled the instability that could be found in the operation of membrane reactors for ammonia decomposition to hydrogen. We think that the trends obtained using the simple mixed model can be useful to shed light at least qualitatively on the complex stability behavior expected in such reactors.