2.2. Computational Model

The membrane permeator cell was modeled using CFD by considering its real dimensions: internal diameter of 30 mm, 30 mm in length, and a total volume of 21,200 mm³. As the geometry presented a symmetry plane, only half of the device was considered in the model. The tubular composite membrane had inner and outer diameters of around 11 mm and 13 mm, respectively, dividing the model geometry into three different regions, as shown in Figure 1. Depending on the flow direction through the membrane module, two different configurations were modeled. In the in–out configuration (Figure 1a), the gas was fed into the reactor through the inner region of the tubular membrane, collecting the retentate stream from the opposite side along the z-direction and the permeate stream in the annular zone. Thus, the H₂ first reached the porous media (modified PSS with intermediate layers) and then the fully dense Pd film. A contrary permeation scheme was considered for the out–in configuration (Figure 1b), where the feed stream was introduced in the cell through the outer region. Therefore, the fed H₂ was in direct contact with the outer surface of the Pd film and the permeate had to go through this layer and then the porous media before reaching the lumen bulk gas phase.

Schematic representation of the membrane permeator working in two different configurations: (a) in–out and (b) out–in. PSS: porous stainless steel.

A conformal mesh was established for the geometry interfaces to reduce the computational cost. The mesh was carefully refined close to the PSS and dense palladium region. This mesh slightly varied for the two different configurations considered in the present study, with 74,935 and 78,694 elements for the in–out and out–in alternatives, respectively. This number of cells was high enough to give mesh-independent results. An example of the established mesh for the out–in configuration is presented in Figure 2, with this being analogous to the other configuration, namely, in–out.

Illustrative mesh details for the membrane permeator for the out–in configuration.

The porous zone, corresponding to the PSS support, was implemented using the Darcy–Forcheimer model. This formulation provides accurate predictions for relatively low-pressure drops by adding a sink term (S_i) to the fluid flow equations [32]. Its implementation in the CFD platform allows for calculating the variation of the velocity due to the effect of the membrane in each spatial direction:

where Si is the sink term, v the velocity, μ the viscosity, ρ the density, and RV and RI are, respectively, the viscous and the inertial vectors of resistance coefficients (m⁻¹) for the ith spatial direction.

Viscous and inertial resistance coefficients for each direction were calculated using:

where ε is the stainless steel’s porosity, Dp is the particle diameter, and Φ is the particle sphericity.

According to commercially available data, the porosity of commercial 0.1 μm PSS supports is around 20%. This value was confirmed using SEM image segmentation and analysis with Digital Micrograph^® (Gatan, Pleasanton, CA, USA), obtaining values of ε = 21.48% and D_p = 21 µm, as previously reported elsewhere [30]. Finally, the PSS particles’ sphericity was set to Φ = 0.9. As shown in Figure 3, the membrane was constituted by the agglomeration and syntherization of large steel particles with variable sizes and morphologies. The average particle diameter and porosity were directly calculated from the segmentation of SEM images, while an average sphericity value of Φ = 0.9 was estimated by matching experimental data of nitrogen flow through the membrane to the model predictions.

Morphology of the considered composite Pd-membrane for modeling: surface and cross-sectional SEM images at diverse stages of the membrane fabrication. BSE: back-scattering electrons, SE: secondary electrons, SS: stainless steel.

With these parameters, viscous and inertial resistance coefficients of 2.61 × 10¹³ m⁻¹ and 3.15 × 10⁶ m⁻¹ were calculated, respectively, for each spatial direction. A factor of 10³ was applied in the axial direction as the permeate flux crossed the porous steel in the driving force direction, which was radial to the geometry. These parameters were validated using permeation experiments of PSS supports before incorporating the Pd film with pure nitrogen at two different feed pressure values.

Assuming that the Pd film was placed on the outer surface of the porous support (as previously described in Section 2.1) and the feed stream was introduced by the lumen of the tubular support (in–out configuration), the complete H₂ permeation mechanism involved the following fundamental steps [33]: (i) external mass transfer in the gas phase, (ii) transport through the porous media, (iii) hydrogen adsorption on the internal metal surface, (iv) hydrogen dissociation into atomic hydrogen, (v) diffusion through the bulk palladium, (vi) hydrogen recombination on the outer side of the palladium film, (vii) hydrogen desorption from the metal surface, and finally, (viii) external mass transfer in the gas phase of the permeate side. For the out–in configuration, in which the gas was fed to the annular region, similar steps could be distinguished, although in a different sequence.

The H₂ solution and diffusion through the bulk palladium is typically the controlling step when the Pd thickness is thicker than a few micrometers [34]. In these conditions, the total permeate H₂ flux (JH2) can be described using the well-known Sieverts’ law [35], in which it is defined as a function of the membrane´s permeance K (mol·s⁻¹·m⁻²·Pa^−0.5) and the difference between the permeate (PH2,per) and retentate (PH2,ret) pressures to the exponent n = 0.5 acting as a driving force between both sides of the Pd film:

In this study, the Pd film was modeled as a theoretically zero-thickness wall. The introduction of source and sink terms (Si) in the continuity and species transport equations simulated hydrogen’s disappearance from one side and its appearance on the palladium wall’s opposite side. Source term formulation was considered for each cell near this layer. Thus, Sieverts’ law could be reformulated to transform the surface flow into a local volumetric source–sink pair (Si), as follows:

where Ac and VC are the cell area and volume, respectively, and Yi is the species mole fraction. Each parameter included in Equation (5) was taken from its respective cell near the Pd film for each simulation iteration. As the hydrogen concentration in the permeation side was always maintained, in theory, equal to 1 (assuming perfect selectivity to hydrogen), Equation (5) could be rewritten as follows:

The user-defined functions (UDFs) for the source and sink terms can be found as Supplementary Information. The experimental K values obtained for each considered H₂ feed molar fraction ranging from 1 to 0.5 were used to evaluate the model performance at diverse pressure driving forces. Despite being considered a zero-thickness wall in the model geometry (for the sake of simplicity), the thickness of the Pd membrane, and its impact on the H₂ permeation, was implicitly considered by the inclusion of the specific K values in the calculations. The same approach was followed with the CeO₂ intermediate layer, whose effect was directly included in the experimental K values.

Here, it should be noted that most of the Pd composite membranes prepared using ELP-PP exhibited a deviation from the ideal Sieverts’ law, meaning that it was necessary to consider a certain value for the y-intercept instead of directly fitting to the origin (0,0). This deviation does not yet have a clear physical interpretation, although an additional resistance to the H₂ transport due to the partial penetration of the palladium in some of the PSS pores was proposed in previous studies [18].

The y-intercept and K values for each configuration used in the CFD model are presented in the result section.

Fluid dynamics calculations have been carried out through the combination of the continuity equation (Equation (7)) and the classical Navier–Stokes equation (Equation (8)):

where ρ, v¯, P, τ=, and g→ are the fluid density, velocity vector, pressure, viscous stress tensor, and gravitational acceleration, respectively.

The standard k–ϵ model was used to model the turbulent flow. This two-equation turbulence model considered the turbulent length and time scale by solving two separate transport equations for the turbulent kinetic energy (k) and its dissipation rate (ϵ). This standard model assumes an isotropic scalar turbulent viscosity, treats the flow as completely turbulent, and neglects the molecular viscosity. Therefore, it is only valid for fully turbulent fluids. The energy transport, the dissipation rate, and the turbulent viscosity are described by Equations (9)–(11), respectively [36]:

In these expressions, σk and σϵ are the turbulent Prandtl numbers for k and ϵ, respectively; C1 and C2 are the model empiric constants that can be found elsewhere [36]; Sk and Sϵ are the terms for user-defined contributions.

Simulations with the turbulent model were compared with those using the laminar one. Complete correspondence between both models was observed except in a few cases, where very slight differences were found.

Hydrodynamics calculations were carried out by considering a three-dimensional, steady-state, and turbulent flow after solving both the continuity and classical Navier–Stokes equations. The fluid was assumed to be Newtonian, compressible, and isothermal (648 K). The thermophysical properties of the mixture were determined using the Peng–Robinson equation of state. The hydrogen diffusion coefficient was calculated using a correlation with diffusion coefficients values taken from the Aspen Plus^® V11 software (Aspentech, Bedford, MA, USA) at similar operating conditions to the experiments (pressure and temperature) as follows:

where DH2 is the hydrogen diffusion coefficient (m²/s) and P is the pressure (Pa).

The viscosity of the mixture at the experimental temperature was calculated as a function of the hydrogen molar fraction according to:

where μ is the dynamic viscosity (Pa·s) and Y is the hydrogen molar fraction.

The inlet feed was set with a direction normal to the boundary, where the values ranged from 1.25 × 10⁻⁴ g/s to 5.57 × 10⁻³ g/s as a function of the hydrogen concentration in the mixture. The retentate pressure was set at the experiment feed pressure in both configurations. On the other hand, the permeate was fixed at atmospheric conditions. A no-slip boundary condition was imposed at all model walls. The interface between the porous substrate and geometry domains was set as an interior region.

The segregated steady-state solver was used to resolve the governing equations. A second-order upwind discretization scheme was used, except for pressure, for which the standard method was selected. The SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm was chosen for the pressure–velocity coupling. The numerical solution’s convergence was ensured by monitoring the scaled residuals below 10⁻⁶ for the continuity and momentum variables. Additionally, the variables of interest were monitored at different computational domain surfaces as an indicator of convergence (at least 50 iterations without changes). The complete model was solved using ANSYS Fluent software v.17.0 (Canonsburg, PA, USA).