CHEMICAL ENGINEERING TRANSACTIONS

VOL. 57, 2017

A publication of

The Italian Association

of Chemical Engineering
Online at www.aidic.it/cet

Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš, Laura Piazza, Serafim Bakalis
Copyright © 2017, AIDIC Servizi S.r.l.

ISBN 978-88-95608- 48-8; ISSN 2283-9216

Progress on Modeling and Design of Membrane Reactors for
Hydrogen Production

Maria Anna Murmura*, Stefano Cerbelli, Maria Cristina Annesini
Dip. di Ingegneria Chimica, Materiali Ambiente, Università di Roma “Sapienza”, Via Eudossiana 18 – 00184 Rome (Italy)
mariaanna.murmura@uniroma1.it

This paper presents an overview of recent research carried out by the authors on the development and
analysis of mathematical models describing hydrogen production in membrane reactors. The case considered
is that of methane steam reforming (SR) in a reactor with the typical double pipe configuration, in which a
hydrogen-permeable membrane is present on the outer wall of the innermost tube. The model developed
accounts for the rate of reaction, convective and dispersive transport in the axial and radial directions, and
hydrogen permeation across the membrane. Density variations with pressure and gas composition have been
accounted for, leading to a full coupling of mass and momentum transport. Different geometric aspect ratios
have also been studied to assess the influence of catalyst volume on the overall performance of the system.
The presence of two distinct transport regimes, in which hydrogen permeation is limited either by transport
within the packed bed or permeation across the membrane, has been identified, along with the operating
conditions that determine their range of existence. This has allowed the development of a simplified model,
valid under the hypothesis that the reaction is fast compared to transport. In the permeation-controlled regime,
the permeate flow rate and recovery may be found by solving a set of two PDEs, whereas an analytical
solution is available for the transport-controlled regime. The main steps and observations that have brought to
the development of the simplified model are presented, along with a guide to its implementation.

1. Introduction

The study of integrated membrane reactors, in which one of the products is selectively removed from the
reacting volume as it is being produced, has been attracting much interest. In particular, the use of membrane
reactors for the production of H2 allows to shift the equilibrium of the reaction to high values of feed conversion
while maintaining a fairly low operating temperature even when endothermic reactions are carried out. In order
to adequately design a membrane reactor, several different phenomena must be taken into account, namely
the reaction itself, mass transport within the reactor packed bed, and hydrogen permeation across the
membrane. These mechanisms depend on temperature, pressure, feed flow rate and composition, reactor
dimensions, membrane permeance, catalyst activity, and mass dispersion. The high number of operating
parameters and the wide range of values that they can assume, give an indication of the complexity of the
problem and of the difficulty in designing the system. Previous work carried out by the authors of the present
article has shown that it is possible to describe the system in terms of four dimensionless parameters
(Murmura et al. 2015, 2016a, 2017a) and that integral quantities, such as the H2 permeate flow rate and yield,
can be obtained through a simplified model, valid when the characteristic time of the reaction is low compared
to that of mass transport (Murmura et al. 2016b). In the present work a procedure for a “first approximation”
design of membrane reactors using a simplified model is described. The results of this work are valid when the
temperature of the reactor may be considered to be uniform and when the characteristic time of the reaction is
low compared to the characteristic time of mass transport.

2. Problem description and main equations

The system considered consists of a membrane reactor formed by two co-axial cylinder. The feed enters the
annular volume between the two cylinders, in which the catalyst bed is packed. H2 produced inside the reactor

DOI: 10.3303/CET1757154

Please cite this article as: Murmura M.A., Cerbelli S., Annesini M.C., 2017, Progress on modeling and design of membrane reactors for
hydrogen production, Chemical Engineering Transactions, 57, 919-924 DOI: 10.3303/CET1757154

919

permeates across the selective membrane, placed on the outer wall of the innermost tube, and is collected in
the inner cylindrical volume. The feed composition has been evaluated by considering the equilibrium
conditions of methane SR at a fixed temperature and steam to carbon ratio. In other terms, the presence of
pre-reforming unit before the membrane reactor has been considered. The system has been described
through momentum and mass balance equations. The former has been expressed through Darcy’s law for

flow in porous media, accounting for both radial and axial components of velocity. Mass balance equations
have instead been written by considering both axial and radial convection and dispersion, and the rate of the
SR reaction. This has been expressed through the relation proposed by Wei and Iglesia (2004). Density
variations with pressure and gas composition have been considered, leading to full coupling between mass
and momentum transport. With regards to H2 permeation, Sieverts’ law has been applied, according to which
the permeating flux is proportional to the difference between the square roots of the partial pressure of H2 at
the membrane surface in the retentate and permeate sides, respectively (Ward and Dao, 1999). The model
has been developed under the assumptions of uniform temperature, negligible H2 partial pressure in the
permeate side, and infinite selectivity of the membrane towards H2. The justification of the hypotheses and a
thorough explanation of the equations used are reported in previous papers (Murmura et al., 2017); here, the
dimensionless formulation is reported. The parameters used to develop the dimensionless form of the problem
are the inner reactor radius, R1, as reference length; the inlet gas velocity, U, as characteristic velocity; the
atmospheric pressure, Patm, as reference pressure, the molecular weight of hydrogen, Wh, as reference mass;
the molecular diffusivity, Dm, as reference diffusion constant. The continuity and mass balance equations thus
become

∇ ∙ (−𝑓�̃�∇̃�̃�) = 0 ; ∇ ∙ (−𝛽𝜔𝑖 𝑓�̃�∇̃�̃� −
1

𝑃𝑒
𝑓�̃��̃� ∙ ∇̃𝜔𝑖 ) = 𝜈𝑖

𝑊𝑖

𝑊𝑚
𝐷𝑎𝑓�̃�𝜔𝑚 (1 − 𝜂) (1)

where �̃� = 𝑃 𝑃𝑎𝑡𝑚⁄ , 𝑓 = 𝑓 𝑊ℎ⁄ , and ∇̃= 𝑅1∇ and f is the average molecular weight of the mixture. The other
variables appearing in Eq. (2) are 𝜔𝑖 , Wi, and 𝜈𝑖 which represent the mass fraction, molecular weight, and
stoichiometric coefficient, respectively, of the i-th component. The dimensionless groups introduced are

Pe =
𝑈𝑅1

𝒟
, Da = ℛ𝑇𝑘𝑅1

𝑈
, 𝛽 = 𝑘𝑃𝑎𝑡𝑚

𝜇𝑅1𝑈
(2)

1 − 𝜂 = 1 −
𝑊𝑤

2𝑊𝑚

𝑊ℎ
2𝑊𝑐

�̃�2𝑓2

𝐾𝑒𝑞

𝜔ℎ
4 𝜔𝑐

𝜔𝑚 𝜔𝑤
1

(3)

The parameter (1 − 𝜂) accounts for the distance of the system from equilibrium conditions; 𝛽 measures
pressure drops in packed bed, which, in common reactor configurations, are found to be always negligible
(Murmura et al., 2016a) and whose influence has therefore not been studied extensively. The boundary
conditions are those of fixed inlet gas velocity and composition and fixed outlet pressure (Eqs. (4) and (5)),
zero flux in correspondence of the outer reactor wall (Eq. (6)) and hydrogen permeation in correspondence of
the membrane (Eq. (7))

in �̃� = 0: 𝜔𝑖 = 𝜔𝑖
0 ; �̃�𝑧 = 𝑣𝑧 𝑈⁄ = 1 (4)

in �̃� = 𝐿/𝑅1: �̃� = 𝑃𝐿 𝑃𝑎𝑡𝑚 = 𝛼⁄ ;
𝜕𝜔𝑖

𝜕𝑧
= 0 (5)

in �̃� = 𝜎 ∇̃�̃� ∙ 𝐧 = 0 ; ∇̃𝜔𝑖 ∙ 𝐧 = 0 (6)

in �̃� = 1 (𝜔ℎ 𝛽∇̃�̃� +
1

𝑃𝑒
�̃� ∙ ∇̃𝜔ℎ ) ∙ 𝐧 = −𝛾√

𝜔ℎ

𝑓�̃�
; −�̃� ∙ 𝐧 = −𝛾√

𝜔ℎ

𝑓�̃�
(7)

where 𝜎 is the ratio between the outer and inner radii, n is the unit vector directed towards the reactor’s axis, 𝛼
is the dimensionless operating pressure, and 𝛾 is a parameter accounting for membrane permeability

𝛾 = ℙ𝑚
ℛ𝑇𝑃𝑎𝑡𝑚

−1/2

𝑊ℎ

𝑈
(8)

For a given reactor geometry, the system may therefore be studied through the parameters Pe, Da, γ, and α.
The performance of the system has been studied in terms of

1. overall permeate flow rate of H2 (Πℎ )
2. separator-based yield (from here on “yield”): the ratio between the permeate and inlet flow rate of H2

(Ψ𝑠)
3. reactor-based recovery (from here on “recovery”): the ratio between the permeate flow rate of H2

and the sum of the inlet flow rate and production rate of H2 (Ψ𝑟 ).

920

3. Overview of main results

The model presented above has been studied over a wide range of operating conditions. It was found that
density variations, often neglected, are significant and cause non-negligible convective fluxes in the radial
direction. This phenomenon is at the basis of the concentration-polarization often observed when operating
membrane reactors. The shape of the concentration profiles of H2, and consequently of methane and the
other components, switches between two regimes: one in which H2 permeation is limited by transport across
the packed bed, and one in which it is limited by the membrane itself. When working with a fixed inlet gas
velocity, the following observations can be made:

1. At low P, permeation is limited by transport within the packed bed. As pressure increases above a
threshold value permeation becomes limited by the membrane.

2. The pressure value in correspondence of which the switch occurs, referred to as “critical pressure”,
depends on Pe and membrane permeability, defined through the parameter γ. In particular, the
critical pressure increases with increasing values of γ and Pe (Murmura et al., 2017a).

3. When all other parameters are fixed, the critical pressure represents the value in correspondence of
which it is convenient to operate, since further increases in pressure bring only a negligible increase
in hydrogen recovery and a decrease in the yield (Murmura et al., 2017a).

4. In the transport-controlled regime (low P) the H2 concentration at the membrane reaches values
close to zero and a strong concentration gradient is observable in proximity of the membrane. When
the characteristic time of the reaction is low, CH4 is almost immediately converted to compensate H2
permeation and its concentration also approaches zero (Murmura et al., 2016b).

5. In the membrane-controlled regime (high P) the radial concentration profiles remain almost uniform,
and the behavior of the reactor approaches that of a plug flow reactor (Murmura et al., 2016b).

6. If the curves depicting the yield as a function of pressure are normalized with respect to the critical
pressure and the maximum yield, they all collapse on the same master curve, regardless of the
operating parameters or reactor dimension (Murmura et al., 2017b).

Figure 1: Permeate flow rate and hydrogen yield at Pe=10, Da=10, and γ=1. The continuous and dashed lines

indicate the transport-controlled and membrane-controlled regimes, respectively. The critical pressure, αc has

been indicated.

Figure 1 shows the behavior of the permeate flow rate and H2 yield. The transition between the two regimes
has been evidenced by the use of continuous lines for the transport-controlled regime and dashed lines for the
membrane controlled-regime. Although the analysis carried out at fixed velocity has been useful to identify and
understand the two transport regimes, in most cases of practical interest the inlet flow rate is a fixed
parameter. In these conditions, when pressure is increased the density of the gaseous mixture increases and,
for continuity, the velocity must decrease. The membrane-controlled regime thus appears at low P whereas
the transport controlled regime appears at higher P. The transition is less evident when studying the behavior

921

of the permeate flow rate, but is very noticeable when looking at the recovery. In addition, when the flow rate
is constant, the identification of the critical pressure as the optimal operating pressure becomes very visible,
as increases of this parameter beyond 𝛼𝑐 bring no benefit to the H2 recovery. Figure 2 shows hydrogen
concentration profiles for three axial coordinates (z = 2, 5, and 8) evaluated in same conditions reported in
Figure 1 for α values of 1 (left) and 20 (right).

Figure 2: Hydrogen mass fraction along the reactor’s cross-section and for three axial coordinates in the

transport-controlled (left) and membrane-controlled (right) regimes.

4. Simplified model: development and use

From this starting point, two distinct sets of equations have been developed to evaluate the overall permeate
flow rate in the transport and membrane controlled regimes (Murmura et al., 2016b). In the following
paragraphs, all variables are intended in their dimensionless forms. The main assumptions made in the
development of the simplified model are:

1. Axial dispersion is negligible with respect to axial convection;
2. Radial convection is negligible with respect to radial dispersion;
3. Average molecular mass is constant in the entire domain;
4. Pressure drops are negligible.

The mass balance equations therefore become

−
𝜕𝜌𝑖
𝜕𝑧

+
𝒟𝑟𝑟
𝑃𝑒

𝜌𝑡𝑜𝑡 (
1

𝑟

𝜕𝜔𝑖
𝜕𝑟

+
𝜕2𝜔𝑖
𝜕𝑟 2

) + 𝜈𝑖 𝐷𝑎𝜌𝑚 (1 − 𝜂) = 0 (9)

where 𝜈𝑖 is the stoichiometric coefficients of the i-th component, corrected to account for the fact that balance
equations are expressed in terms of mass.
Additional assumptions have been tailored to the limiting regime. In conditions in which the membrane limits
hydrogen permeation, the radial concentration gradients may be neglected and the system may be described
through a plug-flow reactor model. The mass balance equations for hydrogen and methane thus become

−
𝜕𝜌𝑚

𝜕𝑧
+ 𝜈𝑚 𝐷𝑎

∗𝜌𝑚 (1 − 𝜂) = 0 ; −
𝜕𝜌ℎ

𝜕𝑧
+ 𝜈ℎ 𝐷𝑎

∗𝜌ℎ (1 − 𝜂) = √𝜌ℎ (10)

where 𝐷𝑎∗ = 𝐷𝑎
2𝛾

(𝜎 2 − 1).

The CO2 and H2O densities can be expressed in terms of CH4 density, taking into account the stoichiometry of
the reaction

𝜌𝑐 = 𝜌𝑐
0 −

𝜈𝑐

𝜈𝑚
(𝜌𝑚

0 − 𝜌𝑚 ) ; 𝜌𝑤 = 𝜌𝑤0 −
𝜈𝑤

𝜈𝑚
(𝜌𝑚

0 − 𝜌𝑚 ) (11)

and the total amount of H2 permeating across the membrane may be determined by an overall mass balance

922

Πℎ = 𝜋(𝜎
2 − 1) [𝜌ℎ

0 − 𝜌ℎ
𝑜𝑢𝑡 −

𝜈ℎ
𝜈𝑚

(𝜌𝑚
0 − 𝜌𝑚

𝑜𝑢𝑡)] (12)

On the other hand, when hydrogen permeation is limited by transport within the packed bed a simplified model
has been developed for the condition, often verified, in which the characteristic time of reaction is significantly
lower than the characteristic time of transport. In this case, the problem can be reduced to a pure-transport
problem, in which the mass fractions of all components are tied by the reaction equilibrium conditions. It
should be noted that the imposition of this constraint, along with the obvious consideration that the sum of the
mass fractions must be equal to 1, results in a problem that can be fully described by determining the behavior
of the mass fractions of only two of the four components present in the system.
When the membrane resistance is negligible, the hydrogen concentration in correspondence of the membrane
can be assumed to be equal to zero. Given the additional constraint of equilibrium conditions, the methane
concentration must also be equal to zero. In these conditions the problem may be studied in terms of an
auxiliary variable, resulting from the linear combination of the hydrogen and methane densities

Ω = −𝜈𝑚 𝜌ℎ + 𝜈ℎ 𝜌𝑚 (13)

−
𝜕Ω

𝜕𝑧
+ 𝜀 (

𝑟

𝜕Ω

𝜕𝑟
+

𝜕2Ω

𝜕𝑟 2
) = 0

(14)

where 𝜀 = 𝒟𝑟𝑟
𝑃𝑒

and the boundary conditions are the same as those reported in Eqs. (4)-(6) for the inlet, outlet,
and outer wall. In correspondence of the membrane wall, on the other hand, one can set

Ω = 0 (15)

The solution to this problem may be expressed through Bessel’s function as

Ω = Ω0 ∑
1

𝑁(𝜆𝑛)
𝑒𝑥𝑝(−𝜀𝜆𝑛

2 𝑧)𝑅0(𝜆𝑛, 𝑟) ∫ 𝑟𝑅0(𝜆𝑛, 𝑟)

𝜎

𝑟=1

𝑑𝑟

∞

𝑛=1

(16)

where

𝑁(𝜆𝑛)
=

𝜋 2

𝜆𝑛
2 𝐽0

2(𝜆𝑛)

𝐽0
2(𝜆𝑛) − 𝐽1

2(𝜆𝑛𝜎)
(17)

𝑅0(𝜆𝑛, 𝑟) = 𝐽0(𝜆𝑛𝑟)𝑌′0(𝜆𝑛𝜎) − 𝐽′0(𝜆𝑛𝜎)𝑌0(𝜆𝑛𝑟) (18)

and 𝜆𝑛 are the roots of

𝐽0(𝜆𝑛)𝑌′0(𝜆𝑛𝜎) − 𝐽′0(𝜆𝑛𝜎)𝑌0(𝜆𝑛) (19)

These expressions are valid when studying an annular domain and when the boundary conditions are of zero
concentration on the inner wall and zero flux on the outer wall.
The local permeating hydrogen flux can be determined through

𝝅𝒉(𝒛) = 𝜺
𝝏𝝆𝒉
𝝏𝒓

|
𝒓=𝟏,𝒛

= 𝜺𝛀𝟎 ∑
𝟏

𝑵(𝝀𝒏)
𝒆𝒙𝒑(−𝜺𝝀𝒏

𝟐 𝒛)𝑹′𝟎(𝝀𝒏, 𝒓) ∫ 𝒓𝑹𝟎(𝝀𝒏, 𝒓)

𝝈

𝒓=𝟏

𝒅𝒓

∞

𝒏=𝟏

(20)

Given that the inner radius of the reactor is equal to 1, in dimensionless terms, the total mass flow rate of
permeating hydrogen is given by

𝚷𝒉 = 𝟐𝝅𝛀
𝟎 ∑

𝟏

𝑵(𝝀𝒏)

𝟏

𝝀𝒏
𝟐

[𝟏 − 𝒆𝒙𝒑(−𝜺𝝀𝒏
𝟐 𝒛)]𝑹′𝟎(𝝀𝒏, 𝒓) ∫ 𝒓𝑹𝟎(𝝀𝒏, 𝒓)

𝝈

𝒓=𝟏

𝒅𝒓

∞

𝒏=𝟏

= 𝟐𝝅𝛀𝟎𝒈(𝑷𝒆) (21)

where 𝒈(𝑷𝒆) is a function that accounts for the dispersive contribution to mass transport and depends
exclusively on the Peclet number and the geometric dimensions of the reactor. Once the Bessel function to
use has been identified, it may be used to find the solution to a given problem, with a good degree of
approximation. The remaining difficulties lie in the identification of the first root of the Bessel function of Eq.
(23), 𝝀𝟏, and in solving the integral introduced in Eq. (20). With regards to the roots, their value depends on
the 𝝈 and has been reported in Table 1 for 𝝈 ∈ [𝟏. 𝟏, 𝟐]. It should be emphasized that the correct identification
of the first root is of paramount importance for the correct evaluation of the permeate.

923

Table 1: First roots of the Bessel function reported in Eq. (19) for 𝜎 ∈ [1.1,2]

𝜎 𝜆1
1.1 15.378
1.2 7.5220
1.3 6.5000
1.4 4.0500
1.5 2.9000
1.6 2.3755
1.7 2.0118
1.8 1.7388
1.9 1.5420
2.0 1.3068

The integral in Eq. (20) is given by

∫ 𝒓𝑹𝟎(𝝀𝒏, 𝒓)

𝝈

𝒓=𝟏

𝒅𝒓 = −
𝟏

𝝀𝒏
(𝒀𝟏(𝝀𝒏)𝑱𝟏(𝝀𝒏𝝈) − 𝒀𝟏(𝝀𝒏𝝈)𝑱𝟏(𝝀𝒏)) (22)

In order to verify that the problem has been solved correctly, Eq. (23) shown below should be verified.

∑
𝟏

𝑵(𝝀𝒏)
𝑹𝟎(𝝀𝒏, 𝒓) ∫ 𝒓𝑹𝟎(𝝀𝒏, 𝒓)

𝝈

𝒓=𝟏

𝒅𝒓 = 𝟏

∞

𝒏=𝟏

(23)

which corresponds to the condition 𝛀 = 𝛀𝟎 in 𝒛 = 𝟎, from Eq. (16). It has been observed that an error lower
than 2% must be achieved to guarantee the correctness of the solution. Equations (12) and (21) yield the
permeate flow rate obtained in conditions in which hydrogen permeation is limited by the membrane and by
transport, respectively. In order to determine, for a given operating condition, the actual permeate flow rate,
the smallest of the two values should be taken. In correspondence of the critical pressure, the two expressions
will result in the same value.

5. Conclusions

This work presents an overview of studies carried out on modeling and design of membrane reactors for
hydrogen production. In particular, the evolution from a fully coupled model to a simplified, pure transport,
model is explained and justified through observations made on the permeate flow rate, hydrogen yield and
recovery, and concentration profiles. In addition, the detailed procedure for the use of the simplified model has
been reported.

Reference

Murmura, M.A., Turchetti, L., Augelletti, R., Annesini, M.C., Cerbelli, S., 2015, How does radial convection
influence the performance of membrane module for gas separation processes? Chem. Eng. Trans. 43,
1063-1068.

Murmura, M.A., Cerbelli, S., Turchetti, L., Annesini, M.C., 2016a, Transport-permeation regimes in an annular
membrane separator for hydrogen purification. J. Membr. Sci. 503, 199-211.

Murmura, M.A., Cerbelli, S., Annesini, M.C., 2016b, An equilibrium theory for catalytic steam reforming in
membrane reactors. Chem. Eng. Sci. 160, 291-303

Murmura, M.A., Cerbelli, S., Annesini, M.C., 2017a, Transport-reaction-permeation regimes in catalytic
membrane reactors for hydrogen production. The steam reforming of methane as a case study. Chem.
Eng. Sci. 162, 88-103

Murmura, M.A., Cerbelli, S., Annesini, M.C., 2017b, Modeling and optimization of hydrogen yield in membrane
steam reforming reactors. Can. J. Chem. Eng. Accepted for publication.

Ward, T.L., Dao, T., 1999, Model of hydrogen permeation behavior in palladium membranes. J. Membr. Sci.
153, 211-231.

Wei, J., Iglesia, E., 2004, Isotopic and kinetic assessment of the mechanism of reactions of CH4 with CO2 or
H2O to form synthesis gas and carbon on nicked catalysts. J. Catal. 224, 370-383.

924