Effect of Asymmetric Membrane Structure on Hydrogen Transport Resistance and Performance of a Catalytic Membrane Reactor for Ethanol Steam Reforming

Abstract

The performance of catalytic membrane reactors (CMRs) depends on the specific details of interactions at different levels between catalytic and separation parts. A clear understanding of decisive factors affecting their operational parameters can be provided via mathematical simulations. In the present paper, main results of numerical studies of ethanol steam reforming, followed by downstream hydrogen permeation through an asymmetric supported membrane, are reported. The membrane module consists of a thin selective layer supported on a substrate with graded porous structure. One-dimensional isothermal reaction–transport model for the CMR has been developed, and its validation has been carried out by using performance data from a lab-scale reactor with a disk-shaped membrane. Simulations demonstrate the model’s capabilities to analyze local concentrations gradients, as required to provide accurate estimates of the relationship between structure–property–performance. It was shown that transport properties of multilayer asymmetric membranes are highly related to the structural properties of each single layer.

1. Introduction

With the growing concerns about environmental issues, catalytic reforming of fuels to hydrogen-rich gas synergistically coupled with membrane technology has become a huge focus of attention, being used increasingly in a broad range of applications. Integrated chemical reaction and separation options facilitate process miniaturization, continuous operation and energy saving. This intensification technique is expected to be a promising route in creating the sustainable green chemistry-driven energy technologies for small-scale applications. Catalytic membrane reactors are able to provide higher fuel conversion with the advantage of producing a very pure hydrogen stream supply [1,2,3,4,5].

In fact, the aim of the membrane reactors technology is to separate hydrogen from feed streams that consist of various species, depending on feedstock and catalytic chemistry. Bioethanol as a renewable source is considered a promising candidate for hydrogen generation in relation to relatively high hydrogen content, being nontoxic, safe and easy to transport, as well [6,7,8,9]. While hydrogen is primarily formed from ethanol, the product gas can also include CO, CO₂, CH₄ and/or C (coking), etc. The reaction conditions, such as the fuel-to-oxidant (steam, CO₂, etc.) molar ratio, nature of catalyst and temperature, have a significant influence on the composition of reaction products [10,11].

In the practical design and operating decisions, such multifunctional reactor configuration consists of two reactor volumes, which are the reaction (feed-side) compartment followed by a permeate sweep-side zone. A functionally selective diffusion barrier separating the reactor compartments enables hydrogen produced by a catalytic reaction to migrate into the side of a lower hydrogen concentration and be carried out by a sweeping gas, while the retentate stream leaves the feed-side compartment. The active area of a suitable catalyst in the feed-side must be much larger than the geometric area of the functional layer of a membrane module. From the point of view of catalyst effectiveness, packed beds consisting of catalyst particles are preferably used. In addition, placing the membrane immediately downstream from the catalyst bed results in a simpler overall design. Other important advantages of such a catalyst bed arrangement are its simplicity in construction and well-established and validated models for its design and scale-up. Additionally, when the catalyst is kept in a fixed position, any damage of the membranes can only happen when loading and unloading the catalyst from the reactor [4,12,13]. However, the catalyst packed bed in CMRs makes possible some limitations to hydrogen transport between the bulk of the catalytic bed (where hydrogen-rich gas is produced) and the membrane surface [14,15].

The membrane, itself, may also be of potential benefit in a chemical transformation. Indeed, membrane characteristics may affect the overall rate of a given reaction, products selectivity and yield. Many efforts have been dedicated to developing high-flux membranes showing long-term stability at practically relevant conditions. The membrane must not rupture, split, crack, creep or otherwise develop defects that allow the significant, nonselective flow of gases from the feed stream into the permeate stream. Further optimization and new developments in membrane modules architecture are required for integration and intensification of processes in a CMR. In particular, intensive research focused on the relationship between structure–property–performance of the reactors is necessary [16,17].

The desired mechanically robust membrane modules may be obtained via an asymmetric structure, which includes a thin, dense, functional (permselective) layer supported on a thick, porous layer. The efficiency of gas permeation through the membrane is typically determined by two main parameters: selectivity and permeance. The dense layer symbolizes the selective barrier in the membrane structure. Next, gas diffusion layers of a support also serve to conduct the permeating gas to be swept away at the outgoing side of the membrane module. The porosity-graded, multilayered structure of gas-permeable supports with a defect-free interfaces between layers is claimed to possess good thermal stability, chemical resistance, and a high compressive strength, all of which are favorable for use in a wide variety of applications [18,19,20,21]. Asymmetric membranes exhibit significant commercial potential due to a high permeation combined with the excellent strength and durability [22,23,24].

However, beyond the sufficient mechanical strength, a membrane module assembly should be adapted to provide a low resistance to mass transport to improve the driving force across the membrane surface. In turn, transport properties of multilayered asymmetric membranes are highly related to the structural properties of each single layer. If the combination of porous substrates with a thin dense layer may overcome a reduction in the driving force due to polarization, proton diffusion across the membrane may no longer become the rate limiting step for the overall hydrogen flux, and catalytic reaction may play a more essential role to further increase the hydrogen flux through the membrane module [25].

The overall hydrogen flux depends on, and is limited by, reactor design, reaction kinetics, the partial pressure driving force and the mass transfer resistance distribution. The driving force for permeation is also increased by using an inert sweep gas on the permeate membrane side, thus reducing the H₂ partial pressure at the outlet of the membrane module assembly. Concentration profiles are of great interest to determine local resistance to hydrogen permeation through the membranes of an asymmetric configuration. To a great extent, overall permeate flux driven by concentration gradients is controlled by parameters of asymmetric membrane module such as thickness of layers, pore size, pore volume, pore distribution and other performance requirements [26,27,28,29].

The asymmetric supported membrane based on Ni + Cu/Nd_5.5WO_11.25−δ mixed proton–electron-conducting nanocomposites have been developed and integrated in a lab-scale tubular reactor. Selective separation of hydrogen from the product flow of ethanol steam reforming in the 5 wt.% Ni + 1 wt. % Ru/Sm_0.15Pr_0.15Ce_0.35Zr_0.3O_2−δ catalyst layer has been successfully tested at ambient pressure and the temperature range of 700–900 °C. The best results have been obtained at 900 °C and feed of ethanol/H₂O mixture in Ar at steam-to-carbon ratio of 2. An overall hydrogen flux was achieved to be about 1.31 N mL min⁻¹ cm⁻² [30,31,32,33,34].

A disk-shaped membrane is the most popular design commonly used for laboratory studies due to the ease of fabrication and convenience for investigating the fundamental features of membrane performance. In the case of the disk-shaped membrane, a membrane module can be mounted between two vertical ceramic or quartz tubes and then placed in a bigger quartz tube. Pyrex or gold gaskets are usually used to obtain an effective seal between the disk and the walls of the inner tubing at high temperatures by placing the assembly in compression with the use of spring clamps. The catalyst is usually packed on the membrane or coated on the feed-side surface [35]. By using computational fluid dynamics, polarization effects in the vicinity of the disc-shaped ionic–electronic conductor were shown to be very important for analysis of dependence of the permeation flux on the operating conditions and should be taken into account in numerical simulations [36,37].

The scope of this paper is to describe the catalytic reaction of ethanol steam reforming and hydrogen permeation through the asymmetric supported membrane in a catalytic membrane reactor by considering reaction rates and transport through the catalyst–membrane assembly, as well as boundary layer effects. The formulated one-dimensional reaction–transport model for the constituent layers of the catalyst–membrane assembly, together with a Sieverts’ equation for the dense functional layer, has been implemented in COMSOL Multiphysics. The mathematical model was then extensively verified with experimental data for wide range of operating conditions and shown to provide good agreements with the data with reasonable parameter values. The simulations demonstrate the model’s capabilities to study resistance features on hydrogen permeation. The performance of the CMR, in terms of efficiency of catalytic process and hydrogen recovery, was also studied.

Most significant results are detailed hereafter.

2. Materials and Methods

Experimental studies were carried out using a specially built setup. Detailed information about the experimental setup and disk-shaped reactor with axis-symmetric flows has been reported in previous papers [31,32,33]. In the experimental reactor, the assembly of a catalyst—an asymmetric membrane based on Ni + Cu/Nd_5.5WO_11.25−δ mixed proton–electron conducting nanocomposites [30] was placed between the feed and sweep-side compartments. The catalyst bed, consisting of particles 5 wt.% Ni + 1 wt.% Ru/Sm_0.15Pr_0.15Ce_0.35Zr_0.3O_2−δ, was located above the membrane. Before the test procedure, the reactor assembly was heated in Ar from room temperature to 900 °C, followed by purging for 1 h. By using a dual bubbler system, the ethanol–steam–argon gas mixture was introduced into the feed-side compartment in a direction substantially perpendicularly to the catalyst surface, while the exiting retentate stream went upward through an external annulus. Argon was applied as a sweep gas at the permeate-side compartment to obtain the partial pressure driving force required for hydrogen permeation at the ambient pressure. Argon entered from an inlet tube placed in the vertical axis, and hydrogen in sweep gas exited the permeate zone at the external area upward. Composition of the outlet streams was monitored by the TEST-1 gas analyzer (Bonair, Russia).

2.1. Membrane Morphology Characterization

A membrane module consisting of a gas-tight nanocomposite functional coating deposited on a gas-permeable substrate was used in the experiments. The procedure of integration of the active NiCu (30 wt.%)—Nd_5.5WO_11.25−δ dense layer with a support requires multiple coating processes when gas-tight coatings are desired, and it is clearly described elsewhere [31,38]. Apart from mechanical strength, supports need to satisfy some other performance requirements. Thus, such parameters as thickness, pore volume, pore size distribution, etc., affect the overall flux through an asymmetric membrane to a great extent [21,26].

Powder Metallurgy Institute (Minsk, Belarus) provided functionally graded foam substrates. The method for making gas-permeable structures involves applying two thin, low-porosity layers onto a thick Ni/Al foam. The tough foam is made from alumina–silica ceramics. In this process, 30 PPI polymeric foam is used as template. Macroporous ceramics obtained by 3 times impregnation of foam polyurethane, followed by centrifugation, drying and sintering at 1350 °C. Moreover, Ni covering of alumina–silica ceramics is done by electrolysis, followed by drying and sintering at 1000 °C in cracked ammonia. The Ni:Al₂O₃ weight ratio obtained was 2–2.2:1 [39]. Formation of two thin layers onto the Ni/Al substrate was performed by the use of dual doctor blades in series. Here, in the doctor blading precision coating process, a well-mixed slurry consisting of a suspension of Ni spherical particles along with polyvinyl alcohol as binder and glycerol as plasticizer (<100 µm and 100–200 µm for first and second layer, respectively) is applied. The compaction of layers was performed via sintering of the specimen at 1000 °C in cracked ammonia, followed by in-pack aluminizing, as described in [40].

As it has been shown by the use of stereology-based image analysis, the resultant gas-permeable substrate exhibits composite structure, in which the top porous layer has smaller pores in diameter than those in the intermediate porous layer on the foam support (Figure 1).

Indeed, the structural differences among porous layers would affect the distribution of hydrogen transfer resistance through the membrane module. Epoxide resin casting method is used to allow porous structures to be easy visualized and quantified [41]. A specimen is subjected to embedding procedure in the Struers’ EpoxiFix-20 epoxide resin by using Struers filling machine. The epoxy preparation is then poured over the sample and allowed to cure for 24 h. Finally, after hardening, the sample is removed from the form. Grinding is performed by using Struers manual grinding machine, and then the surface was subsequently smoothed by 200, 600 and 1200 diamond grinding wheels and polished by 9 and 3 µm diamond suspensions using MD-Mol canvas on StruersRotopol-35 automatic grinding–polishing machine. Finally, the sample is washed with distilled water and ethanol and then dried.

The morphology of the prepared membrane module with graded porosity and pore sizes is visualized in two magnification scales (Figure 2). Its layer-by-layer assembled asymmetrical structure is demonstrated by a cross section photo on the Figure 2a, while the upper membrane layers (a reference area on the photo) are shown in a scanning electron microscope (SEM) micrograph on the Figure 2b. For preparing SEM images, Inkscape software was used. As can be seen from Figure 2b the permselective layer with a thickness of around 150 µm is a gas-tight layer at which only few closed pores can be seen, thus not allowing a significant amount of nonselective flow of gases from the feed stream into the permeate stream. A good adherence between the permselective layer and the next porous layer is observed.

In terms of pore size distributions and porosity, the three-layer structure of the porous support is clearly quantified and designated in

Table 1. Morphological and structural characteristics of the prepared asymmetric membrane module.

Layer	Composition	Thickness (µm)	True Density (g cm−3)	Particle Size b (µm)	Pore Diameter b (µm)	Porosity c (%)
Dense layer	NiCu/Nd5.5WO11.25−δ	93.3–115 (center); 194–256 (edge)	6.6	0.045 for NiCu, 0.1–1 for Nd5.5WO11.25−δ	15 (x) 42 (y)	~4
Powder layer	Ni-Al	380–440	~7	65 (x) 81 (y)	12 (x) 11 (y)	12–14
Intermediate layer	Ni-Al	400–1300	5.34	45 (x) 50 (y)	27 (x) 27 (y)	27–32
Foam layer	Al2O3-SiO2 foam with NiAl coating	4500–5000	4.63	2400(x) a 1800 (y) a	1000 (x) 1100 (y)	38–40 83 d

^a Cell diameter. ^b x axis is parallel, and y axis is perpendicular to the membrane surface. ^c Quantification with vector program and by image analysis. ^d Overall porosity, including pores in Al₂O₃.

. Quantification of the constituent particles and pore systems was performed via the image analysis with a Joyce–Loebl Mini Magiscan (Joyce–Loebl, Ltd., Gateshead, UK) computerized image analysis system [42,43] and with a vector program by applying an appropriate spatial calibration in and a series of parallel and perpendicular lines along the image axes in the reference area [44,45]. In this method, extracting numerical porosity data in a selected zone includes the image calibration, fixation, segmentation and selection of necessary measurements. The calibration was carried out by using a standard test object for zooming. The images were converted into electrical digital image signals and stored. The particles were then segmented from background of the image to form a binary image. Once a segmented digital object (particles) was available, porosity, pore and particles sizes could be computed by counting the area of particles exposed to the pore space in the porous media image.

2.2. Approach to Mathematical Modeling and Simulations

Hydrodynamics of flow and effect of operational factors on fluid dynamics in a similar reactor design were studied by using the computational fluid dynamic simulations in [36]. If the thickness of the membrane disk module is small compared to its diameter, deflection of flow streamlines and concentration gradients in radial direction, at the proximity of the membrane, was shown to be rather low, while the gradients in the axial direction are sufficient. Therefore, a significant hydrogen–partial pressure gradient exists only across the membrane, and, due to the radial symmetry, the flow can be considered as one-dimensional. With the purpose to capture concentration distributions in axial direction, a computationally efficient 1-D modeling approach, neglecting the diffusion along radius, was applied. Ideal gas, negligible radial gradients and negligible axial pressure drop are assumed, while variations in total molar density were not neglected in the developed model. The flow field in the feed compartment was coupled with the membrane module in its one-dimensional form.

The overall transfer rate of hydrogen depends both on reaction kinetics in the catalyst layer and on the mass transfer resistance through catalyst–membrane assembly. Chemical reactions in the catalyst bed, a hydrogen flux through the dense layer following Sievert’s law and diffusive mass transfer in the constituent layers of gas-permeable support are supposed to occur in series. In this study, an approach is used where fluxes through all constituent layers are coupled by defining interface concentrations and flows on the boundaries inside the catalyst–membrane assembly and combined by means of these interface concentrations. Molar flow change due to reaction and membrane transport is accounted for. The developed model is validated by confronting its predictions with data from experimental studies with disk-shaped reactor [30,31,32]. In addition, experimental data was used to obtain and estimate all necessary parameters to quantify the characteristics for the developed mathematical model.

Dimensions, structural parameters of the CMR under consideration, as well as operating conditions and experimental data considered in simulations, are given in

Table 2. Structural parameters of the constituent layers considered in simulations.

Layer Property	Catalyst	Dense Layer	Powder Layer	Intermediate Layer	Foam Layer
Thickness (mm)	5	0.15	0.4	0.6	4.5
Particle size (mm)	1		0.072	0.061	2.2 Cell diameter
Hydraulic pore diameter (mm)	1.75		0.012	0.027	1.006
Porosity (-)	0.42		0.2	0.4	0.75
Tortuosity (-)	1.37		4.2	3.4	1.42
Volumetric surface area (m2 m−3)	3480		66667	59259	1395.4

and

Table 3. Operating conditions and experimental data considered in simulations.

Parameter	Experiment No
	1	2	3	4	5
Feed-side compartment
Feed composition (NL h−1) a
EtOH	0.3	0.132	0.264	0.396	0.66
Steam	1.2	3.41	3.41	3.41	3.41
Ar	1.5	5.2	5.2	5.6	6.0
H2; CO; CO2 in the dry retentate, (vol.%), at T (°C)
700	23.7;5.3; 4.1	2.4; 0.4; 1.0	4.8; 1.0; 1.6	7.9; 2.3; 2.2	8.7; 2.5; 2.5
800	24.6; 8.3; 3.5	3.0; 0.5; 1.1	6.1; 1.3; 1.9	9.1; 3.1; 2.3	10.7; 3.8; 2.6
900	24.4; 8.1; 4.0	3.4; 0.5; 1.1	6.0; 1.3; 1.8	10.6; 3.8; 2.3	12.1; 4.7; 2.2
Sweep-side compartment
Ar flow rate (NL h−1) a	2	10	10	10	10
Hydrogen, (vol.%), at T (°C)
700	17.7	0.9	2	3.3	3.9
800	19.9	1.2	2.6	4.4	4.9
900	20.8	1.5	2.7	5.0	5.7
Hydrogen permeation flux (mol m−2 s−1) at T (°C)
700	0.01005	0.00203	0.00451	0.00743	0.00879
800	0.01161	0.00255	0.00553	0.00935	0.01104
900	0.01227	0.00319	0.00574	0.01126	0.01284

^a Letter “N” stands in the unit for a flow at STP conditions of temperature of 273.15 K (0 °C) and an absolute pressure of exactly 1 atm (101.325 kPa).

. Model parameters are applied to be closely related to the characterization of the structural properties of the constituent layers (Table 1).

2.2.1. Feed-Side Compartment Model

The steam reforming of ethanol is assumed to include a stepwise reaction scheme [46,47]:

• r₁—Ethanol decomposition (cracking reaction):

$${\mathrm{C}}_2{\mathrm{H}}_5\mathrm{OH}\to {\mathrm{CH}}_4+{\mathrm{H}}_2+\mathrm{CO},{\mathsf{\Delta }}_r{H}_{298\mathrm{K}}^o=49.0\mathrm{kJ}{\mathrm{mol}}^{-1};$$

(1)

• r₂—Steam reforming:

$${\mathrm{CH}}_4+{\mathrm{H}}_2\mathrm{O}\leftrightarrow 3{\mathrm{H}}_2+\mathrm{CO},{\mathsf{\Delta }}_r{H}_{298\mathrm{K}}^o=206.3\mathrm{kJ}{\mathrm{mol}}^{-1};$$

(2)

• r₃—Water–gas shift:

$$\mathrm{CO}+{\mathrm{H}}_2\mathrm{O}\leftrightarrow {\mathrm{CO}}_2+{\mathrm{H}}_2,{\mathsf{\Delta }}_r{H}_{298\mathrm{K}}^o=-41.2\mathrm{kJ}{\mathrm{mol}}^{-1};$$

(3)

• r₄—Complex shift reaction:

$${\mathrm{CH}}_4+2{\mathrm{H}}_2\mathrm{O}\leftrightarrow 4{\mathrm{H}}_2+{\mathrm{CO}}_2,{\mathsf{\Delta }}_r{H}_{298\mathrm{K}}^o=165.1\mathrm{kJ}{\mathrm{mol}}^{-1};$$

(4)

The basic reactions (1)–(4) are considered as dependent on the concentration of the reactants, the temperature and the catalyst external surface area. Reaction kinetics based on per unit catalyst surface were taken from an earlier experimental study of ethanol steam reforming performed over the 5 wt.% Ni + 1 wt.% Ru/Sm_0.15Pr_0.15Ce_0.35Zr_0.3O_2−δ catalyst [47]. Details for observed reaction kinetics with effective parameters and variables and their corresponding descriptions are provided in the Appendix A.

It is believed [48,49] that, in reactors with the separation option, reduced Reynolds numbers (Re) are desired, indicating highly ordered laminar flow and low Péclet number (Pe of less than or around 1), when diffusive transports between phases are preponderant. Such is the case with the experimental CMR operating conditions (Table 3). Thus, the Reynolds number varies up to a maximum of 550, assuring laminar flow conditions in all experiments. Chemical reactions result in production or consumption of species (R_i, see Appendix A), which is modeled as a molar source or sink for the i-th specie. A continuity equation for species molar balance in the feed-side compartment contains partial derivatives of molar flow rates and species concentrations with respect to position, and it is given by the following equation:

$$\epsilon {\rho }_{tot}\frac{\partial x_{i,f}}{\partial t}+\frac{1}{A_m}\frac{\partial (F_f·x_{i,f})}{\partial z}-D_{i,eff,f}{\rho }_{tot}\frac{{\partial }^2{x}_i}{\partial z^2}=S_{V,f}{R}_i\left(\begin{array}{l}i=1(\mathrm{EtOH}),i=2({\mathrm{H}}_2\mathrm{O}),i=3({\mathrm{CH}}_4),\\ i=4(\mathrm{CO}),i=5({\mathrm{CO}}_2),i=6({\mathrm{H}}_2)\end{array}\right)$$

(5)

At steady-state formulation, the term $\epsilon {\rho }_{tot}\frac{\partial x_{i,f}}{\partial t}$ vanishes in Equation (5). By substitution of the conservation equation for a change in the total molar flow rate between no conversion and complete conversion (Appendix B, Equation (A15)) in Equation (5), a molar balance can be defined for each of the components i except for argon:

$$\epsilon {\rho }_{tot}\frac{\partial x_{i,f}}{\partial t}+\frac{1}{A_m}{F}_f\frac{\partial x_{i,f}}{\partial z}+2·S_{V,f}·x_{i,f}·\left(r_1+r_2+r_4\right)-D_{i,eff,f}{\rho }_{tot}\frac{{\partial }^2{x}_i}{\partial z^2}=S_{V,f}{R}_i,$$

(6)

$$x_{\mathrm{Ar}}=1-{\displaystyle \sum _6{x}_i}.$$

(7)

If pressure gradients through the catalyst–membrane assembly are considered to be not significant, the total molar density of the feed mixture is adjusted to ensure ${\rho }_{tot}=\frac{P}{RT}=\mathrm{const}$ at the given temperature. Diffusion hydrogen molecular flux goes from the regions of higher concentration into the regions of lower concentration. Effective area of the permeating surface A_m = 5.3066 × 10⁻⁴ m². It is known [48] that the complexity of the molecular transport processes does not allow a purely theoretical fundamental approach in analysis of the diffusive transport. In the numerical simulations, transport coefficients (the diffusion coefficients, viscosity coefficients, etc.) are calculated from the transport coefficients for the individual species. Therefore, fitting the model to the experimental data (No.1 and Nos.2–5, Table 3) results in different sets of effective diffusivities in the catalyst layer. Molar fluxes are described, then, based on a mixture averaged diffusion coefficient approach.

It is known that, for any surface in which the species being consumed, a concentration gas boundary layer is formed. Because the membrane restricts the passage of species except of hydrogen, their fluxes are zero at the membrane surface: ${D_{i,eff}{\rho }_{tot}\frac{\partial x_{i,f}}{\partial z}|}_{z=h_{cat}}=0$ $i\ne 6({\mathrm{H}}_2)$. The Direchlet boundary condition is applied at inlet of the catalyst layer:

$${x_{i,f}|}_{z=0}=x_{i,f}^0;F_f=F_f^0$$

(8)

The concentration of H₂ at the dense membrane surface becomes also somewhat lower than that in the retentate gas. The difference between the concentrations, namely $x_{{\mathrm{H}}_2,f}$ and $x_{{\mathrm{H}}_2,dm}$, is a driving force of the molecular hydrogen transport through the boundary layer to membrane surface. The hydrogen flux at the catalyst/membrane interface relates to the concentration difference and a mass transfer coefficient ${\beta }_f$, according to:

$${\left(\frac{1}{A_m}·F_f·x_{{\mathrm{H}}_2,f}-D_{{\mathrm{H}}_2,eff}{\rho }_{tot}\frac{\partial x_{{\mathrm{H}}_2,f}}{\partial z}\right)|}_{z=z_{h_{cat}}}=-{\beta }_f·{\rho }_{tot}·\left(x_{{\mathrm{H}}_2,f}-x_{{\mathrm{H}}_2,dm}\right).$$

(9)

Consequently, the boundary conditions are implemented in Comsol package in the following form:

$$\begin{gathered} -\overrightarrow{n}·\left(-c\frac{\partial x\_{\mathrm{H}}_2}{\partial x}-ax\_{\mathrm{H}}_2+\gamma \right)=g-qx\_{\mathrm{H}}_2, \\ c=D_{{\mathrm{H}}_2,eff}{\rho }_{tot}, \\ a=0, \\ \gamma =0, \\ g={\beta }_f·{\rho }_{tot}·x_{{\mathrm{H}}_2,dm}, \\ q=\left({\beta }_f·{\rho }_{tot}+\frac{1}{A_m}·F_f\right). \end{gathered}$$

The effective gas-dense membrane surface mass transfer coefficient ${\beta }_f$ can be calculated by its relationship to the dimensionless Sherwood number, Sh, as in the equation:

$${\beta }_f=\frac{Sh·D_{{\mathrm{H}}_2,eff}}{d_h}.$$

(10)

The Sh number is a measure of the ability to mass transfer. With regard to fluid transport properties in packed beds, the Sherwood number Sh is estimated by its relationship with the dimensionless Schmidt number, Sc, and Reynolds number, Re, via the equation that is valid for Pr = Sc = 0.6–10 [50,51,52]:

$$Sh=0.515R{e}_{eq}^{0.85}S{c}^{0.333}.$$

(11)

Generally, characteristic length scale d_h is the distance at the proximity of confining surface corresponding to the essential change of flow velocity. This definition is equivalent to characteristics of the diffusion boundary layer at the catalyst/membrane interface, which is dependent on the neighboring particles. To all appearance, porosity or void fraction of the catalyst layer is increased at the proximity of the confining membrane surface. Thus, a definite nonrandom trend up to distances of about 2 particle diameters from the wall was observed in the fixed beds with different aspect ratios [53]. They found that the void fraction was increased by factor two (${\epsilon }^{\prime }=2\epsilon$≈0.8). In the simulations, the value of an equivalent hydrodynamic diameter d_h is calculated based on the equivalent diameter of the particles:

$$d_h=\frac{{\epsilon }^{\prime }}{\left(1-{\epsilon }^{\prime }\right)}{d}_p.$$

(12)

The effective molecular diffusion coefficient for hydrogen $D_{{\mathrm{H}}_2,eff}$ (m² s⁻¹) depends on the properties of gas and porous media. For more details on all available features and model parameters, see Appendix C.

2.2.2. Through-the-Membrane Transport Model

Dense Permselective Layer

Permselectivity of a dense layer is derived from the intrinsic properties of the active material, while the driving force for hydrogen transport across the dense membrane is caused by the difference of hydrogen partial pressure on each side.

The Sieverts’ law, which identifies the difference of the hydrogen partial pressure square roots as the permeation driving force, is a temperature-activated phenomena. It is valid with an underlying assumption that the surface coverage is low when interfacial equilibrium is achieved, and the rate-limiting step is an atom diffusion through the membrane. These criteria are satisfied in most cases when the temperature is relatively high [54]. Obviously, the steady state hydrogen flux through the dense layer of thickness h_dm being described by Sielverts’ law is equal to permeate flux through the asymmetric membrane (Equation (9)):

$$J_{{\mathrm{H}}_2}={\beta }_f·{\rho }_{tot}\left(x_{{\mathrm{H}}_2,f}-x_{{\mathrm{H}}_2,dm}\right)=Q_{dm}\left(\sqrt{P_{{\mathrm{H}}_2,dm}}-\sqrt{P_{{\mathrm{H}}_2,pl}}\right).$$

(13)

In the Equation (13), $P_{{\mathrm{H}}_2,dm}=P·x_{dm}$ and $P_{{\mathrm{H}}_2,pl}=P·x_{pl}$ are the hydrogen partial pressures on opposite sides of the dense layer, upstream, on the feed-side and downstream, at the interface with powder layer of the gas-permeable support (Table 2), respectively. The dense membrane permeance Q_dm (mol m⁻² s⁻¹ Pa^−0.5) is evaluated as an Arrhenius-like equation:

$$Q_{dm}=\frac{\Theta }{h_{dm}}\exp \left(-\frac{E_{dm}}{RT}\right).$$

(14)

The permeance and the permeation constant, permeability, $\Theta$ (mol m⁻¹s⁻¹Pa^−0.5), are commonly used as a measure indicating permeation ability. For asymmetric membranes, where the thickness of the selective layer is not well defined, the permeance is more useful in terms of membrane comparison than permeability.

In the following simulations, the apparent activation energy for permeability E_dm is taken as 60 kJ mol⁻¹ (Figure 3). The temperature dependence of the proton conductivity $\left(\sigma \right)$ for the dense NiCu (30 wt.%)—Nd_5.5WO_11.25−δ nanocomposite in moist atmosphere of hydrogen shows three distinct regimes (Figure 3a). This is an evidence for the existence of different predominant proton transfer mechanisms occurring in each temperature range [55].

The corresponding apparent activation energies (Figure 3b) can be calculated formally as in Equation (15) using temperature-dependent conductivity, which follows the Arrhenius relation. The slope of the high-temperature curve (490–700 °C) results in an activation energy of ≈60 kJ mol⁻¹.

$$E_{dn}=-R\frac{\partial \ln \left(\sigma T\right)}{\partial \left(\frac{1}{T}\right)}.$$

(15)

It is known that the hydrogen flux and the hydrogen permeation coefficient can be estimated precisely only for given conditions. Therefore, for predicting the hydrogen flux, the apparent gas permeability is useful to quantify for certain operating regime [56,57]. The catalyst–membrane assembly in the CMR has been tested under many different flow rates (Table 3, experiments No.1 and Nos.2–5), and actual fluid dynamics may affect an overall permeability considerably. Thus, $\Theta$ = 0.014879 (mol m⁻¹s⁻¹Pa^−0.5) is determined for the set No.1, while the magnitude of the apparent gas permeability for the set Nos.2–5 is estimated to be 2.0924 × 10⁻⁴ (mol m⁻¹s⁻¹Pa^−0.5).

Explicit solution for $x_{{\mathrm{H}}_2,dm}$ is attained from the parity Equation (13):

$$x_{{\mathrm{H}}_2,dm}={\left(\frac{-Q_{dm}\sqrt{P}+\sqrt{P{Q}_{dm}{}^2+4·{\beta }_f{}^2·{\rho }_{tot}{}^2·x_{{\mathrm{H}}_2,f}+4·{\beta }_f·{\rho }_{tot}{Q}_{dm}\sqrt{P}\sqrt{x_{{\mathrm{H}}_2,pl}}}}{2{\beta }_f·{\rho }_{tot}}\right)}^2.$$

Powder Layer

The gaseous mixture inside the micro/mesoporous powder layer of 0.4 mm thick (Table 2) is a binary mixture of argon and hydrogen. The penetration of the sweep gas, argon, into the substrate can only be achieved by diffusion. Therefore, diffusion is the only mass transport mechanism in the powder layer, so the mass balance equation for hydrogen transport reads as Equation (16), while Equation (17) is the relation between mole fractions in the gas phase.

$${\epsilon }_{pl}{\rho }_{tot}\frac{\partial x_{{\mathrm{H}}_2,pl}}{\partial t}-{\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},pl}^{eff}\frac{{\partial }^2{x}_{{\mathrm{H}}_2,pl}}{\partial z^2}=0,$$

(16)

$$x_{\mathrm{Ar},pl}=1-x_{{\mathrm{H}}_2,pl}.$$

(17)

Boundary conditions, which are applied at inlet of the powder layer, at z = h_cat + h_dm:

$${-{\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},pl}^{eff}\frac{\partial x_{{\mathrm{H}}_2,pl}}{\partial z}|}_{z=h_{cat}+h_{dm}}=Q_{dm}P\left(\sqrt{x_{{\mathrm{H}}_2,dm}}-\sqrt{x_{{\mathrm{H}}_2,pl}}\right).$$

(18)

The analytically calculated value for $x_{{\mathrm{H}}_2,dm}$ (Equation (15)) is used in Equation (18). The linkage between the powder layer and next, the intermediate layer, requires equality of the fluxes as a boundary condition at z = h_cat + h_dm + h_pl:

$$-{\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},pl}^{eff}\frac{\partial x_{i,pl}}{\partial z}={\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},il}^{eff}\frac{\partial x_{i,il}}{\partial z}.$$

(19)

Here, in Equation (19), $D_{{\mathrm{H}}_2-\mathrm{Ar},pl}^{eff}$ and $D_{{\mathrm{H}}_2-\mathrm{Ar},il}^{eff}$ are effective binary diffusion coefficients in the powder and intermediate layers, respectively, accounting for both free-molecular $D_{{\mathrm{H}}_2-\mathrm{Ar}}$ and Knudsen diffusion $D_{{\mathrm{H}}_2,pl}^{kn}$, $D_{\mathrm{Ar},pl}^{kn}$ via Wilke–Bosanquet’s-type relation (a harmonic mean approximation), which is written as being dependent on the value of obstruction factor $\frac{{\epsilon }_{pl}}{{\tau }_{pl}}$ in the following form:

$$D_{{\mathrm{H}}_2-\mathrm{Ar},pl}^{eff}=D_{\mathrm{Ar}-{\mathrm{H}}_2,pl}^{eff}=\frac{{\epsilon }_{pl}}{{\tau }_{pl}}·\frac{1}{2}\left(\frac{1}{1/D_{{\mathrm{H}}_2,pl}^{kn}+1/D_{{\mathrm{H}}_2-\mathrm{Ar}}}+\frac{1}{1/D_{\mathrm{Ar},pl}^{kn}+1/D_{{\mathrm{H}}_2-\mathrm{Ar}}}\right).$$

(20)

Diffusion coefficient for binary mixture $D_{{\mathrm{H}}_2-\mathrm{Ar}}$ is calculated as described above [58]. Knudsen diffusivity is a property of a single component and depends inversely on the molecular dimensions and on mean pore diameter [59]:

$$D_{i,pl}^{kn}=\frac{d_{pore,pl}}{3}{\left(\frac{8RT}{\pi M_i}\right)}^{1/2}=48.5d_{pore,pl}\sqrt{\frac{T}{M_i}}.$$

(21)

The void fraction, porosity, ${\epsilon }_{pl}$, restricts the cross-sectional area available for transport in a porous medium. The tortuosity, ${\tau }_{pl}$, accounts for the increase in path length which the molecules must follow. In the current model, mean pore size d_pore,pl (hydraulic pore diameter in Table 2) is applied to be in the interdependent relation between void fraction and volumetric surface area per unit volume, S_V,pl:

$$d_{pore,pl}=\frac{4{\epsilon }_{pl}}{S_{V,pl}}.$$

(22)

The geometric nature of the ${\tau }_{pl}$ factor depends on the void space topology. A linear function being obtained for five structures with void fractions between 0.10 and 0.42 is defined as follows [60,61]:

$${\tau }_{pl}=5.0-4{\epsilon }_{pl}.$$

(23)

Intermediate Layer

Equations and parameters that describe the hydrogen gas transport inside the meso/macroporous intermediate layer of 0.6 mm thick (Table 2) are similar to those for the powder layer:

$${\epsilon }_{il}{\rho }_{tot}\frac{\partial x_{{\mathrm{H}}_2,il}}{\partial t}={\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},il}^{eff}\frac{{\partial }^2{x}_{{\mathrm{H}}_2,il}}{\partial z^2},$$

(24)

$$x_{{\mathrm{H}}_2,il}=1-x_{\mathrm{Ar},il}.$$

(25)

Boundary conditions are as follows:

$$atz=h_{cat}+h_{dm}+h_{pl}-{\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},il}^{eff}\frac{\partial x_{{\mathrm{H}}_2,il}}{\partial z}={\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},pl}^{eff}\frac{\partial x_{{\mathrm{H}}_2,pl}}{\partial z},$$

(26)

$$atz=h_{cat}+h_{dm}+h_{pl}+h_{il}-{\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},il}^{eff}\frac{\partial x_{{\mathrm{H}}_2,il}}{\partial z}={\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},foam}^{eff}\frac{\partial x_{{\mathrm{H}}_2,foam}}{\partial z}.$$

(27)

Foam Substrate I Layer

The manufactured open-cell foams are typically characterized by their pore size being expressed in terms of pores per linear inch (PPI). The pore density in PPI is categorized as low (<20 PPI), medium (<50 PPI) or fine pore (100 PPI). In the quantitative simulation of the porous solid–fluid systems, the basic morphological and geometrical parameters of foam structures, namely cell and window diameter, strut (a bar-like part of the cell network to resist the compression) diameter as the characteristic lengths, surface area-to-volume ratio and porosity are greatly important [62]. The cell diameter d_cell is the most reliable quantity to be measured with simple optical techniques. Correlation for the cell diameter ($0.6<\epsilon <0.95$) with the pore diameter and the open porosity is given in [63].

The hydrogen gas transport inside foam layer of the h_foam = 4.5 mm thickness (Table 2) is also governed by Fick’s diffusion:

$${\epsilon }_{foam}{\rho }_{tot}\frac{\partial x_{{\mathrm{H}}_2,foam}}{\partial t}={\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},foam}^{eff}\frac{{\partial }^2{x}_{{\mathrm{H}}_2,foam}}{\partial z^2},$$

(28)

$$x_{\mathrm{Ar},foam}=1-x_{{\mathrm{H}}_2,foam},$$

(29)

$$D_{{\mathrm{H}}_2-\mathrm{Ar},foam}^{eff}={\epsilon }_{foam}\left(\frac{D_{{\mathrm{H}}_2-\mathrm{Ar}}}{{\tau }_{foam}}+0.5d_{p,foam}{u}_{foam}\right),$$

(30)

$$d_{p,foam}=\frac{d_{cell}}{\left(3.7033-2.5516{\epsilon }_{foam}+0.7054{\epsilon }_{foam}^2\right)},$$

(31)

$${\tau }_{foam}=1+\frac{4.867\left[1-0.971{\left(1-{\epsilon }_{foam}\right)}^{0.5}\right]}{4{\epsilon }_{foam}{\left(1-{\epsilon }_{foam}\right)}^{0.5}}\left(1-{\epsilon }_{foam}\right).$$

(32)

Here, the effective pore diffusivity $D_{{\mathrm{H}}_2-\mathrm{Ar},foam}^{eff}$ is a function of the structural parameters of the heterogeneous porous medium Equation (30), viz., d_p,foam pore diameter, the average diameter of the windows which connect the cells, Equation (31); tortuosity ${\tau }_{foam}$, Equation (32); and u_foam, the interstitial (average pore) velocity of gas (argon and hydrogen) [64].

The molar fluxes boundary conditions are written as follows:

$$atz=h_{cat}+h_{dm}+h_{pl}+h_{foam}-{\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},il}^{eff}\frac{\partial x_{{\mathrm{H}}_2,il}}{\partial z}={\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},foam}^{eff}\frac{\partial x_{{\mathrm{H}}_2,foam}}{\partial z},$$

(33)

$$atz=h_{cat}+h_{dm}+h_{pl}+h_{il}+h_{foam}-{\rho }_{tot}{D}_{{\mathrm{H}}_2-\mathrm{Ar},foam}^{eff}\frac{\partial x_{{\mathrm{H}}_2,foam}}{\partial z}={\beta }_{foam}{\rho }_{tot}\left(x_{{\mathrm{H}}_{2}foam}-x_{{\mathrm{H}}_2,sw}\right).$$

(34)

In the present simulations, the effective sweep gas–foam surface mass transfer coefficient ${\beta }_{foam}$ is expressed in the form of Sherwood numbers and accounted for by the following specific mass transfer correlations for foam configurations [63]:

$$S{h}_{d_{s,avg}}={\epsilon }_{foam}^{-2}\left(0.566R{e}_{d_{s,avg}}^{0.33}+0.039R{e}_{d_{s,avg}}^{0.8}\right)S{c}_{foam}^{1/3}$$

(35)

The characteristic length in the dimensional Reynolds Re_s,avg, Sherwood $S{h}_{d_{s,avg}}$ and the Schmidt Sc_foam numbers is the average cylindrical strut size [65].

$$d_{s,avg}\cong \frac{2.85·\left(1-{\epsilon }_{foam}\right)}{S_{V,foam}},$$

(36)

$$R{e}_{s,avg}=\frac{{\rho }_{foam}·u_{foam}·d_{s,avg}}{{\mu }_{foam}},$$

(37)

$$S{h}_{d_{s,avg}}=\frac{{\beta }_{foam}·d_{s,avg}}{D_{{\mathrm{H}}_2-\mathrm{Ar},foam}^{eff}},$$

(38)

$$S{c}_{foam}=\frac{{\mu }_{foam}}{{\rho }_{foam}·D_{{\mathrm{H}}_2-\mathrm{Ar},foam}^{eff}},$$

(39)

$$S_{V,foam}=\left(\frac{2\sqrt{3\pi }}{d_{cell}}\right){\left(1-{\epsilon }_{foam}\right)}^{0.5}.$$

(40)

Sweep Compartment

The unknown variable in Equation (34) is $x_{H_2,sw}$ molar fraction of hydrogen in the stream leaving the sweep compartment. At steady state, assuming a perfectly mixed gas, mass balance equations for hydrogen Equation (40) and argon Equation (41) in the volume V_sw of sweep compartment at z|₊ = h_cat + h_dm + h_pl + h_il + h_foam are written down as follows:

$$\begin{gathered} \frac{V_{sw}{\rho }_{sw}·}{A_m}\frac{d{x}_{{\mathrm{H}}_2,sw}}{dt}={\beta }_{foam}{\rho }_{tot}\left(x_{{\mathrm{H}}_2,foam}-x_{{\mathrm{H}}_2,sw}\right)+\frac{{\rho }_{tot}}{A_m}\left(G_{sw}^0·x_{{\mathrm{H}}_2,sw}^0-G_{sw}·x_{{\mathrm{H}}_2,sw}\right), \\ \frac{V_{sw}{\rho }_{sw}·}{A_m}\frac{d{x}_{{\mathrm{H}}_2,sw}}{dt}={\beta }_{foam}{\rho }_{tot}\left(x_{{\mathrm{H}}_2,foam}-x_{{\mathrm{H}}_2,sw}\right)-\frac{{\rho }_{tot}{G}_{sw}}{A_m}·x_{{\mathrm{H}}_2,sw}. \end{gathered}$$

(41)

$$\begin{gathered} \frac{V_{sw}{\rho }_{tot}}{A_m}·\frac{d{x}_{\mathrm{Ar},sw}}{dt}=\frac{{\rho }_{tot}{G}_{sw}^0}{A_m}{x}_{\mathrm{Ar},sw}^0-\frac{{\rho }_{tot}{G}_{sw}}{A_m}·x_{\mathrm{Ar},sw}, \\ \frac{V_{sw}{\rho }_{tot}}{A_m}·\frac{d{x}_{\mathrm{Ar},sw}}{dt}=\frac{{\rho }_{tot}{G}_{sw}^0}{A_m}-\frac{{\rho }_{tot}{G}_{sw}}{A_m}{x}_{\mathrm{Ar},sw}, \end{gathered}$$

(42)

$$x_{{\mathrm{H}}_2,sw}^0=0;G_{sw}^0=G_{\mathrm{Ar},sw}^0;x_{\mathrm{Ar},sw}^0=1$$

(43)

$$G_{sw}=G_{sw}^0+{\beta }_{foam}{A}_m\left({x_{{\mathrm{H}}_2,foam}|}_{z=z_{h_{foam}}^{-}}-x_{{\mathrm{H}}_2,sw}\right).$$

(44)

The molar flow rate for the species H₂ and Ar in the sweep gas compartment is related to the volumetric flow change using the ideal gas law. Summing Equations (40) and (41) and taking into account that $x_{{\mathrm{H}}_2}+x_{\mathrm{Ar}}=1$, the volumetric gas flow accounting for a change in the total molar flow rate of the sweep gas at a given temperature Equation (43) can be obtained. Moreover, the explicit formula for the hydrogen molar fraction in the stream leaving the sweep compartment is yielded at steady state:

$$x_{{\mathrm{H}}_2,sw}=\frac{\left({\beta }_{foam}+\frac{G_{sw}^0}{A_m}+{\beta }_{foam}\right)-\sqrt{{\left({\beta }_{foam}+\frac{G_{sw}^0}{A_m}+{\beta }_{foam}·x_{{\mathrm{H}}_2,foam}\right)}^2-4{\beta }_{foam}{x}_{{\mathrm{H}}_2,foam}}}{2{\beta }_{foam}}.$$

(45)

The scaled reactor model equations for catalyst–membrane assembly layers and the sweep compartment with the boundary conditions have been implemented in COMSOL Multiphysics 3.5 software to execute the numerical model and solved simultaneously in the five domains along the full length. The finite element method for numerical solutions of differential equations being used employs a uniform (fine) mesh and error control in the domains. A numerical problem arises from the initial values in the feed, which generates a division by zero in the reaction rate equations. This problem has been sorted out by using very small values for the mole fractions of the generated species at the inlet. The solver, an implicit time-stepping scheme, is well suited for solving stiff and nonstiff, nonlinear boundary value problems. At the end of the solving process, species concentrations and fluxes are known at each axial location within all reactor’s domains.

2.3. Performance Parameters of the Reactor

The resistance concept proposed by Henis and Tripodi [66] is widely used in research to have insight and understanding about which of the several layers of the membrane assembly is controlling the total hydrogen flux. Irrespective of transport mechanisms, a total transport resistance consists of the layers’ resistance in series and also includes resistance by gas film on both sides of the membrane assembly [67,68,69].

With respect to modeling, the overall permeation hydrogen flux should comprehend the mass transport across dense, powder, intermediate and foam layers. The apparent resistance for permeation through a layer (n = 1−4) is influenced not only by a partial pressures difference $\mathsf{\Delta }{p}_{{\mathrm{H}}_2}$ across each layer but also by the total flux and is defined as

$$R_{S,n}=\frac{\mathsf{\Delta }{p}_{{\mathrm{H}}_2,n}}{J_{{\mathrm{H}}_2}RT}=\frac{P\mathsf{\Delta }{x}_{{\mathrm{H}}_2,n}}{J_{{\mathrm{H}}_2}RT}=\frac{{\rho }_{tot}\mathsf{\Delta }{x}_{{\mathrm{H}}_2,n}}{J_{{\mathrm{H}}_2}}.$$

(46)

Resistance to mass transfer by gas film on both sides of the membrane assembly (feed- and sweep-side compartments) is calculated as the inverse to mass transfer coefficients Equation (46). The overall hydrogen transport resistance is calculated from Equation (47).

$$R_{S,f}=\frac{1}{{\beta }_f}\mathrm{and}{R}_{S,sw}=\frac{1}{{\beta }_{foam}}$$

(47)

$$R_S=R_{S,f}+R_{S,dm}+R_{S,pl}+R_{S,il}+R_{S,foam}+R_{S,sw}.$$

(48)

Performance characteristics describing the overall efficiency of the integrated reaction–separation process in the lab-scale CMR are expressed in terms of hydrogen recovery and yield of ethanol conversion. A hydrogen recovery factor is calculated as a molar ratio of hydrogen permeated through the membrane to hydrogen produced from the ethanol steam reforming reaction in the feed-side Equation (48).

$${\mathrm{RecoveryH}}_2=\frac{{\mathrm{H}}_2\mathrm{permeated}\left(\mathrm{mol}{\mathrm{s}}^{-1}\right)}{\left({\mathrm{H}}_2\mathrm{permeated}+{\mathrm{H}}_2\mathrm{retentate}\right)\left(\mathrm{mol}{\mathrm{s}}^{-1}\right)}·100\%,$$

(49)

$$\mathrm{Yield}{\mathrm{H}}_2=\frac{\left({\mathrm{H}}_2\mathrm{permeated}+{\mathrm{H}}_2\mathrm{retentate}\right)\left(\mathrm{mol}{\mathrm{s}}^{-1}\right)}{6·\mathrm{ethanol}\left(\mathrm{mol}{\mathrm{s}}^{-1}\right)}·100\%.$$

(50)

In regard to evaluation of the efficiency of catalytic process in the CMR, the percent yield is calculated (Equation (49)) as the extent to which the reaction theoretical yield is achieved. Both reactants (water and ethanol) contain hydrogen atoms and contribute to the hydrogen yield. The theoretical yield of the reaction is found to be six moles hydrogen per mole of reacted ethanol, once the ethanol steam reforming reaction is complete.

3. Results and Discussion

3.1. Verification

The most important question in any simulation is the reliability and validity of an assessment tool and whether we can get more out of our model by tuning the model parameters. Since the quality and significance of experiments play a crucial role, the applicability of the one-dimensional model formulation and computational simulations are verified and validated by direct comparison of model results with existing measurements of interest. Certainly, there are some shortcomings of using only one validation data set. Aiming to improve the usefulness and predicting performance of the model, training on the data derived from the reactor operating at quite different conditions is applied. Two experimental datasets (Table 3, No.1 and Nos.2–5) were used to thoroughly evaluate and quantify all necessary parameters for the developed model and the simulation procedure. In Figure 4, the results of the simulation for the hydrogen flux and concentrations of species are compared to the known values in the experimental study being carried out under specified operating conditions, Table 3, No.1. Good agreements with the experimental data on hydrogen permeation flux and on species concentration in both feed and sweep sides are observed. In particular, by considering effect of the temperature shown in Figure 4, the average error between experimental and modeling data for hydrogen fluxes is about 2%, while a highest deviation of about 6% is observed for hydrogen concentration in the retentate gas at 900 °C. As a result, model parameter values produce reasonable fits to the data regarding the considered range of experimental study.

Comparison between experimental data and numerical modeling in the case of high flowrates of both the feed and sweeping gas is shown in Figure 5. Experiments Nos.2–5, Table 3, were carried out in order to evaluate the influence of ethanol concentration [34]. It can be seen that the model provides more accurate quantitative predictions for the reactor performance at 6.6% EtOH in the feed, with a steam-to-carbon ratio of 2.55 (No.5, Table 3). This is most probably due to the failure of the applied kinetic model in predicting the kinetic behavior of ethanol steam reforming reaction at higher steam-to-carbon ratios. Effect of the temperature on species concentrations in the retentate and sweep streams is shown in Figure 6a. In particular, for the hydrogen concentrations in retentate and in sweep gas streams, minimum values of an average error are estimated to be about 3–5%.

In Figure 6b we compare catalytic and separation functions for this reactor operating at different flowrates (Nos.1 and 5, Table 3). The definitions of hydrogen recovery (Equation (48)) and yield of hydrogen (Equation (49)) have been used for comparing the reactor performance. It can be seen that, at temperature of 700 °C, the characteristics for both operating regimes, so far considered, are nearly identical. An increase in the temperature favors the overall efficiency of the integrated reaction–separation process in the reactor. However, for the case study No.5, increase in the temperature from 700 °C up to 900 °C results in improving the yield of hydrogen by nearly 40%, while hydrogen recovery carried out with the membrane becomes lower than that for the comparative case study No.1. In the case of higher flowrates, the yield of hydrogen increases, while permeate activity decreases. It is evident that, despite the fact that the total hydrogen flux is nearly the same in both cases (see Figure 4a or Figure 5b), the lab-scale catalytic membrane reactor operating at lower flowrates (No.1) shows better performance in terms of hydrogen recovery as compared with that at high flow rates (No.5).

Therefore, the developed one-dimensional isothermal reaction–transport model for the catalytic membrane reactor is reliable in description of the experimental research. Direct validation of the model formulation and computational simulations by using the data recorded in the experiments with the disk-shaped CMR shows that an acceptable agreement between experiment and simulation is achieved with the error less than 6%. This ensures applicability of the computational model to examine the effects of the operating conditions and membrane structural properties on the performance of the catalytic membrane reactor.

3.2. Resistance to the Hydrogen Mass Transfer

The concentration polarization resistance due to boundary layer effects and internal concentration polarization within the composite multilayer membrane structure is an inherent functional property of any membrane module. The phenomenon of polarization arises due to nonuniform species transport. The concentration distribution along the catalyst–membrane assembly module affects the overall hydrogen permeation flux, which, in turn, depends on mass transport parameters: reaction kinetics, Peclet number and effective diffusivity in the constituent layers [70,71]. Understanding all these phenomena is important when studying factors that direct the hydrogen removal.

It was shown [31] that thermodynamic imposes constraints concerning the hydrogen concentration at temperatures higher than 650 °C. However, under the operating conditions, of which Peclet numbers are in the range of 0.2–0.8, transport of reactants to the interface of catalytic particles on which the chemical reaction takes place is limited by diffusion that mitigates the thermodynamic driving force, and the catalytic process is far from equilibrium. The validated numerical model with the assumptions made, in terms of the reaction pathway, kinetic expressions and parameters, as well as the solution of the partial differential equations using the COMSOL Multiphysics^® package, have been used to elucidate the concentration distribution and to study the reactor–catalyst–separation system in detail. A layer with the highest concentration gradient would provide the major mass transfer resistance to the overall permeation process.

Typical concentration profiles as those calculated for the case No.1, Table 3 (see Figure 7a or Figure 8a) evolve throughout the catalyst–membrane assembly. According to the reaction network (Equations (1) and (2)), ethanol fully decomposes at the inlet part of catalytic layer and generates intermediate products. Steam is consumed in the reforming reactions (Equations (2)–(4)), and its molar fraction continuously decreases with the axial length. Along with CO and CO₂, methane is also formed as a byproduct of the catalytic reforming reaction, up to a few percent.

The hydrogen concentration increases up to maximum value due to the high ethanol reforming reaction rate. After its maximum value, the concentration gradually decreases down toward the membrane surface so that concentration gradient adjacent to the membrane is established. In addition, hydrogen removal from the reaction zone leads to an increase in CO mole fraction and lowers CO₂ and CH₄ fractions, which corresponds to the higher amount of syngas produced. As to the experimental case No.5, Table 3 (Pe = 0.805 at 700 °C), the concentration profile is becoming flatter due to increased convective velocity (plots are not being displayed in the paper for brevity). Concentration gradients on both sides of the membrane indicate that mass transfer resistance in the gas phase, indeed, contribute to the overall mass transfer resistance.

The polarization resistance curve in the case under consideration (Table 3, No.1) is shown in Figure 7b. High temperature is beneficial to syngas formation (Figure 8a). When the operating temperature increases from 700 up to 900 °C, maximum concentration of the formed hydrogen increases by ≈14%. The dependence of hydrogen flux on the temperature shows the same trend and increases by ≈22%, as can be seen in Figure 4a. At the same time, the overall resistance across the catalyst–membrane module decreases by ≈25% as the temperature increases (Figure 8b or Figure 9a). The values are, at least, in the same order of magnitude. This shows that, all these phenomena considered, the reaction and permeation rates, and also the resistance to hydrogen transport, are important in improving the overall efficiency of the integrated reaction–separation process. The derived model allows us to describe resistance to the hydrogen mass transfer as a function of operating conditions and the structural properties of each single layer. For operational conditions with lower flow rates (Table 3, No.1), the contribution of the permselective dense layer to the overall resistance is about 9%.

Plotting the calculation results on the impact of operating temperature on the resistance at various feed and sweep gas flow rates produces Figure 9. The model predictions suggest that the overall resistance is affected by temperature and is consistently higher at low temperature in the case of high flow rates, due to the lower permeance of the permselective dense layer at these operating conditions. The contribution of the permselective layer resistance drops from nearly 40% at 700 °C to 20% at 900 °C. A conjecture on the explanation for variable resistance of the dense layer is the nonideal membrane behavior at high operating flow rates, in which external phenomena may also play a role during permeation [72,73]. This indicates that the real driving force does not actually obey the Sieverts-type empirical law at the high operating flow rate.

The role of the support layers is demonstrated in Figure 9b. It is interesting that the hierarchically structured asymmetric support, which is required for desirable industrial applications of membranes, offers a similar extent of resistance to hydrogen transport under these different operating conditions (Figure 9b). The results obtained clearly evidence that the support contributes up to 70% to the overall resistance across the membrane module. As a consequence, in order to improve the separation efficiency of the CMR, a thorough evaluation of morphological and structural properties of the asymmetric support is a fundamental task, since it controls internal concentration polarization within the composite multilayer membrane structure.

3.3. Effect of the Asymmetric Support

Changes in a layer morphology could, in principle, change the mass transport parameters. The following subsection reports the results of numerical research on the effect of structural parameters on the performance enhancement due to creating higher driving force for hydrogen permeation. Following of the fact that the more profound negative impact of the concentration polarization in the support is observed for the low flowrate, simulation results for the operating conditions of the set No.1, Table 3, have been considered in detail. Based on the chosen structural parameters of the constituent layers (Table 2), the major mass transfer resistance to the overall permeation process is found to be provided by the powder layer (see Figure 8b or Figure 9b), that is, the transport in support is mainly controlled by the characteristics of the layer immediately followed by the permselective dense layer.

The plots in Figure 10 elucidate the effect of the structural parameters of the powder layer on the process parameters, such as hydrogen concentrations and permeation flux. It is possible to observe that a powder/intermediate layer ratio is one of the most important factors affecting the process performance. The situation in which the ratio is about 0.66 corresponds to the actual membrane module. The influence of the powder/intermediate layer ratio is appraisable up to the ratio of about 2. Pore dimensions have a lower contribution to the overall permeation if compared to that for porosity of the powder layer. Changing the thickness and porosity of the powder layer mainly affects the concentration in the sweep flow (Figure 10a) and the amount of hydrogen permeated (Figure 10b) through a given gas-permeable supporting structure. Such behavior is attributed to higher catalytic reaction rates compared to the hydrogen removal rate through the membrane module.

In contrast, the curves in Figure 11 show a weak dependency of the hydrogen flux on the intermediate layer porosity. Thus, for the reference case (a powder/intermediate layer ratio of 0.66), the increase of porosity from 0.3 up to 0.5 results in the increase of hydrogen flux by about 2.6 %.

As for the foam layer that contributes significantly (25–26%) to the overall resistance across the membrane module, the influence of the thickness and porosity on the process performance becomes remarkable (Figure 12). Simulation illustrates the key role of the foam thickness for the hydrogen flux (Figure 12a). Results referring to the hydrogen concentration and flux versus foam porosity are shown in Figure 12b. The results are presented at three different pore densities, namely 20, 30 and 40 PPI. With increasing PPI, the solid phase content in the foam structure per unit volume increases, resulting in an increase of the surface area per unit bulk volume of the foams, and, therefore, in improving the mass transfer inside the foam structure. According to the graphs in Figure 12b, both concentration and hydrogen removal do not change significantly. As demonstrated by the simulations, at ${\epsilon }_{foam}=0.75$ the hydrogen flux differs by about 1.2% between the 20 PPI and 40 PPI foams.

The analysis carried out above allows us to understand the influence of such factors as structural parameters of the constituent layers of the membrane module on the separation option of the reactor and provides the general trends of improving performance of the asymmetric membrane module with the variation of the morphology of the support layers. It is evident that a more advantageous structure can be predicted by the developed model equations, thereby creating a higher driving force for the hydrogen permeation. This consideration is of a great practical importance.

4. Conclusions

To obtain a better fundamental understanding of interrelated catalytic and permeation phenomena in the catalytic membrane reactor for ethanol steam reforming, numerical experiments have been performed. One-dimensional isothermal reaction–transport model for the constituent layers of the catalyst–asymmetric membrane assembly, together with a Sieverts’ equation for the functional dense layer, as well as with boundary layer effects, was developed and implemented in COMSOL Multiphysics. The mathematical model was experimentally validated with performance data from a lab-scale reactor with disk-shaped membrane operating under different flow rate conditions. Good agreements with the experimental data with reasonable parameters values were provided. The simulations demonstrated the model’s capabilities to efficiently capture the concentration distribution and hydrogen mass transfer resistance throughout the catalyst–asymmetric membrane assembly. The performance of the CMR, in terms of efficiency of catalytic process and hydrogen recovery, was also studied.

The numerical results comprise the following major findings: firstly, the asymmetric support contributes up to 70% to the overall resistance across the membrane module; secondly, the transport through the support is mainly controlled by the characteristics of the layer immediately followed by the permselective dense layer. In summary, a strong impact of the structural parameters of constituent layers of composite membranes on the performance enhancement due to a higher driving force for hydrogen permeation has been quantified.

The developed model can be used as a basis for optimal design of asymmetric porous supports with a purpose to minimize mass transfer limitations and performance losses introduced by the gas-permeable support that is very important for different technological applications of membranes. Through numerical analysis of catalysis–permeation process and simulation study, the integrated design and process control aiming at the maximum driving force can be performed.