A Reynolds mass flux model for gas separation process simulation:I.Modeling and validation☆

Wenbin Li,Kuotsung Yu ,Xigang Yuan ,*,Botan Liu

1 Collaborative Innovation Center of Chemical Science and Engineering(Tianjin),Tianjin University,Tianjin 300072,China

2 State Key Laboratory for Chemical Engineering and School of Chemical Engineering and Technology,Tianjin University,Tianjin 300072,China

Keywords:Mathematical modeling Computational fluid dynamics(CFD)Computational mass transfer(CMT)Anisotropic turbulent mass diffusion Packed bed Absorption

ABSTRACT Separation process undertaken in packed columns often displays anisotropic turbulentmass diffusion.The anisotropic turbulent mass diffusion can be characterized rigorously by using the Reynolds mass flux(RMF)model.With the RMF model,the concentration and temperature as well as the velocity distributions can be simulated numerically.The modeled Reynolds mass flux equation is adopted to close the turbulent mass transfer equation,while the modeled Reynolds heat flux and Reynolds stress equations are used to close the turbulent heat and momentum transfer equations,so that the Boussinesq postulate and the isotropic assumption are abandoned.To validate the presented RMF model,simulation is carried out for CO2 absorption into aqueous NaOH solutions in a packed column(0.1 m id,packed with 12.7 mm Berl saddles up to a height of 6.55 m).The simulated results are compared with the experimental data and satisfactory agreement is found both in concentration and temperature distributions.The sequel Part II extends the model application to the simulation of an unsteady state adsorption process in a packed column.

1.Introduction

Turbulent mass diffusion in a packed column is found to be anisotropic since many experimental data[1–5]have shown that the turbulent axial mass diffusivity is far greater than the radial one.Therefore,it is suggested[3,6]that anisotropy of the turbulent mass diffusion should be taken into account for precise simulation.

Due to the complexity of fluid flow in a packed column coupled with reaction and heat effect,many researchers focused on linking engineering models with rigorous simulation tools.Among them,computational mass transfer(CMT)simulation is appropriate to be used to give more visualized predictions of hydrodynamic,heat and mass transfer for design and optimization of process equipment[7–10].The current CMT methods to simulate the separation process in a packed column can be classified into two categories.One is called the digital packing algorithm,which features direct solution of the governing equations of fluid flow without any simplifications,in a box or cylinder containing the packing elements[11,12].This method can supply the details of flow behavior and the results may be helpful for investigating the mass and heat transfer mechanism,but it needs more computational meshes and complicated boundary conditions.Among the modeling works based on this method so far,simulation was done under the condition that at most 1250 balls were structurally housed in a tube[13],and obviously,this number is quite small compared to any industrial column.The other method is simplifying the packed column as a porous medium model with consideration of radial distribution of porosity[14–19].Here,the packing is simplified,and the geometrical shape of the packing neglected.In this way the computation resource needed is greatly reduced and the performance of the whole column can be conveniently predicted.In such a model,the special condition of mass transport between phases can be described by adding a specific source term to the turbulent mass transfer equation.However,since turbulent mass diffusion is the dominating factor behind mass transport in the separation process, finding a proper method for determining the turbulent mass diffusivity is of great challenge.

In conventional CMT models,the turbulent mass diffusivity is obtained either by employing empirical correlation[17,20]or by guessing a constant Peclet number[21,22].Nevertheless,the right choice of Peclet numbers is difficult.The value of turbulent mass diffusivity obtained by using empirical correlation from an inert tracer experiment is different from that with mass transfer proceeding.To avoid the use of empirical correlation or Peclet number,Sun et al.[23,24]developed the-εcmodel for simulating the chemical engineering processes in packed and tray columns.By this model,the axial turbulent mass diffusivity is determined by solving auxiliary model equations.Although it is a sort of progress,the radial turbulent mass diffusivity is assumed to be equal to the axial one.In other words,the turbulent mass diffusivity Dt,iis assumed to be isotropic,but as mentioned before many experimental investigations[1–5]indicate that the turbulent mass diffusion is anisotropic.Therefore,the assumption of isotropic turbulent mass diffusion leads to simulation error.

In fact,there have been many attempts[2,3,6,25–27]at correlating the anisotropic turbulent mass diffusivity through an inert tracer experiment.However,a wide variation exists in the value of the turbulent mass diffusivities from different experimental correlations.Macías-Salinas and Fair[28]attribute such discrepancies to the sensitivity of the effective diffusivity to experimental errors and the use of different model concepts for estimating the turbulent mass diffusivity.Although the anisotropy of turbulent mass diffusivity can be captured to some degree by using such experimental correlations,it can be expected that use of different correlations would yield divergent predictions for concentration distributions.To reduce uncertainty of the simulation on concentration field,Yu et al.[29,30]proposed a new CMT model,so-called the hybrid Reynolds mass flux model,for the simulation of distillation and obtained the concentration distribution in the packed column.The feature of the hybrid RMF model is that the modeled Reynolds mass flux equation is used to close the turbulent mass transfer equation so that the anisotropic turbulent mass diffusion can be characterized.To the best of our knowledge,it is the first attempt at predicting and investigating the anisotropic turbulent mass diffusion using the CMT model.However,in the work by Yu et al.[29,30],the conventional k–ε model was used to evaluate the Reynolds stress for closing the turbulent momentum transfer equation.That is to say,the Boussinesq postulation along with the isotropic turbulent viscosity was still retained,which contradicts the highly anisotropic turbulence existing in most of chemical processes.Therefore,the hybrid RMF model needs to be further improved so that the Reynolds stress can be obtained more rigorously.

The aim of the present study is to propose a standard RMF model which combines the modeled Reynolds mass flux equation with the modeled Reynolds stress equation.Then,the proposed RMF model is applied to the simulation of gas separation processes,including the chemical absorption and adsorption processes.The work accomplished is presented in two parts.In Part I,the hybrid RMF model[29,30]is modified and the modeled Reynolds stress equation is incorporated to establish the standard RMFmodel.To validate the proposed model,simulation is carried out for CO2absorption into aqueous NaOH solutions in a packed column[21].Part II extends the application of the proposed standard RMF model to simulating the unsteady adsorption process in a packed column for methylene chloride adsorption on activated carbon,and compares the predictions with published experimental data[31].

2.Model Equations

The established standard Reynolds mass flux model consists of a turbulent mass transfer equation along with the modeled Reynolds mass flux equation,and the accompanied computational fluid dynamics(CFD)and computational heat transfer(CHT)equation sets.

2.1.Turbulent mass transfer equation with closures

Turbulent mass transfer equation for the concerned species:

The foregoing equation should be modeled in order that it can be suitable for computation.In this paper,the exacttransport equation is deduced and then modeled as below[29,30]:

where the constants were found to be[32]:CC1=0.09,CC2=3.2,and CC3=0.55;δlmjkis the Kronecker delta,and δlmjk=1 only if l=m=j=k;for solving the Reynolds stress,the following simplified equation based on the Boussinesq postulation was adopted in the hybrid RMF model[30]:

However,in this paper the following modeled Reynolds stress equation proposed by Chen and Jaw[32]is used so that the isotropic turbulent viscosity is no longer needed:

where ρ is the density of the fluid phase;and the model constants are[32]:Cu1=0.09,Cu2=2.3,and Cu3=0.4.

The auxiliary equations of k and ε are as follows:

k equation:

ε equation:

where the constants are[32]:Cε1=0.07,Cε2=1.45,and Cε3=1.92.

2.2.CFD equation set

Continuity equation:

Turbulent momentum transfer equation:

2.3.Accompanying CHT equation set

Turbulent heat transfer equation for the fluid phase:

where the constants are[32]:CT1=0.07,CT2=3.2,and CT3=0.5.

Energy balance equation for the static phase:

where ρsis the particle density;ksis the thermal conductivity of the solid phase;cpsis the specific heat of the solid phase;χ is the local porosity of the random packed bed,and(1-χ)is the solid phase volume fraction;andis the thermal source term.

3.Model Implementations

A number of issues are related to the implementation of the Reynolds mass flux model in chemical absorption simulation.They include evaluation of various terms,boundary conditions and numerical procedures.

3.1.Evaluation of various terms

The object of simulation is CO2absorption into aqueous NaOH solutions and the experimental data are taken from the published report by Tontiwachwuthikul et al.[21].Absorption of CO2into aqueous NaOH solution can be regarded as a gas absorption control process accompanied by a relatively slow second-order reaction[33].Detailed information about the reaction mechanism could be found in Liu et al.[34].

For the simulation of turbulent gas–liquid two-phase flow,the liquid phase model(with interaction with the gas phase)is employed,which has been shown to be more convenient and computer-time saving with sufficient accuracy provided if the gas–liquid interacting source terms can be properly expressed[35].In addition,the following assumptions are made:

(1)The absorption operation is at steady state,and the liquid is incompressible;

(2)The heat of solution and that of reaction are all absorbed by the liquid phase.The heat transferred from liquid to the gas phase is neglected.The absorption process is adiabatic,which means the heat transferred from the column to the environment is also neglected.

(3)Only the CO2component in the gas phase is absorbed by the aqueous NaOH solution;the water content in the solutions does not transfer to the gas phase.

3.1.1.Net mass increase of liquid and OH-

where aedenotes the effective area for the gas–liquid mass transfer;MCO2is the molar mass of the species CO2;is the molar concentration of CO2at the interface;and XCO2is the molar concentration of CO2in the bulk liquid,which is assumed to be zero.

The molar concentration of CO2at the interfacecan be expressed in terms of Henry's law:

where ptis the total pressure of the gas phase;yCO2is the volume fraction of CO2in the gas phase;and H is Henry's constant,which can be determined from the correlation by Danckwerts[36].The liquid phase coefficient of mass transfer kLis calculated by

The molecular diffusivity of CO2in the gas phase DGis calculated from the empirical correlation by Poling et al.[37].The molecular diffusivity of CO2in the MEA aqueous solution DCO2,Lis also determined by using the N2O analogy[38].And the enhancement factor E can be determined from the correlation by Wellek et al.[39].The effective mass transfer area aein Eq.(13)is considered equal to the wetted surface area aw,which can be calculated by the following correlation[33]:

3.1.2.Force acted on liquid phase

where FLG,ican be expressed by

where the pressure drop Δpt,caused by the interaction between the gas and liquid phases,is the sum of dry bed pressure drop Δpdand the additional pressure drop due to the liquid adhered to the packing surface ΔpL;Δpdand ΔpLare calculated by the correlations reported by Robbins[40];and Uslipis the slip velocity between the gas phase and liquid phase,which is defined as

where UGis the superficial velocity of the gas phase.

The body force FLS,icaused by the random packing is determined by the Ergun equation[41].

The liquid volume fraction β can be calculated by

where the static holdup hsof ceramic Berl saddles is 0.0317[42];for the metal Pall rings,hscan be calculated from the correlation by Engel et al.[43];the operating holdup hopis calculated by the method given by Sater and Levenspiel[44];the local porosity χ of a random packed bed is calculated by the equation of Liu et al.[45],giving a local porosity of the random packing varying throughout the packed bed.

3.1.3.Thermal effects on liquid

where HAis the heat of physical absorption and HRis the heat of chemical reaction.

The thermal source term of the static phaseis neglected in the case of absorption.

3.2.Boundary conditions

Fig.1 shows the computational domain and boundary arrangement.The boundary conditions for model equations are specified as follows when the column works in the counter-current flow mode.

Fig.1.Computational domain and boundary arrangement of the packed column.

(1)Inlet conditions

The inlet boundary conditions for the liquid phase are[34]

where the hydraulic diameter of random packing dHcan be calculated by[46]

The surface area per unit volume of packed column apcan be calculated by[41]

The inlet condition for the Reynolds stress is set to be[47]

The inlet conditions for the gas phase are

In the present study,the generalized Fick-like law[48]is adopted forandat the inlet:

where Dtand αtare recommended to be Dt= νt/0.7 and αt= νt/0.9[49],thus,

and the inlet turbulent viscosity νt,in=with constant Cμ=0.09.

(2)Outlet conditions The outlet of flow is considered to be close to fully developed,so that zero normal gradients are chosen for all variables except pressure.

(3)Wall conditions The no-slip condition of flow is applied to the wall,and the flow behavior in the near wall region is approximated by using the“standard wall function”.The zero flux condition is applied to the turbulent mass and heat transfer equations.

3.3.Numerical procedures

The model equations were solved numerically by using the commercial software Fluent 6.3.26 with the finite volume method.The SIMPLEC algorithm[49]is employed for resolving the pressure–velocity coupling in the momentum equations.Considering the flow is symmetric,only half of the vertical cross section passing the axis is simulated.

The grid arrangement of the pilot scale column of 6.55 m in height and 0.05 m in radius comprises totally 160000 quadrilateral cells.There are 2000 nodes uniformly distributed along the column height and 80 nodes non-uniformly distributed along the radial direction with higher grid resolution in the near wall region.

4.Results and Discussion

12 sets of experimental data for absorption of CO2into aqueous NaOH solutions were reported by Tontiwachwuthikul et al.[21].In the present study,the simulated results by the RMF model will be compared not only with the experimental data but also with the results obtained by using the c2-εcmodel[29,34].Since Liu et al.[34]have given the simulated results on experiments T1,T2,T7,T8,T11 and T12,in this paper only the simulations on experiments T1,T2 and T12,which cover the typical experimental conditions,are taken for comparison.

4.1.Model validations

To validate the presented Reynolds mass flux model,comparisons are made between the simulated results and the experimental data on concentration and temperature as well as the velocity profiles.

4.1.1.OH-concentration and liquid temperature profiles

In Figs.2–4,circular symbols denote the experimental measurements,while the solid lines are the simulated results obtained by using the present model.Satisfactory agreement is seen between the predicted distributions of radially averaged OH-concentration and liquid temperature by using the present model and the experimental data,which confirms the validity of the presented Reynolds mass flux model.

Figs.2(a)and 3(a)also give the predicted OH-concentration by the c2-εcmodel[29,34](dash lines),for which the Boussinesq postulation was adopted and the turbulent mass diffusivity was assumed to be isotropic.It is found that the present model provides better predictions(solid lines)on OH-concentration,which supports the statement that the anisotropic turbulent mass diffusion should be taken into account for more precise simulation.

Fig.5 shows the profiles of the simulated OH-concentration and the liquid temperature in T12.From the enlarged column section in Fig.5(a),the OH-concentration is found to be lower atthe column center than at the column wall at a fixed axial position.It is due to the reason that the velocity of the liquid phase is decreased sharply to the column wall surface,resulting in worse contact with the gas phase,and consequently less CO2would be absorbed,producing less exothermic heat of solution and reaction.It is noted from Fig.5(b)that the liquid temperature is higher at the column center than at the column wall at different axial positions.

Fig.2.Profiles of OH-concentration and temperature in T1.

Fig.3.Profiles of OH-concentration and temperature in T2.

Fig.4.Profiles of OH-concentration and temperature in T12.

4.1.2.Liquid velocity profiles

The shapes of the velocity profile at different axial positions are similar as shown in Fig.6(a),thus an example of the relative axial velocity profile at x=3.25 m is given in Fig.6(b).The velocity profile shows a significant variation which is due to the non-uniform packing porosity especially at the near wall region(see Fig.6(c)),and then becomes relatively uniform at the region of 2dpapart from the wall.This phenomenon has been confirmed by many investigators[45,50].

4.2.Discussion on anisotropic turbulent mass diffusion

Fig.7 gives the profiles of the Reynolds mass fluxpredicted by the present model in terms ofand.It is seen that although the Reynolds mass fluxes in two directions are practically in the same order of magnitude,difference in their values exist.

As mentioned in the Introduction section,the present model characterizes the anisotropic turbulent mass diffusion.Based on the Fick-like law,the turbulent mass diffusivity could be obtained by

where Dt,xand Dt,rare the axial and radial turbulent mass diffusivities respectively,and the concentration gradients in two directions are shown in Fig.8(a)and(b).

Fig.9 shows the distributions of Dt,xand Dt,ralong the radial direction.No analogy of the distribution pattern is found between the diffusivities in the axial direction and that in the radial direction,which implies that the isotropic assumption of the turbulent mass diffusivity in the-εcmodel[29,34]is unreasonable.

Fig.5.Profiles of OH-concentration and temperature in T12.

Fig.6.Profiles of liquid velocity and packing porosity at axial position x=3.25 m in T12.

Table 1 gives the comparisons of the turbulent mass diffusivity obtained by the present model,the-εcmodel[29,34]and the literature correlations.The evaluated turbulent mass diffusivities from various experimental correlations show disagreement.As mentioned above,such discrepancies are attributed to the sensitivity of turbulent mass diffusivity to experimental errors and the use of different model concepts for estimating the turbulent mass diffusivity.It is seen from Table 1 that the axial turbulent mass diffusivity is larger than the radial one no matter what correlation is used.Yin et al.[17]and Sherwood et al.[5]attribute this to the reason that the axial component of flow is larger than the radial one in the packed column.This anisotropy of the turbulent mass diffusivity is characterized by the present model.Nevertheless,isotropic turbulent mass diffusivity is used in the-εcmodel[29,34].Difference exists between the turbulent mass diffusivity obtained by the present model and that by the-εcmodel[29,34]and empirical correlations because the-εcmodel[29,34]and empirical correlations are only approximations.

Fig.7.Profiles of Reynolds mass flux at different axial positions in T12.

Fig.8.Profiles of concentration gradient at different axial positions in T12.

Fig.9.Profiles of anisotropic turbulent mass diffusivity at different axial positions in T12.

4.3.Grid independence test

To ensure that the simulated results presented in this paper are independent of the number of grid points considered,we mesh the column with different grid points.The simulated radial OH-concentration profiles with 40×1000,80×2000 and 160×4000 grid points are given in Fig.10.It is found that difference of the simulated results is significant when the grid points are lower than 80×2000;while difference is small when the number of grid points is beyond 80×2000,which means that the present study with 80×2000 grid points is sufficient for the simulation of concentration distribution.

Table 1 Comparisons between the turbulent mass diffusivity evaluated by using experimental correlations,the -εc model[34]and the Reynolds mass flux model

Table 1 Comparisons between the turbulent mass diffusivity evaluated by using experimental correlations,the -εc model[34]and the Reynolds mass flux model

Note:Peclet number based on molecular dispersion coefficient Pe m=|u|d p/D,for the present study Pe m is about 1×105;axial Peclet number for mass dispersion Pex=|u|d p/D t,x;radial Peclet number for mass dispersion Per=|u|d p/D t,r;|u|is the averaged interstitial velocity of liquid;Reynolds number Re p=ρd p U/μ,for the present study Re p is about40;U is the superficial velocity of liquid;and(av.)means the volume averaged value.

Reference Experimental correlation Value Value Pex Per Remark Pex Per D t,x× 10-4/m2·s-1 D t,x× 10-4/m2·s-1 Sater and Levenspiel[51]Pex=7.58×10-3(Re p)0.703– – 0.11– 4.38 –Saffman[52]■ – Valid for Pe m?1 0.45– 6.84 –Fetter[26] –P1 ex=1 6 ln（ 3Pe m）-■1 4■Pe m+0.055 – –17.9– 0.17 Wen and Fan[25] – Per=P1 er=■■1 2 Re p+11.4 Valid for high values of Pe m – 11.7– 0.26 Delgado[6]17.5■Pe m+0.025 Valid for 300＜Pe m＜105 0.56 38.1 5.48 0.08 Liu et al.[34]– Isotropic diffusivity–9.31(av.)Present model – – Anisotropic diffusivity – –2.31(av.) 0.60(av.)）Valid for 10＜Re p and Pe m＜106 0.65 13.9 4.69 0.22 Sahimi[53]P1 ex=25Sc1.14/■■2■Pe m+0.5 1 P1 er=45.9-33.9 exp-21Sc/■■2 1（ ■Pe m P1 ex=■Pe m+2.2■■1 P1 er=1 2■■2

Fig.10.The predicted radial profiles of OH-concentration with 80 and 150 radial grid points at x=3.25 m.

Generally,the grid to be used in the simulation of anisotropic turbulence needs to be much finer.It is necessary to test whether the use of the present grid scale is sufficient for capturing the concerned anisotropic turbulent mass diffusion.Taking Run T12 as an example,the simulated anisotropic turbulent mass diffusivities with different grid scales are given in Table 2.

Table 2 Simulated anisotropic turbulent mass diffusivity with different grid scales

5.Conclusions

A Reynolds mass flux model is presented to simulate the separation process in a packed column so that anisotropic turbulent mass diffusion can be characterized rigorously.The presented Reynolds mass flux model consists of the turbulent mass transfer equation,closing with the recently modeled Reynolds mass flux equation,and the accompanied equation sets of computational fluid dynamics(CFD)and computational heat transfer(CHT).To validate the proposed model,simulation is carried out for a randomly packed column for CO2absorption into aqueous NaOH solutions.The simulated results are compared with the experimental data and satisfactory agreement is found between them in both concentration and temperature distributions,which confirms the validity of the present model.Moreover,the anisotropy of the turbulent mass diffusivity is discussed.Comparison of simulated results with experimental data shows that the present model performs better than the model[34]with isotropic turbulent mass diffusivity.It also confirms the statement that the anisotropy of the turbulent mass diffusion should be taken into account for precise simulation and analysis of the absorption process in a packed column.

Nomenclature

ae,aweffective and wetted surface area per unit volume of packed column,m-1

apsurface area per unit volume of packed column,m-1

CC1,CC2,CC3model constants for the concentration field

Cininlet mole concentration,kmol·m-3

cpf,cpsspecific heat of fluid phase and solid phase,respectively,J·kg-1·K-1

CT1,CT2,CT3model constants for the temperature field

D molecular diffusivity,m2·s-1

dH,dphydraulic and nominal diameters of packing,m

E enhancement factor

FLG,Iinterface drag force between gas and liquid phases,N·m-3

FLS,Ibody force created by the random packing,N·m-3

g gravity acceleration,m·s-2

HAphysical absorption heat of mol CO2absorbed,J·kmol-1

HRchemical reaction heat of mol CO2absorbed,J·kmol-1

hopoperating liquid holdup for absorption process

hsstatic liquid holdup for absorption process

httotal liquid holdup for absorption process

k turbulent kinetic energy,m2·s-2

kLliquid phase coefficient of mass transfer,m·s-1

L flow rate of liquid,kg·m-2·s-1

M molecular weight of species,kg·kmol-1

PemPeclet number based on molecular dispersion(Pem=|u|dp/D)

PeiPeclet number based on turbulent mass diffusivity(Pei=|u|dp/Dt,i)

RepReynolds number based on particle diameter(Rep=ρd|u|/dp)

r radial distance from the axis of the column,m

T0ambient temperature and initial temperature of the solid phase,K

t time,min

xi,xjcoordinates in axial and radial directions

yCO2mole fraction of CO2in the gas phase

α,αt,αeffmolecular,turbulent and effective thermal diffusivities,respectively,m2·s-1

β volume fraction of the concerned phase

ε turbulent dissipation,m2·s-3

δijKronecker delta

ρ,ρsfluid phase and solid phase densities,kg·m-3

μ gas molecular viscosity,kg·m-1·s-1

νtturbulent kinematic viscosity,m2·s-1

σ surface tension of aqueous solutions,N·m-3

σctcritical surface tension of the packing material,N·m-3

σεturbulence model constants for diffusion of ε

χ local porosity of the random packings

Subscripts

e effective

G gas phase

I interfacial

in inlet

L liquid phase

Acknowledgments

The authors acknowledge the assistance by the staff in the State Key Laboratories for Chemical Engineering(Tianjin University).