On concentration polarization in fluidized bed membrane reactors

A R T I C L E I N F O


Introduction
Currently, hydrogen is mainly produced on large scale via steam reforming of methane (SMR) [1]. In this process, methane is first reformed with steam (Eq. (1)) in high temperature multi-tubular packed bed reactors. In a second step the carbon monoxide is converted via the water gas shift (WGS) reaction (Eq. (2)) in packed bed reactors. Typically, a two stage WGS is used to take advantage of fast reaction rates at high temperatures (450°C) and higher equilibrium conversions at lower temperatures (200°C). Finally, the hydrogen produced is further purified using pressure swing adsorption (PSA).
Steam methane reforming reaction (SMR): Water gas shift reaction (WGS): The equivalent hydrogen efficiency of the whole process is approximately 80% thanks to steam/electricity export [2]. The heat integration between the different stages becomes more complicated at smaller scales, while heat export cannot be realized in distributed hydrogen production applications. For this reason the system becomes inefficient and uneconomical at smaller scales. The cost of the hydrogen produced at large scale is around 0.2 €/Nm 3 , while it increases up to 0.4-0.5 €/Nm 3 at smaller scales [2].
The efficiency of the hydrogen production via methane reforming can be increased by integrating hydrogen production and separation in a single multifunctional reactor. This can be achieved by using perm- area (m 2 ) c c , 1 2 constants in frictional stress model (  Recovering the hydrogen during the reaction, results in a shift of the equilibrium towards the products, thus allowing achieving much higher conversions at lower temperatures. The equilibrium displacement (Le Chatelier's principle) allows to minimize the reactor volume up to 80% for WGS [3] and maximize the efficiencies, as total conversion can be achieved already at lower temperatures [4]. In literature, both packed bed and fluidized-bed membrane reactor configurations have been proposed for SMR and WGS reactions. The latest developments in the fabrication of ultra-thin membranes with high permeation rates [5], have once more sparked the debate on the inherent bed-to-membrane mass transfer limitations (concentration polarization) in packed bed membrane reactors [6,7].
From an experimental point of view, Hara et al. [8] studied the decline of hydrogen permeation in a packed bed membrane reactor by injecting the reactor with H 2 -Ar and H 2 -CO mixtures. It was found that the reduction in hydrogen permeation was caused by CO poisoning of the Pd based membrane and concentration polarization near the membrane wall. It was concluded that in order to fairly predict the membrane reactor performance, concentration polarization needs to be taken into account.
Mori et al. [9] investigated the influence of concentration polarization on hydrogen production via SMR in a packed bed membrane reactor with a highly permeable membrane. They performed experiments and compared them with a simple model that did not take into account the effect of concentration polarization. By increasing the reactor pressure, they found that the experimental methane conversion was lower than the simulated conversion. This implies that concentration polarization is occurring in the reactor and affects the methane conversion. The presence of concentration polarization was confirmed with experiments with a binary mixture of hydrogen and nitrogen.
Caravella et al. [6] made a model predicting the permeance of hydrogen in a hydrogen-nitrogen mixture including the effect of concentration polarization in an empty annular tube. It was found that the effect of polarization is relevant not only for the very thin membranes (1-5 µm) with high fluxes but also for the thicker ones (100 µm) at certain operating conditions.
In a Computational Fluid Dynamics (CFD) study by Nekhamkina et al. [10], the mass transfer processes in two configurations were studied: an empty reactor with (i) the membrane at the wall, and (ii) an annular cylinder with the membrane as the inner tube. A model was developed to predict the membrane flux considering the effect of concentration polarization. A parameter Γ was defined which represents the ratio of the diffusion to the permeation flux. It was concluded that only when Γ > 6 the effect of concentration polarization can be neglected.
To circumvent the mass transfer limitations typical of empty or packed bed membrane reactor configurations, fluidized bed membrane reactors were suggested, because of their improved heat and mass transfer characteristics. Patil et al. [11] and Gallucci et al. [12] successfully demonstrated this membrane reactor concept for the SMR reaction with relatively low flux membranes. No concentration polarization effects were reported, but the flux of the membrane used was 5-10 times lower than recently available highly permeable membranes.
More recently, Helmi et al. [13] successfully demonstrated the long term (> 900 h) performance of a fluidized bed membrane reactor utilizing very high flux membranes for ultra-pure hydrogen production via WGS. Although the long term stability of this membrane reactor has been confirmed (with CO content in the permeate side < 10 ppm), no independent study has been performed yet on the concentration polarization effects in fluidized bed membrane reactors where a large amount of gas is extracted via the membranes.
This work focuses on the quantification of the extent of concentration polarization in fluidized bed membrane reactors. A simple onedimensional phenomenological model is developed which can capture the effect of concentration polarization in fluidized bed membrane reactors. The Two-Fluid Model, an Euler-Euler model using the Kinetic Theory of Granular Flow to describe the solids phase rheology, is used to estimate the mass transfer boundary layer thickness required by the 1D model. The predictions by both models are compared with results obtained from experiments showing good agreement, confirming that concentration polarization can also prevail in fluidized bed membrane reactors.

1D phenomenological model
In this work a one-dimensional, two-phase flow model for bubbling fluidized beds was used. This model was originally proposed by Kato and Wen [14] for standard fluidized beds without internals, and was later extended to membrane assisted fluidized bed reactors by Deshmukh et al. [15]. In this model the fluidized bed is divided into a number of Continuous Stirred Tank Reactors (CSTR) in series for both the emulsion and bubble phases, while mass transfer limitation from the bulk of the bed to the surface of the membranes is not accounted for. The number of CSTRs determines the degree of gas back mixing inside the reactor (1 CSTR stands for total back mixing (i.e. ideal mixing) and an infinite number of CSTRs represents the plug flow conditions). In this work, 100 CSTR's in series were chosen for both emulsion and bubble phases assuming a negligible gas back mixing in the axial direction, which is appropriate for the considered small lab scale set-up with a relatively high bed aspect ratio.
In a fluidized bed membrane reactor, the extraction of a large amount of gas may induce the formation of a densified zone (a region with higher solids hold-up) around the membranes [16,17]. This is caused by the high flux through the membrane relative to the fluidization velocity, resulting in a drag force of the particles towards the membranes, which is more pronounced for smaller particles [16]. In the densified zone the mixing is worse than in the bulk of the bed because the particles have much lower velocities (eventually behaving like a slowly moving bed). Because of the good mixing properties of fluidized beds, it is expected that the concentration gradient mainly resides in the densified zone around the membrane, and not in the bulk of the fluidized bed. However, this assumption needs to be validated and implemented in the phenomenological model.
Several theories exist to describe the mass transfer from the bulk of the bed to the surface of the membranes, where the simplest model is the film layer model [18]. In this model it is assumed that the concentration gradient resides entirely in a thin stagnant film (mass transfer boundary layer) with a thickness δ around the membrane. Fig. 1 shows a schematic representation of the one-dimensional fluidized bed membrane reactor model considering a film layer around the membrane.
Furthermore it is assumed that: -Bulk concentration is constant in the radial direction -There is no axial convection in the film layer -Dispersion only occurs in the radial direction (axial dispersion is neglected) -The flux is independent of the radial distance from the membrane -Isobaric and isothermal conditions.
Using the before-mentioned assumptions and Fick's law to describe the H 2 flux through the boundary layer while accounting for the induced net (convective) drift flux, the H 2 flux towards the membrane can be written as: with the mass transfer coefficient from the bulk to the membrane wall defined as which can be determined from a Sherwood correlation: where d H is the hydraulic diameter of the reactor (d H = d reac − d m ). Thus, for the thickness of the film layer: In literature, no Sherwood correlation was found that can describe the gas mass transfer from the bulk of a fluidized bed to an immersed wall inside the bed. Moreover, also no generally applicable correlation for the radial gas dispersion in fluidized beds is available. Most of the proposed correlations in literature are derived to predict the solids dispersion inside a gas-solid fluidized bed. They are 1D type equations that can describe the axial/radial movement of the solids inside a fluidized bed at the investigated operating conditions [19][20][21][22][23][24][25].
On the other hand, the equations found for gas dispersion in fluidized beds were derived for risers, circulating systems or fast fluidized beds, and for operating conditions that were often orders of magnitude higher than the superficial gas velocities used in our experiments [26][27][28][29][30][31][32]. Furthermore, these equations are not useful for CFD models, because they do not contain local and instantaneous gas and solids properties. Therefore, the radial dispersion in the densified zone of the mass boundary layer close to the membranes is estimated using the correlation by Tsotsas and Schlünder [33] for the dispersion coefficient in packed beds (see Table A.1 in Appendix A). It should be noted that the dispersion coefficient is likely to be somewhat under-predicted.
The mass transfer of hydrogen through the selective dense Pd/Ag layer of the supported membrane is described with the solution diffusion mechanism. Following Sieverts' law [34], the flux through the dense layer is proportional to the difference between the square-root of the H 2 partial pressure at the retentate side (reaction zone) and the permeate side (inside the tubes) of the membrane. The diffusion through the selective dense layer is considered as the rate limiting step for H 2 permeation and it is assumed that there is no concentration gradient (nor pressure gradient) across the porous ceramic support layer of the membrane, and also mass transfer limitations at the permeate side are assumed to be negligible (Fig. 2). These assumptions are valid in this work because the selective Pd-Ag layer was applied on the outer side of the asymmetrical porous tube. Thus, there is no concentration gradient over the porous support, since on the permeate side only virtually pure H 2 is present (the ideal perm-selectivity was in the order of 5000), and there can only be a very small pressure gradient over the porous support following the Dusty Gas model. It will be shown that the 1D model can already well describe the experimental flux when using the experimentally determined pressure at the permeate side, so that the pressure drop over the porous support can indeed be neglected, which corresponds very well with the findings by Caravella et al. (2016) [35]. If necessary, the model could be extended to account for these factors.
The membrane flux is thus described by Sieverts' law [34]: in which P m is the membrane permeability, P m0 is the permeation constant, E a is the membrane activation energy and t m is the membrane selective layer thickness. Furthermore, it is assumed that there is only mass transfer from the bubble phase to the emulsion phase, not directly to the membrane, because of the relatively small bubble hold-up in bubbling fluidized beds and especially near the vertically immersed membrane tubes. From the emulsion phase the hydrogen transfers to the film layer and from there it permeates through the membrane. Therefore the component mass balance for the bubble phase reads: The rise velocity of bubbles in a swarm (u b rise , ), the bubble fraction bed membrane reactor model. It is assumed that the gas flow rate in the bubble phase (u b ) is equal the excess gas flowrate above that is required to keep the emulsion phase at minimum fluidization velocity (u mf ). There is only mass transfer between the bubble phase and the emulsion phases (no direct mass transfer from the bubble phase to the membrane surface is considered). ( f b ) and the bubble to emulsion phase mass transfer coefficient K be are determined from correlations reported in [36]. The total superficial velocity in CSTR number n (u tot n , ), is calculated by subtracting the flow through the membrane from the axial flow in CSTR number n − 1: where A m is the surface area of the membrane and A reac is the cross sectional area of the reactor.
The emulsion phase exchanges hydrogen with the bubble phase and transports it via the film layer to the membrane wall. This can be described as: For each CSTR, one value for X m and X e will be calculated representing the average concentration of H 2 in that CSTR. The flux entering the film layer should be equal to the flux through the membrane, thus: An overview of all the hydrodynamic parameters is provided in Table A.1 (Appendix A). For a detailed discussion on the model equations and assumptions the interested reader is referred to [37].

Two-fluid model
To supplement the one-dimensional phenomenological model with an estimate of the thickness of the mass transfer boundary layer, simulations using the Two-Fluid model (TFM) have been performed, using OpenFOAM twoPhaseEulerFoam version 2.3.1. This solver has been extended with gas-phase species balance equations and realistic membrane models to simulate the selective extraction of hydrogen.
The TFM considers the gas and solids phases as interpenetrating continua. The governing and constitutive equations are presented in Table B.1 (see Appendix B). The gas phase is described as an ideal gas with Newtonian behavior, whereas the rheology of the solids phase is modeled with the Kinetic Theory of Granular Flow (KTGF). Extraction of mass via the membrane is accounted for with a source term (S m ) in the gas phase continuity equation.
The drag between the solids and the gas phase is calculated with the Gidaspow drag model [38], which combines the drag model of Ergun [39] and Wen & Yu [40]. Ergun's model is valid for high solids hold-ups (20% and higher) and Wen & Yu's model is valid at lower solids holdups (below 20%). The drag coefficient C d is determined based on the particle Reynolds number.
To approximate the rheological properties of the particulate phase in a fluidized bed, the KTGF closure equations are used. The closure equations used in this work are presented in Table B.2 [41]. A number of closure equations were not available in the original OpenFOAM TFM, so they were added to the model. Further details on the TFM and KTGF can be found in literature [38,[42][43][44][45][46]. Detailed information on the OpenFOAM TFM specifically has also been published by other authors [47,48].
To model mass transfer phenomena and extraction of hydrogen via membranes, a hydrogen species balance was added to the TFM (Eq. (14)). The effect of the membranes on the system was taken into account via the source term, S m , which is applied to the computational cells adjacent to a membrane boundary (illustrated by the red cells in Fig. 3). The source term in Eq. (14) is the membrane flux calculated with Sieverts' law, multiplied by the boundary cell's area A cell , divided by the cell volume V cell , see Eq. (15). This approach to simulate permselective membranes was also used by Coroneo et al. [49].
Previous research has shown that extraction of gas from a fluidized bed can create densified zones which may affect the flow patterns of the solids [16,50]. In the case of selective hydrogen extraction, the removal of momentum from the system is expected to have a limited effect due to low molecular weight of hydrogen. However, when modelling extraction or addition of a component with a higher molecular weight, the extraction of momentum may become more significant. Therefore, a boundary condition for the momentum balances was modified in the TFM which accounts for the extraction of momentum due to the membrane permeation. The boundary condition effectively imposes a velocity u m , whose magnitude is correlated to the extracted mass source term S m , and in the normal direction to the membrane boundary as given in Eq. (16). The velocity is always imposed normal to the membrane surface and ensures that momentum is extracted from the momentum equations in Table B The experimental setup is a cylindrical fluidized bed reactor with a single submerged membrane in the center of the reactor. This system was approximated with a 2D simulation where the bottom and top of the membrane are the axial height of the inlet and outlet of the TFM simulations. A sketch of the experimental set-up and how it has been approximated with the model is presented in Fig. 4. Hydrogen was extracted via the left boundary, to which the membrane velocity boundary condition described by Eq. (16) was applied. On the right boundary a no-slip condition was imposed. For the solids phase, a Johnson & Jackson partial slip boundary condition with a specularity coefficient of 0.50 was applied on both the left and right walls (see Table B.3 in Appendix B).
The settings for the vertical membrane simulations are presented in Tables B.4 and B.5 (in Appendix B). The domain width is equal to the radius of the experimental reactor and the domain height is equal to the membrane length. The selected grid was 0.5625 by 0.5625 mm, which is sufficiently fine to yield converged solutions. Temporal discretization was done with the second order Crank-Nicolson scheme. A combination of two second order schemes, the Gauss linear scheme and the Van Leer scheme, were used for spatial discretization. The TFM simulations were performed at three different hydrogen molar fractions at the inlet and four different reactor pressures at the outlet. The 2D Cartesian approximation of the cylindrical fluidized bed membrane reactor was applied because of its simplicity and reasonable simulation times required to obtain the results. It is an often used approximation when simulating fluidized beds [51][52][53]. However, a number of phenomena will not be simulated fully realistically with this 2D approach. To take into account all the hydrodynamic effects that occur in the fluidized bed with an immersed membrane, such as bubbles passing by around all sides of the membrane, 3D cylindrical grids are required with approximately 0.5·10 6 -1·10 6 computational cells, which is not in the scope of this study. Nonetheless, fast X-ray analysis of cylindrical fluidized beds with inserted permeating internals performed by Helmi et al. (2017) [54] showed that most of the bubbles are pushed away from the internals by the solids, which is similar to what was found in the 2D simulations.
Furthermore, in the Cartesian 2D approach, the increase in radial area when moving from the membrane surface towards the reactor wall is not taken into account. Therefore, the concentration difference between the membrane and bulk will be overestimated in the 2D Cartesian simulations compared to a 3D cylindrical system. To take into account the dependency of the hydrogen flux on the radial position, the molar balance in cylindrical coordinates can be integrated from the membrane surface (r m ) to the radial positions where the bulk concentration (r bulk = r m + δ). The result can be found in Eq.
This means that the film layer thickness, δ, can be estimated with Eq. (18), which relates the 2D Cartesian TFM film layer thickness to the actual film layer thickness with radial dependence. The TFM film layer thickness, δ TFM , can be obtained by using Eq. (4).
In this work, the 2D approach serves as learning model to qualitatively understand how concentration polarization can manifest itself in cylindrical fluidized bed membrane reactor systems. A full 3D approach is required to capture all the details of the phenomena occurring in cylindrical fluidized beds with immersed membranes. A possible alternative approach to full 3D cylindrical fluidized bed reactor simulations, while still partly accounting for the geometrical shape of cylindrical beds, is the 2.5D approach proposed by Li et al. (2015) [55]. This 2.5D approach would be computationally much more efficient than the 3D approach, but it remains to be investigated whether the 2.5D approach fully and quantitatively resolves all the required hydrodynamic features that prevail in a 3D cylindrical fluidized bed with immersed membranes. A second alternative approach to full 3D cylindrical simulations is simulating the full diameter of the fluidized bed in 2D and placing the membrane in the middle of the bed. This approach has been tested before adopting the approach used in this work, and it resulted in unrealistic gas flow profiles, where the gas regularly flows downwards around the membrane and drags bubbles down past both sides of the membrane. This gas and bubble down flow in turn causes the concentration profiles between the wall and the membrane to become more diffuse, so a stable bulk concentration far from the membrane is never reached.

Experimental
Pt/Al 2 O 3 particles with an average particle size of 200 µm and density 1400 kg m −3 were used for the experiments (provided by JM®). Detailed information on particle size measurement, minimum fluidization velocity and Geldart classification can be found in [13]. In a previous study it was ensured that particles do not chemically interact with the Pd membrane surface [56].
A 365 mm long cylindrical stainless steel tube with an inner diameter of 45 mm was used for the experiments. The gas distributor was a porous stainless steel plate of 40 µm pore size. Two thermocouples were placed inside and outside of the membrane (close by the surface of the membrane). In the fluidized bed experiments, 180 g of Pt/Al 2 O 3 particles were integrated inside the reactor to ensure full immersion of the membrane at minimum fluidization conditions. For more detailed information on the experimental setup, see [13].
In the center of the reactor a 113 mm long Pd 0.85 -Ag 0.15 based membrane supported on porous Al 2 O 3 (100 nm pore size at the surface) was placed. The outer diameter of the membrane was 1 cm and was placed 3 cm above the distributor plate. The supported membrane was fabricated with an electroless plating technique with an average selective layer thickness of 4.5 μm along the membrane. The membrane was first integrated into the reactor module without catalyst particles to activate the membrane (see [56]). Subsequently, the membrane permeation properties (P m0 andE a ) were characterized at different retentate side pressures (1.2-1.6 bar) and at different temperatures (350, 400, and 450°C) under pure hydrogen flow (P m0 = 1.76·10 −8 mol m -1 s −1 Pa −0.5 , E a = 7.1 kJ mol -1 ). The obtained permeation properties from the experiments were used in the models (phenomenological model and TFM) to describe the H 2 flux through the membrane. The permeate side was kept at 1 bar for the entire characterization period. To ensure that the membrane was leak tight, the nitrogen leakage rate was monitored during the experimental work at identical operating pressures (measured average ideal H 2 /N 2 selectivity was 5000).
After the characterization procedure, experiments were performed with binary gas mixtures of N 2 and H 2 to quantify the concentration polarization effect. Experiments with pure hydrogen at the inlet were used to monitor the stability of the membrane over time. All the experiments were performed at 400°C and hydrogen mole fractions of 0.1, 0.25 and 0.45 were used. Finally, the relative fluidization velocity (u/u mf ) was varied between 1.3 and 3.3 and the membrane performance was measured at constant pressure of 1.3 bar at the retentate side.

Results and discussion
The experiments reported hereafter were performed for various H 2 mole fractions. For every mole fraction experiments were carried out at different H 2 partial pressure differences across the membrane. First the experimental results will be compared with the results obtained with the TFM to obtain a proper estimation of the radial dispersion coefficient for the fluidized suspension. As described in Section 2, in the phenomenological model the concentration polarization is modelled by assuming a mass transfer film layer with thickness δ around the membrane with an external mass transfer coefficient of k d . The thickness of the film layer will be determined using the optimized dispersion coefficient in the TFM.
The effect of densified zones on concentration polarization is looked into and the effect of the hydrogen mole fraction and reactor pressure on the boundary film layer thickness will also be investigated. Subsequently, results from the phenomenological model without considering concentration polarization (referred to as 1D) and with accounting for concentration polarization (indicated by 1D/k d ) will be compared with experimental results for identical conditions. Finally, it will be discussed whether the bubble-to-emulsion, emulsion-to-membrane or the mass transfer across the membrane is the rate limiting step for gas extraction in fluidized beds.

Film layer thickness
In order to use the 1D/k d model which is developed in this work, the thickness of the film layer needs to be estimated. To help estimating the magnitude of δ, and the corresponding radial dispersion coefficient D r , the Two-Fluid Model (TFM) was used. To the authors' knowledge, there are currently no generally accepted relations that describe the local radial dispersion in fluidized beds with internals as a function of (amongst others) the local solids hold-up. To estimate a minimum value for the radial dispersion in the densified zones close to the membranes, we have used a correlation for packed beds. The equation of Tsotsas and Schlünder in Table A.1 was used to calculate an average radial dispersion coefficient at solids hold-ups between 0.05 and 0.60 and at an (interstitial) gas velocity of 0.1 m/s, yielding for the considered conditions a dispersion coefficient in the order of 5·10 −5 to 1·10 −4 m 2 s −1 . Variation of the gas velocity and particle diameter within the experimental range had only a very limited effect on the estimated dispersion coefficient. Thus, we have found that using a single constant dispersion coefficient based on the binary gas diffusion coefficient estimated with Fuller's equation was sufficiently accurate for the simulations performed in this work. Fig. 5 shows the computed averaged hydrogen fluxes versus the difference in the square-root of the hydrogen partial pressures across the membrane at various combinations of inlet compositions and reactor pressures for the experiments, the TFM with a dispersion coefficient of 5·10 −5 and 1·10 −4 m 2 s −1 , and the 1D model. According to Fig. 5, a very good match between TFM and experimental observations was obtained when using a radial dispersion coefficient of 1·10 −4 m 2 s −1 . Therefore, and due to the fact that no validated correlation exists for the radial dispersion coefficient inside fluidized bed membrane reactors, this value was used. The difference between the 1D model predictions and the experiments/TFM shows that the 1D model over predicts the hydrogen flux and does not take concentration polarization into account, because the 1D model does not simulate the concentration drop near the membrane surface.
Each set of experimental data (separated by the dashed boxes) in Fig. 5 was measured at different moments in time. On each day and before each experiment, the membrane performance (permeability and ideal H 2 /N 2 selectivity) was checked to ensure it was stable throughout the entire experimental work. The data obtained from the daily tests shows that the permeability of the membrane was slightly fluctuating around an average value. On the other hand, the TFM results were obtained using one average experimentally obtained permeation rate.
Due to these minor fluctuations in membrane behavior throughout the experimental program, a slight difference between model predictions and the experimental observations can be observed as both over and under predictions.
Investigating the concentration profiles in the vicinity of the membrane as computed by the TFM, the concentration of hydrogen significantly decreases from a bulk concentration to a minimum value close to the membrane for all the cases. This confirms the existence of a mass transfer boundary layer near the membrane imposing a mass transfer resistance from the bulk of the fluidized suspension to the surface of the membranes. Simulations were performed for different inlet H 2 mole fractions of 0.1, 0.2, 0.45 and 1 to investigate the thickness of this boundary layer for various operating conditions (see Figs. 6 and 7). The computed results clearly show a thinner film layer at the bottom of the membrane that increases significantly as a function of the axial position. This shows that the assumption of a film layer with a constant thickness is obviously a simplification. The description could be extended using boundary layer theory to account for this, but the results shown later will show that the assumption of a constant film layer thickness is sufficient for this system.
The TFM film layer thickness was calculated with Eq. (4). The fluxes for the three inlet mole fractions at 1.5 bar pressure are presented in Fig. 8 and were used for the calculations. The flux strongly reduces within the first 2 cm of the membrane, and then stabilizes, which shows that using a single value for the film layer thickness is a good initial estimation. When using all the flux values to calculate the film layer thickness and averaging the δ TFM values, the TFM film layer thickness is 1.55 cm. When using the average of the stable flux values above z = 4 cm, the average TFM film layer thickness was found to be 1.94 cm. The corrected film layer thicknesses δ, calculated via Eq. (18), are then 0.96 cm and 1.14 cm respectively. In the 1D/k d model, a δ of 1 cm was used. The results for all film layer thickness values are summarized in Table 1.

Densified zones
The two-dimensional TFM simulations were used to investigate the formation of densified zones near the membrane and their effect on concentration polarization. Non-consecutive snapshots of the instantaneous hydrogen mole fractions and gas bubbles (defined as regions with a gas porosity above 0.85) show that the bubbles do not  come close to the membrane (Fig. 9). The hydrogen molecules therefore have to depend on diffusion to reach the membrane, the distance of the bubbles to the membrane is too large to convectively refresh the hydrogen at the membrane. Fig. 10 presents the spatial average in axial direction of all timeaveraged solids hold-up profiles for a fluidized bed injected with binary gas mixtures with 10, 25 and 45 mol% hydrogen. The 25 mol% case was also performed for the same bed without hydrogen extraction. When hydrogen is extracted, the solids shift more towards the membrane. At higher hydrogen molar fractions, the solids hold-up near the membrane increases slightly compared to lower hydrogen molar fractions, because the momentum flux of the hydrogen towards the membrane is higher.
The solids hold-up at the wall opposite to the membrane simultaneously decreases by about 1 to 2%, indicating that the solids shift slightly more towards the membrane at higher extraction fluxes. However, when no extraction takes place, the solids hold-up near the membrane and right wall is higher than for the case with extraction. The extraction of hydrogen thus did not significantly alter the extent of the densified zones, and these small changes in the solids hold-up near    To better elucidate the effects of extraction, in Fig. 11 the gas velocity profiles in an empty membrane tube with and without hydrogen extraction are compared, where all the conditions were kept the same as for the FBMR case with 25 mol% hydrogen at 1.5 bar pressure but without solids. A clear shift of the gas velocity profile towards the membrane was observed, which explains the slight reduction in the solids hold-up near the membrane and the overall shift of the solids hold-up profile towards the membrane.

1D/kd model verification
In an independent experimental observation by Patil et al. [57] for a FBMR with a Pd-based membrane with much lower permeation properties (P m = 1.35·10 −12 mol m −1 s −1 Pa) in comparison with the membrane used in this work, no influence of concentration polarization was reported. Therefore, it was investigated whether the 1D/k d model is able to describe these experiments with negligibly small film layer thickness. Fig. 12 compares the reported experimental observation for a case at 400°C, with a superficial gas velocity between u 0 = 3 and u 0 = 5 cm s −1 , and the H 2 partial pressure drop was between ΔP = 0.5 and ΔP = 3 bar, with the predictions from the phenomenological model using 100 CSTR's in series for both emulsion phase and the bubble phases.
The model without considering concentration polarization (1D) predicts the experimental results very well. To validate the 1D/k d model, a film layer thickness of 1·10 −6 m around the membrane was assumed with a radial dispersion coefficient of 1·10 −4 m 2 s −1 , and the model reduces indeed to the results from the 1D model confirming the absence of concentration polarization for the membrane used in the work by Patil et al. at the specified operating condition.

Model vs. experiments
In this section results from the one-dimensional models (1D and 1D/ k d ) for the fluidized bed will be compared with the experimental observations at identical operating conditions. Experiments in the fluidized bed were performed with an inlet superficial gas velocity of u 0 = 0.05 m s −1 (u/u mf = 3.3). The membrane permeability was determined to be P m = 1.76·10 −8 mol m −1 s −1 Pa −0.5 , the operating temperature was 400°C, the reactor pressure was varied between 1.44 and 1.8 bar and the H 2 mole fraction was varied between 0.1 and 1.0, and the model parameters were set up accordingly. Fig. 13 summarizes the experimental observations in comparison with the obtained results from simulations with the 1D and 1D/k d models for different H 2 partial pressure differences. The figure clearly shows that the 1D model ignoring concentration polarization effects largely overestimates the membrane flux for all transmembrane pressure differences and the discrepancies further increase for smaller hydrogen concentrations, whereas the 1D/k d model that accounts for concentration polarization effects accurately predicts the membrane flux for all the considered cases.
Furthermore, the simulation results of the 1D/k d model for the case without bubble-to-emulsion mass transfer resistance were virtually identical to the results when bubble-to-emulsion mass transfer was taken into account, from which it can be concluded that the bubble-toemulsion mass transfer resistance is negligible compared to the external mass transfer resistance to the membrane wall for the considered cases.
The 1D model in this work uses the standard correlations from Kunii and Levenspiel for the bubble-to-emulsion phase mass transfer (K be ). The Davidson and Harrison correlation for the bubble-to-emulsion phase mass transfer coefficient was derived for single rising bubbles, so it will under predict the bubble-to-emulsion mass transfer for freely bubbling flows. In Fig. 14 the axial H 2 mole fraction profiles along the membrane length are plotted for the bubble phase, emulsion phase and at the surface of the membrane. This is shown for one selected   experiment, but all the other simulations showed the same trend. Since the concentration differences between the bubble and the emulsion phase and along the membrane are small compared to the concentration differences between the emulsion and the membrane wall, it can be concluded that the external bed-to-membrane mass transfer is rate limiting, thus the under prediction in the bubble-to-emulsion mass transfer rate will not have a significant effect on the simulation results. Medrano et al. (2017) [58] have developed a correlation for bubble-toemulsion mass transfer in freely bubbling beds, however only valid for pseudo-2D beds. To more accurately calculate the bubble-to-emulsion mass transfer, improved correlations for K be need to be developed specifically for fully 3D fluidized beds and fluidization in the presence of (permeating) internals.
To figure out the effect of inlet flow velocity on concentration polarization, experiments were performed with different inlet flow velocities for three different inlet compositions. Binary mixtures of H 2 and N 2 were chosen with a H 2 content of 10, 25 and 45% at the inlet. For each inlet flow composition, the inlet velocity was varied to investigate the performance of the model at higher inlet flow rates in the bubbling fluidization regime (Fig. 15).
According to the obtained results, in general the 1D/k d model can predict accurately the flux through the membrane for different inlet flow rates and inlet gas compositions. Considering the fact that the thickness of the film layer was considered with a constant value of 0.01 m, it can be concluded that this constant can be a good estimate for the average thickness of the film layer for a wide range of inlet flow velocities. Investigating the obtained modeling results with and without considering concentration polarization, the effect of concentration polarization becomes more pronounced for higher inlet gas velocities and lower hydrogen inlet concentrations, and the developed 1D/k d model can accurately capture this.  In general, for the phenomenological fluidized bed model a lot more research needs to be done. The influence of the membrane immersion on the hydrodynamic properties of the bed needs to be further investigated. The 1D/k d model gives a very good prediction of the experimental observations when using a radial dispersion coefficient estimated from TFM simulations. However, an accurate correlation for the radial dispersion in fluidized beds (and preferably accounting for the presence of immersed objects) would facilitate the modelling. Another complicating factor is that the film layer thickness is quite large compared to the reactor diameter (δ = 0.01 m and d r = 0.045 m). It should be noted that for this research a lab-scale reactor with a relatively small diameter was used, in larger reactors the effect of the film layer may be less pronounced. On the other hand, often membrane modules are inserted into the bed, where the effect of the presence and permeation through neighboring membranes might need to be accounted for in the estimation of the mass transfer boundary layer thickness.

Conclusions
A simple one-dimensional two-phase phenomenological model was developed that captures the effect of concentration polarization in fluidized bed membrane reactors (1D/k d ). In this model the fluidized bed was divided into a number of CSTR's in series for both the emulsion and bubble phase while accounting for mass transfer limitations from the bed bulk to the surface of the membranes, assuming that this occurs entirely in a thin stagnant film layer with constant thickness around the membrane. The H 2 flux through the membranes was described by Sieverts' law.
A more detailed Euler-Euler CFD model, the Two-Fluid Model, was developed in OpenFOAM where the solver was extended with species mass balance equations and membrane models to simulate the selective extraction of hydrogen. The model was used to quantify the extent of concentration polarization in a lab-scale experimental reactor and to determine the mass transfer boundary layer thickness which is required by the one-dimensional phenomenological model. Comparing the results obtained from experiments with the TFM model a very good agreement was found when an appropriate value for the gas phase dispersion coefficient was selected. The computed concentration profiles near by the membrane, confirmed the existence of a concentration boundary layer in the vicinity of the membrane that imposes a mass transfer resistance from the bulk of the fluidized bed to the surface of the membranes. Although the thickness of the film layer increases with the axial position, and decreases slightly for higher mole fractions, an average film layer thickness was estimated at 0.01 m for all the different operating conditions and was assumed constant. This film layer thickness and gas dispersion coefficient was used in the phenomenological 1D/k d model.
The results of the 1D/k d model for the membrane flux were compared with experimental observations over a wide range of inlet concentrations, operating pressures and inlet gas velocities, and a very good agreement was found, despite the fact that the film layer thickness was assumed constant. It was also found that the bubble-to-emulsion phase mass transfer limitations are much less pronounced relative to the emulsion-to-membrane wall mass transfer resistances for the investigated cases. Comparison with the 1D model results that do not account for concentration polarization, clearly indicates the very pronounced effect of concentration polarization, also for fluidized bed membrane reactors.