A generalized CFD model for evaluating catalytic separation process in structured porous materials

Anshi Hong,Zisheng Zhang,2,Xingang Li,Xin Gao,*

1 School of Chemical Engineering and Technology,National Engineering Research Center of Distillation Technology,Collaborative Innovation Center of Chemical Science and Engineering (Tianjin),Tianjin University,Tianjin 300072,China

2 Department of Chemical and Biological Engineering,University of Ottawa,Ottawa K1N 6N5,Canada

Keywords:CFD (computational fluid dynamics)Hybrid multiphase model SiC foam Mass transfer Catalytic distillation

ABSTRACT A hybrid multiphase model is developed to simulate the simultaneous momentum,heat and mass transfer and heterogeneous catalyzed reaction in structured catalytic porous materials.The approach relies on the combination of the volume of fluid (VOF) and Eulerian-Eulerian models,and several plug-in field functions.The VOF method is used to capture the gas-liquid interface motion,and the Eulerian-Eulerian framework solves the temperature and chemical species concentration equations for each phase.The self-defined field functions utilize a single-domain approach to overcome convergence difficulty when applying the hybrid multiphase for a multi-domain problem.The method is then applied to investigate selective removal of specific species in multicomponent reactive evaporation process.The results show that the coupling of catalytic reaction and interface species mass transfer at the phase interface is conditional,and the coupling of catalytic reaction and momentum transfer across fluid-porous interface significantly affects the conversion rate of reactants.Based on the numerical results,a strategy is proposed for matching solid catalyst with operating condition in catalytic distillation application.

1.Introduction

Reactive distillation (RD) has been known as one of the most important process intensification technologies for decades [1-5].By means ofin situproduct removal,it provides considerable improvements in chemical manufacturing,and substantially large reduction in equipment volume,energy consumption,and waste formation.Although reactive distillation has already been successfully implemented in dozens of industrial applications,like synthesis of methyltert-butyl ether (MTBE),alkylation of butane to isooctane and hydration of ethylene oxide to ethylene glycol,design and scale-up of new reactive distillation system still poses as a great challenge[6-10],which requires the knowledge of reaction kinetics,phase equilibrium and packing characteristics.Here,we focus on reactive distillation in structured packing,like Katapak-S[11],or the promising solid foam catalyst[12,13],where the liquid film is flowing down a porous material and exposed to a countercurrent gas stream.

As an alternative to the conventional equilibrium stage model for RD process,the non-equilibrium stage model tries to include the finite mass transfer rates across the gas-liquid interface into reactive distillation modelling,thus accurately predict conversion rate and selectivity instead of fitting process simulation results to experiment results by adjusting equilibrium relevant variables.However,it is rather complex to predict the mass transfer rate inpacking that is formed from porous materials,because both the gas-liquid flow and the mass transfer must be determined simultaneously.Moreover,the coupling effect between reaction and separation is barely clarified for liquid film flowing down a catalytic porous material.Thus,the simultaneous prediction is indeed desirable for the design and operation of efficient multifunctional distillation columns or for any kinds of application that reactive mass transfer occurs in complex geometries.

Several experimental mass transfer studies have been conducted to develop correlations for catalytic structured packing[14,15].However,in most cases the experimental studies have been carried out with academic fluids like water,ethanol,oxygen,air,etc.In this context,numerical simulations can be a complementary tool in further exploring the complex interactions between gas-liquid equilibrium,transport phenomena and chemical kinetics [16-18].

In recent years,several works can be found in the literature,which dealt with numerical simulations of multiphase flow in regular structured packing,and found good agreement between the experiment and the simulation results.The volume of fluid(VOF) method has been widely exploited in capturing the sharp phase interface in liquid film flow[19-21].Meanwhile,researchers have tried to modify the VOF method to include thermodynamic equilibria at the phase interface[22-24],and extend its application to evaporating flows[25-27].It is shown that even though absorption process can be elegantly modeled by the VOF-based method,the model formulation becomes extremely complicated when it comes to evaporation problem,in which empirical parameter or additional equations must be introduced,due to the shared velocity,temperature and species concentration formulation of the VOF method.

On the other hand,while there is computational fluid dynamics(CFD) study concerning homogeneous reaction in liquid film flow on regular structured packing [22,23],to our knowledge,there is no CFD contributions concerning reactive mass transfer over heterogeneous catalysts.Film flow over porous material has not yet drawn much attention.Since experimental measurement using techniques like PIV is hard to conduct due to the complex structure of porous material,theoretical studies play a significant role in investigating the film stability [28,29] and the momentum transport across the fluid-porous interface [30,31].Meanwhile,only a few researches can be found concerning the evaporation of liquid falling film over porous materials[32-35],which only assume that gas-liquid interface is invariantly located at the surface of the porous layer,and thus neglects the momentum balance at the fluidporous interface and the motion of the gas-liquid interface.However,referring to experimental studies on the catalytic packing Katapak-S [11] and the structured SiC foam packing [12],as well as theoretical studies of thin film on porous materials [28-31],it is more physically appropriate to assume gas-liquid interface as a free surface.

It can be concluded that existing models are not fully formulated for the mass transfer in catalytic structured packing.Comprising and adapting vital concepts from the literature,we have already utilized a hybrid multiphase method,namely the multifluid VOF,to calculate the momentum and mass transport of liquid film flowing over porous material [36].The aim of this work is to further develop the hybrid multiphase method for calculating coupled transport phenomena and catalytic reaction simultaneously.This paper is decomposed as follows.Firstly,the details of physical process and the hybrid model are presented in Section 2.Then we explore some basic characteristics of evaporation without chemical reaction of liquid film flowing over porous material in Section 3.1,before we finally present results of catalytic evaporation in Sections 3.2 and 3.3.

2.Mathematical Models

This study deals with two-dimensional countercurrent gas-liquid flow in a thermally insulated space,and the porous medium associated with the study is displayed in Fig.1(a)and(b).A porous layer of 1 mm constant thickness lies in the left side of the flow field as displayed in Fig.1(c).The widths of the liquid inlet and gas inlet are 1.5,10 mm,respectively,and the total height of the domain is 120 mm.The physical process illustrated in Fig.1(c)can be abstracted out of various industrial applications,like the distillation process happened in the SiC solid foam structured packing [12],and the catalytic distillation process happened in the Katapak-S packing [14].Major steps involved in the process are summarized as follows:

Fig.1.Photos of the open-cell of (a) structured compressed SiC foam sheet [12],(b) pre-compressed SiC foam,and (c) schematic diagram of the computational domain.

1.Species are transferred from the bulk of liquid film to the phase interface.

2.Species is evaporated into vapor at the phase interface if the species is volatile.

3.Species are transferred from phase interface to the bulk of gas phase.

4.If the porous material carries catalyst,the reactants are consumed inside the porous zone,and the products are transferred across the fluid-porous interface into the bulk of liquid phase flowing outside the porous zone.

For the record,the boundary that separates porous zone and non-porous zone is usually referred as ‘fluid-porous interface’ in related studies [30,31].

Traditionally,when the flow field has discontinuous spatial property,the model can be solved with composited computational mesh,namely to develop meshes for porous and non-porous domains separately,and then merge them together to create a multi-domain mesh for the hybrid model.An interior boundary condition between different domains is bound to exist in such mesh arrangement,and the number of equations implicitly doubles since each domain with different spatial property has specifically formulated equations.As the equations are solved individually for each domain,meanwhile the solver has to keep field properties consistent in the interior boundary,the computation is nearly impossible to converge.In order to maintain the robustness of the hybrid model for multi-domain problems,a single domain approach is adopted.Several profile and property functions concerning the spatial properties of the flow field are developed and implanted.Source codes of the functions are shown in Supplementary Material S1.

The conservation equations of mass,momentum,heat and species are in the same form for different phases.By using αqand γ to denote phasic volume fraction forarbitrary phaseqand porosity respectively,the conservation equations of mass and momentum are presented below:

where vqis the velocity of phaseq,and ˙mcharacterizes the interphase mass transfer.

wherepis the pressure shared by all phases,τqis the stress-strain tensor of phaseq.

Here,vintis the interface velocity,defined as follow.If liquid phase mass is being transferred to gas phase,vint=vl;if gas phase mass is being transferred to liquid phase,vint=vg.

The last four items in Eq.(2) are the gravity force,the surface tension force,the interaction force between phases,and the viscous and inertial resistance of porous material,respectively.

The local density and viscosity near the interface are calculated from volume average properties:

The permeabilityC1and the inertial loss coefficientC2of compressed SiC foam are derived from the Ergun equation [37]:

It is concluded from various experiments that solid foam remains isotropic and homogeneous after compression [38-41],so by assuming the strut diameterdsremains constant during compression,γ anddpare defined as follow:

The conservation of energy without species diffusion or viscous dissipation in the Eulerian-Eulerian multiphase model can be written as:

Here,Sqis the volumetric heat sources that don’t include the heat generated by finite-rate volumetric or surface reactions,since the species formation enthalpy is already included in the total enthalpy calculation.˙His the interface enthalpy,in evaporation process,it is the enthalpy of vapor at the temperature of liquid film.˙Qis the heat exchange between phases,which must comply with the local balance conditiontakes the form:

hTis the evaporative heat transfer coefficient,and the NusseltnumberNulis specified by Ref.[42]:

whereKais the Kapitza number.

Aintis the interfacial area concentration,or specific surface area,defined as:

As the interfacial transport phenomena are more related to the liquid film surface than the bulk of the liquid film,and there is literature that indicates the actual velocity at the interface is an adequate parameter for predicting the mass transfer rate in the liquid side [22],by comparing liquid film flowing over porous layer to that flowing over rigid wall,the characteristic dimensiondl,namely equivalent film thickness,is related to the liquid film surface velocity according to the Nusselt’s theory [43,44]:

The treatment of species transfer is based on the interface species mass transfer approach embedded in the Eulerian-Eulerian framework.The transport equation forthe local mass fraction of speciesiin phaseq,is:

The three items in the right side of Eq.(15) are the diffusion flux,the net production rate of homogeneous species through volumetric reaction in phaseq,and the source term for interface mass transfer,respectively.

The cascade reaction evaluated in our study is written as follows:

where the three species in liquid phase are pseudo species of the same physical properties.Detailed mechanism of the above reactions is defined in Supplementary Material S2.The production rate from a single reaction is defined as follow:

Thus,the net production rate of the main product is written as:

wherekmis the interphase volumetric mass transfer coefficient,andis the equilibrium mass concentration of speciesiin phaseq,which can be determined by equilibrium relations,like Henry’s law,Raoult’s law,etc.In our simulation,we assume that only the main product is volatile,and enforce equilibrium only at the phase interface.Equilibrium ratioKiat the interface is expressed by the Raoult’s law:

whereis the saturated vapor pressure in Pa.

Providing the transport resistances depend on the rates of diffusion through the phases on each side of the interface,the interface mass transfer rate is:

Mass transfer coefficientkmis defined as:

The mass transfer coefficients can be computed from Sherwood number correlations that are concluded from experiments with similar hydrodynamic condition to our simulation.The Sherwood correlations for gas [45] and liquid phase [46] are written as follows respectively:

3.Results and Discussion

The described equations were solved by the multi-fluid VOF approach embedded in software package Fluent 17.0 with double precision.The phase coupled SIMPLE algorithm was used for the pressure-velocity coupling,and the Geo-Reconstruct scheme was used for interpolation near the interface.Since we modeled laminar flow with a quadrilateral mesh,and the numerical diffusion would be naturally low,the convection and diffusion terms in all equations were discretized by first order upwind scheme.The adaptive time stepping scheme was used by setting the global courant number to 1.All inlets and outlets in the domain were set as velocity-inlet and pressure-outlet,respectively.In addition,all cases converged more easily if the temperature was limited to a reasonable small range according to the temperatures of liquid and gas phase.Unless otherwise stated,the temperatures for the gas and liquid at the inlets were 303 and 293 K,respectively,the velocity magnitude normal to the gas inlet was 2 m·s-1,the saturated vapor pressure for main product was 100 kPa,the reaction rate constant was 0.08 s-1,and the initial composition for liquid phase was 60% reactant,40% main product.Firstly,the flow and volume fraction equations were solved to obtain an initial solution,and then the species concentrations were adjusted to the same values as the inlets,finally the species equations were turned on for complete simulation.The simulation was considered to reach pseudo-steady state when the species mass fractions in all outlets were approximately constant.

The flow domain is discretized with non-uniform hexahedral mesh,and a grid independence test was conducted to determine a reasonable mesh while maintaining grid independent predictions.Temporal evaluations of the evaporation rate (defined later)are presented in Fig.2(a) for three grid resolutions.Simulations with different grid resolution predicted similar value.Hence,a resolution consisting of 51686 elements was selected as an optimum mesh.The node spaces in both side of the fluid-porous interface grow linearly from 0.02 mm to less than 0.05 mm,as shown in Fig.2(b).

Fig.2.(a) Evaporation rate for three grid resolutions,(b) sketch of local mesh configuration.

As the back side of porous structured packing has no constrain on the liquid phase,the slip boundary condition was defined on the wall,unless the porosity was defined as 1,namely the liquid film flowed over a rigid wall.As shown in Fig.3(a),the existence of the porous layer significantly increases the thickness of liquid film.The velocity profile of the bulk liquid phase flowing inside the porous layer is close to that of plug flow,and due to momentum exchange,a rapid yet continuous change of velocity occurs near the fluid-porous interface.Meanwhile,the velocity profile of the liquid film flowing outside the porous layer is close to that of liquid flowing over a rigid wall,which supports us for using Nusselt’s theory to relate the equivalent film thickness to liquid film surface velocity.In Fig.3(b),it’s found that neither the liquid nor the gas flow regime is closed to ideal plug flow,which indicates that the simultaneous prediction of transport phenomena and reactions is indeed needed to further develop and design efficient catalytic distillation equipment.

Fig.3.Velocity profiles at x coordinate 0.06 m calculated for different configuration (Rel，inlet ≈837): (a) comparison of velocity profiles in liquid phase fordifferent fluid domains;(b) velocity profiles in liquid and gas phase of fluid domain with a porous layer.

3.1.Mass transfer characteristic of evaporation without reaction

Properties of the porous medium used in Section 3 are listed in Table 1.The evaporation rate in Fig.4 is expressed as follows:

Table 1Properties of the porous medium

The influence of entrance effects on the determination of the evaporation and reaction rate is taken into account by taking samples from the height 10 mm and 110 mm instead of the inlets and outlets,as implemented by Schoutenet al.[47].

Fig.4.Comparison of the evaporation rate at different gas F-factor for three porosities (Rel，inlet ≈837).

As shown in Fig.3,a great amount of liquid flows inside the porous material and is far from the gas-liquid interface,where evaporation actually takes place.In result,a decrease in evaporation rate is observed when the liquid film flows over porous material and the porosity is higher.Meanwhile,although the mass transfer coefficient for gas phase increases with increasingF-factor(vgρ0.5g),namely increasing gas flow rates as indicated by Eq.(22),the evaporation rate is restricted by the mass transfer rate in liquid phase at large gas flow rates.At highF-factor,the larger mass transfer coefficient for the liquid phase overcomes the laziness of liquid flowing inside the porous material,so the evaporation rate was higher for high porosity material than the low porosity one at large gas flow rates.

Fig.5.Species profiles of liquid film at x coordinate 0.06 m for porous material with different porosities: (a) 0.76 at Rel，inlet ≈478,(b) 0.76 at Rel，inlet ≈598,(c) 0.82 at Rel，inlet ≈598,(d) 0.82 at Rel，inlet ≈717.

3.2.Effect of liquid flow rate

3.2.1.Coupling of catalytic reaction and mass transfer

Catalytic reaction and interface species mass transfer happen simultaneously yet at different locations of liquid film.While catalytic reaction happens inside the porous material,interface species mass transfer exists only at the gas-liquid interface.The simplification algorithm for equipment design largely depends on the coupling effect of catalytic reaction and mass transfer.As illustrated graphically in Fig.5(b)and(d),the species diffusion from the bulk liquid phase to the phase interface involves only a narrow layer of liquid near the phase interface,thus there exists no interaction between catalytic reaction and mass transfer.In addition,the species profiles of evaporation without chemical reaction and reactive evaporation are identical near the phase interface.Meanwhile,when the phase interface nearly overlaps the fluid-porous interface,the increase of species concentration due to catalytic reaction enhances interface species mass transfer,as shown in Fig.5(a) and (c).

3.2.2.Coupling of catalytic reaction and momentum transfer

As for catalytic reaction,we assume that the volumetric reaction rate doesn’t vary with the porosity,namely the quantity of solid catalyst.The assumption is valid as the reaction rate maybe restricted by the liquid-solid mass transfer rate,or the catalyst activity is variable by controlling the density of active components.The average mass fraction change of species per unit length and the average selectivity of the main product in Figs.6 and 7 are defined as follows:

Fig.6.Average mass fraction change of reactant per unit length at different porosity and liquid flow rate.

Fig.7.Average selectivity of the main product at different porosity and liquid flow rate.

Fig.8.Velocity magnitude at the fluid-porous interface at different liquid flow rate.

Fig.9.Average mass fraction change of species per unit length at different saturated vapor pressure and reaction rate constant(Rel，inlet ≈717，γ=0.82).For xaxis,valueof -1/lgkr decreases withdecreasingreaction rate constant,value of 1/ldecreaseswith increasingsaturated vapor pressure.The bottomx axis links to the black line,and the top x axis links to the red one.For y-axis,value of-1/lg（Δmf）increases with increasing consumption rate of the main product,which is calculated according to Δmfmainproduct in evaporation process,or Δmfbyproduct in reactive evaporation process.

High porosity is advantageous for both production and selectivity of the main product,as illustrated in Figs.6 and 7 respectively.At the same volumetric reaction rate,according to Eq.(15),the consumption rate of reactant is proportional to the porosity,namely the volume of flowing liquid,which makes up for the less reaction time at higher porosity.Meanwhile,the conversion rate per unit length decreases with increased liquid flow rate,which indicates the interaction between catalytic reaction and momentum transfer.As shown in Fig.3,due to the momentum transfer across the fluid-porous interface,the velocity profile inside the porous material is affected by the liquid flowing outside,thus the average velocity of the liquid flowing inside the porous material will increase with increased liquid flow rate,which can be told from the variation of velocity magnitude at the fluid-porous interface shown in Fig.8.In result,the reaction time reduces with increased liquid flow rate,which makes misdistribution in catalytic distillation column a new problem for conversion rate of reactants besides catalyst deactivation,and requires further study on the wetting behavior of liquid film over porous material.

Fig.10.Average mass fraction change of by product per unit length at different saturated pressure (Rel，inlet ≈717，γ=0.82).

3.3.Effect of saturated vapor pressure and reaction rate constant

When the evaporation and cascade reaction happen simultaneously,the main product is evaporated into vapor at the gas-liquid interface,and transformed into the by product inside the porous material.The consumption rates of the above two processes differ a lot,and also significantly affect the detailed design of catalytic distillation equipment.The relative relation of the consumption rate between the evaporation and reaction is illustrated in Fig.9.

As shown in Fig.9,the reaction rate is usually much larger than the evaporation rate,as the specific surface area differs more than an order of magnitude between the solid catalyst and gas-liquid interface.Meanwhile,referring to Fig.9,we can design a catalytic distillation system that requires efficient recovery of speciesin situ,by customizing the catalyst according to the acceptable range of operating pressure,temperature and equilibrium ratio.

On the other hand,it is indicated from Fig.10 that the coupling effect of reaction and interface species mass transfer is barely related to the phase equilibrium in largeRel，inlet,as the production rate of the by product hardly varies with saturated vapor pressure.

4.Conclusions

The behavior of the liquid film is a key aspect to the design of novel porous structured packing utilized in distillation and catalytic distillation processes.The main focus of this work is to develop a multiphase model for computing the detailed information about transport phenomena in and between moving fluid phases,which allows capturing the free motion of the phase interface,and reproduces both the discontinuity of species concentration at the phase interface,and the continuity of velocity at the fluid-porous interface.

The model has been applied to study non-reactive and reactive mass transfer in a falling liquid film over a porous material.The results show that liquid flowing inside the porous material has a negative effect on the interface species mass transfer,and the coupling of heterogeneous reaction and interface species mass transfer only exists when the fluid-porous interface overlaps phase interface.In addition,misdistribution in catalytic distillation column affects the conversion rate of reactants.Finally,the numerical results are helpful for screening operating condition and solid catalyst.

It should be mentioned that when the liquid flow rate is out of range,the parameters at the outlets are turbulent,and the flow field loses its stability to oscillations.Moreover,the wetting behavior remains a rarely mentioned topic for liquid film over porous material.Stability and spreading of the liquid film appear to be the tasks for the future.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors acknowledge financial support from the National Key Resaerch and Development Program of China(2019YFE0123200) and National Natural Science Foundation of China (21776202).The authors also thank the reviewers for their specific comments and suggestions.

Supplementary Material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.cjche.2021.12.004.

Nomenclature

Ainterfacial area concentration,m2·m-3

apore size before compression,m

C1permeability,m2

C2inertial resistance factor,m-1

cbulk phase molar concentration,kmol·m-3

molecular diffusivity of speciesi,m2·s-1

dpequivalent particle diameter,m

dsstrut diameter,m

dlequivalent film thickness,m

Fqmomentum resistance of porous medium,Pa·m-1

Fσ surface tension force,Pa·m-1

F-factor,(m·s-1)·(kg·m-3)0.5

g gravitational acceleration,m·s-2

Hqenthalpy,J·kg-1

hToverall heat transfer coefficient,J·m-2·s-1·K-1

Kequilibrium ratio

KaKapitza number,

Klginterface exchange coefficient,Pa·m-2

kggas phase mass transfer coefficient,m·s-1

klliquid phase mass transfer coefficient,m·s-1

kminterphase volumetric mass transfer coefficient,m·s-1

krreaction rate constant,s-1

local mass fraction of speciesi

NuNusselt number,

n surface normal,?αq

ReReynolds number,

rccompression ratio

ScSchmidt number,

ShSherwood number,

Ttemperature,K

vintinterface velocity,m·s-1

vggas phase velocity,m·s-1

vlliquid phase velocity,m·s-1

vqvelocity vector,m·s-1

xxcoordinate

yycoordinate

αqphasic volume fraction

κTthermal conductivity,W·m-1·K-1

μqviscosity,Pa·s

μlglocal viscosity,μlg=αlμl+αgμg,Pa·s

ρqdensity,kg·m-3

ρlglocal density,ρlg=αlρl+αgρg,kg·m-3

σ surface tension coefficient,N·m-1

γ0original porosity before compression

γ porosity

τqstress-strain tensor,Pa

Subscripts

e equilibrium g gas phase

int interface

lliquid phase

qarbitrary phase

sat saturation pressure

Superscripts

i,jarbitrary species