Bing Yuan,Junbo Xu*,Zaisha MaoYongqiang Zhang,Chao Yang,*

1College of Chemical Engineering,Sichuan University,Chengdu 610065,China

2CAS Key Laboratory of Green Process and Engineering,Institute of Process Engineering,Chinese Academy of Sciences,University of Chinese Academy of Sciences,Beijing 100190,China

3SINOPEC Research Institute of Petroleum Processing,Beijing 100083,China

Keywords:CFD Packed bed Model Simulation Transport

ABSTRACT The traditional fixed-bed reactor design is usually not suitable for the low tube-to-particle diameter ratios(N=D/d ＜8)where the local phenomena of channeling near the wall and backflow in the bed are dominant.The recent“solid particle”meshing method is too complicated for mesh generation,especially for non-spherical particles in large random packed beds,which seriously hinders its development.In this work,a novel high-fidelity mesh model is proposed for simulation of fixed bed reactors by combining the immersed boundary and adaptive meshing methods.This method is suitable for different shapes of particles,which ingeniously avoids handling the complex“contact point”problem.Several packed beds with two different shapes of particles are investigated with this model,and the local flow in the bed is simulated without geometrical simplification.The predicted pressure drop across the fixed bed and heat transfer of the single particle are in good agreement with the corresponding empirical relations.Compared with spherical particles,the packed bed packing with pentaphyllous particles has lower pressure drop and better heat/mass transfer performance,and it shows that this method can be used for the screening of particle shapes in a fixed bed.

1.Introduction

Fixed bed reactors are widely used in the chemical industry for catalytic reactions.The design of such reactors is often done with simplifications,e.g.,the assumptions of plug flow and equally distributed viodage over the whole fixed bed[1],and based on predictions by empirical correlations,such as the Ergun equation[2]for pressure drop and the Dixon equation[3]for voidage.However,these assumptions are too coarse for accurate modeling of fixed beds,especially for low tube-to-particle diameter ratios(2 ＜N ＜8,N=D/d),due to the presence of strong wall effects across the entire radius of the bed[4].

In recent decades,researchers began to consider the reproduction of the void structure in a fixed bed so as to describe exactly the fluid flow.Giese et al.[5]used Laser-Doppler-Velocimetry(LDV)to measure the velocity profile within the packed bed.In a cross-section,they measured about 5000 particles inside and outside a sampling volume.These values were averaged in the circumferential direction to obtain the fluid flow velocity distribution in the bed void along the radial direction.Similarly,nuclear magnetic resonance(NMR)imaging[6,7]and particle tracking methods[8]have also been used to observe the flow phenomena in fixed beds.These methods are subject to certain limitations:LDV requires the fluid to be refractive index matched with the transparent material of the packed bed wall,the flow rate of the NMR method must be very low,and the particle tracking method requires the tag to be observed and counted[9].

With the increase in computational capacity and speed,it is feasible to use computational techniques to simulate the fluid flow in fixed beds in detail.Lattice Boltzmann method has been used by many researchers[10-12]to simulate packed beds because it has a potentially high efficiency in computing somewhat larger scale of packing.However,the shortcomings in handling energy balance significantly affect the heat transfer simulation,which cannot be neglected[13].

As for body-fitted meshes for computational fluid dynamics(CFD)simulations of fixed bed reactors,the cells near particle-particle and particle-wall contact points are either highly skewed,which would lead to convergence problems during the calculation,or the mesh is highly refined in these regions,which increases the number of cells and as a direct consequence the computational time [1].There are four primary meshing methods,i.e.,“Gaps”,“Overlaps”,“Bridges”and“Caps”,to deal with contact points.The“Gaps”approach was used forCFD in early works[9,14],which avoids the contact point problem by shrinking the particles.The“Overlaps”meshing method refers to enlarging the particles,and it was put forward by Guardo et al.[15],in which the contact points become contact lines or faces,thereby eliminating the problem of high skewness of the cells near the contact point.The“Bridges”meshing method means to create connecting bridges at the particle-particle and wall-particle contact points[16].The“Caps”approach is achieved by cutting a plane at the contact point[1].These four methods can all be called“solid particle”meshing methods,and they can also be divided into two types of local meshing methods(“Bridges”and“Caps”)and global meshing methods(“Gaps”and“Overlaps”)[17,18].In global meshing methods,the particle diameter change is 0.5%-5%[4,9,19,20],which would drastically influence the flow patterns in the system and both the convective and conductive heat transfer mechanisms [4].Dixon et al.[21]pointed out that local meshing methods were better than global meshing methods for drag coefficient/pressure drop calculations,and only the Bridges model offered the flexibility for accurate computation of particle-particle and particle-wall heat transfer with the expense of increasing complexity and computing resources.However,in local meshing methods,the positions and orientations of each contact point need to be documented and then the contact points can be processed locally.It is not easy to be used in randomly packed beds with a large number of particles,especially in the case of non-spherical particles.

Table 1 Particle properties

The immersed boundary method(IBM)developed by Peskin[22]was firstly used to simulate cardiac mechanics and blood flow.This method is performed on a Cartesian mesh,which does not conform to the geometry and feedback forcing to represent the effect of boundaries.Compared with the body-fitted mesh method,IBM has the advantages of saving memory and CPU and easy generation of meshes without sacrificing accuracy[23,24].Based on these characteristics,the IBM is suitable for simulating the flow field with complex shape structure and deals with various dynamic boundary problems.In recent years,it has been widely used to simulate the free-fall motion of objects[25],blood flow in the heart[22],and ocean flow[26].

In the present study,the commercial CFD tools combine userdefined IBM functions to create and simulate the packed bed.Theapplication of IBM avoids the contact point problem,and the adaptive meshing method ensures that the total number of cells is not too large while refines the cell near the particle surface.This method is suitable for different shapes of particles and uses an orthogonal hexahedral mesh.The results obtained by this approach are in good agreement with literatures in respect to viodage distribution and pressure drop and single spherical particle heat transfer.A comparison of the flow and heat/mass transfer performance of the packed bed of the two particles shows that this method can be used for the screening of particle shapes in a fixed bed.

Table 2 Geometrical parameters of different particles

2.Mesh and CFD Modeling

2.1.Generation of random packing and meshing

During the random packing of a fixed bed,the particles are usually poured into a confining cylindrical tube,and move randomly to the column bottom under the action of gravity.DEM is an explicit numerical method based on contact mechanics that can be used to simulate the dynamic and static behavior of particles [27].In this study,the in house discrete element method code was used to simulate the random packing process of particles.The properties of the particles include density,shear modulus,Poisson's ratio,etc.,which affect the void ratio of the bed and the packing structure of the particles are shown in Table 1.The values of density,shear modulus,Poisson's ratio and coefficient of restitution are chosen referring to values of typical rock just for reality in the magnitude.Fig.1 shows the schematics of the two types of particles used in this work.Spherical and pentaphyllous particles are shown in Fig.1(a)and (b),respectively.Fig.1(c)shows a schematic crosssection of the pentaphyllous particle.The detailed geometries of these particles are shown in Table 2,and the two particles have equal equivalent spherical diameters.

Fig.2 shows the schematic of several packed beds numerically generated by DEM.Fig.2(a),(b)and(c)show the packed beds with different tube-to-particle diameter ratios formed by spherical particles.Fig.2(d)is the packed beds of random packing with pentaphyllous particles.The tube inner diameter(D)and the packing height(H)are marked in Fig.2.Table 3 shows the parameters of each packed bed sketched in Fig.2.The tube-to-particle diameter ratios of the packed beds and the number of particles in Fig.2(c)and(d)are equal.

Fig.1.Schematic diagram of the geometry of the particles.

Fig.2.Schematic of different packed beds:(a)random packing with spherical particles at N=D/d=2.667,(b)random packing with spherical particles at N=D/d=3,(c)random packing with spherical particles at N=D/d=4,(d)random packing with pentaphyllous particles at D/dESD=4.

Table 3 Packed bed parameters

A schematic of the method is shown in Fig.3(a),the steps of this method are described as follows:(1)An initial non-dense hexahedral mesh is generated;(2)DEM is employed to obtain the stacking structure of a packed bed and distinguish the solid and fluid regions;(3)the meshes are adaptively subdivided at the immersed boundary region;(4)the fixed bed is simulate at given boundary conditions.Taking spherical particles as an example,in the Cartesian coordinate system,the definition of the solid regions and the fluid regions can be summarized as follows.For the solid regions,the cell belongs to the solid region when the distance (li)between the cell center and the center of the spherical particles must satisfy:

where(x,y,z)is the center of the cell to be judged within the domain,Npis the number of particles and(xi,yi,zi)is center of the ith particle.For the fluid regions,it must satisfy:

Using the meshing method of this study,we can avoid the processing of numerous critical points,i.e.,particle-particle and particle-wall contact points,and the whole computing domain can be discretized with structured hexahedral meshes.In addition to spherical particles,this meshing method can also be applied to cylindrical,toroidal ones and so on.In this paper,spherical and pentaphyllous particles are studied.

Fig.3.Schematic diagram of the method(a)and the adaptive mesh(b).

Isometric Cartesian mesh is often used in IBM,which wastes computing resources to some extent.In order to reduce the computational time,an adaptive mesh method with the immersed boundary has been employed.The schematic of adaptive orthogonal hexahedral mesh in this work is shown in Fig.3.The adaptive encryption of the mesh is implemented by the adapt module in the commercial software Fluent 15.0.It is to be noted that the iso-values of the y+(defined with UDF and described in the next section)near the interface must satisfy y+-criterion even when the Reynolds number (Rep)is as large as 15000.The Repis defined as

where u0is the inlet velocity.Compared with the uniform mesh,the application of the adaptive mesh method reduces the computational complexity by about 95%.

2.2.Solver method

As shown in Fig.3,particle surfaces divide the computational domain into a solid regions Rsand fluid regions Rf.In fluid regions,the flow equations(Navier-Stokes equations)and heat transfer equations for incompressible fluid of constant physical properties are as follows:

here xiare the Cartesian coordinates and viare the velocity components.p is the static pressure,t is the time,ρ and μ are the density and dynamic viscosity,respectively.and kTare the constant pressure heat capacity and thermal conductivity,respectively.In solid regions,heat transfer exists only in the form of heat conduction.

As the large amount of computation,the direct numerical simulation method is currently limited to the calculation of fluid mechanics at low Reynolds numbers.Reynolds-averaged Navier-Stokes(RANS)theory assumes that the flow feild variables in turbulent folw are decomposed intoand vi′components.Assuming that the turbulent Reynolds stress is proportional to the stress,it calculates the average components rather than the turbulence fluctuations at all scales,which reduces the computational complexity.When it is used in Eq.(5),the result is

Velocity and other solution variables are now expressed as Reynolds averaged values,and the effect of turbulence is expressed as Reynolds stress,The commercial software Fluent 15.0 combined with user-defined packed bed structure was used for the simulations.Due to the use of a non-body-fitted mesh,user-defined solid wall boundary conditions(no-slip boundary conditions)are required.For turbulent flow,the wall function method is applied to the solid particle surface,which is achieved by adding a source term to the equation[28].Similar to the ghost-cell immersed boundary method (GCIBM)[26],the cell near the solid wall is defined as the boundary region,here we use the standard wall functions [28]to associate the fluid region with the solid region.The expression of the boundary forcingis as follows:

where k and Acellare the turbulent kinetic energy of the control volume and the area at the wall surface,respectively.Δl is the distance from the control volume to the solid wall surface.The constant Cμis equal to 0.0845 in the RNG k-ε model[29].The dimensionless parameter vpis the velocity parallel to the boundary,and the velocity perpendicular to the wall is defined as zero.Γ is the boundary layer region near the solid boundary.v+and y+are the dimensionless parameter,y+=11.63 is usually viewed as the boundary between the viscous bottom layer and the logarithmic law layer[28].The near-wall heat transfer source term STis as follows:

For all geometries,the whole process,including the formation of a fixed bed,meshing and flow simulation,only takes 3 to 8 h depending on the number of cells,which scales with the number of particles and the mesh accuracy(Δl/d).All simulations were done in parallel on 60 threads(Intel Xeon processor,2.1 GHz,31.3 GB RAM).

3.Results and Discussion

3.1.Viodage and pressure drop

The viodage of the packed bed is a very important indicator of the bed.Dixon et al.[21]obtained a rule of thumb that an increment of 1%in the viodage would give rise to a corresponding increment approximately of 3% in pressure drop.Through summarizing the literature data,de Klerk[30]introduced a radial distribution equation of viodage in a packed bed of spherical particles,

expressed separately for the near wall region and the core region,where Z=(D/2 ?r)/d is the normalized wall distance and εbis the infinite diameter average bed viodage.The maximum relative error of this equation is 20%and the error increases with decreasing N values.The local viodage ε(r)can be calculated by

Fig.4.Radial average distributions of viodage for different D/d ratios:(a)1#,N=2.667,(b)2#,N=3 and(c)3#,N=4,(d)4#,N=4.

where A(r)solidrepresents the solid shell area at r.As shown in Fig.4(a)-(c),the radial average viodage of the packed bed of spherical particles with different N values is compared with the empirical equation by de Klerk [30].The viodage distribution with spherical particles agrees well with the trend of Eq.(10),and the error mainly occurs near the axis channel of the packed bed.The formation of such channels is due to the combination of the effect of restrictions by the wall and the high coefficient of friction,which prevents the particles from sliding towards the channel [1].Fig.4(d)is the radial average viodage distributions of packed beds of pentaphyllous particles.(The normalized wall distance is defined as ZESD=(D/2 ?r)/dESD.)For all shapes of particles,the contact between the particles and the wall surface is a point or a line,causing the viodage at the wall surface to approach 1.

Fig.5.Velocity magnitude distributions of the packed beds at inlet velocity u0=5 m·s?1(Rep=3422.9):(a)tube 3#,(b)tube 4#.

Fig.6.Volume fraction of different velocity magnitude in the packed beds at inlet velocity u0=5 m·s?1(Rep=3422.9).

Fig.5(a)and(b)shows the velocity distributions in spherical and pentaphyllous particles packed beds at Reynolds number of 3422.9,respectively.Due to the variety of fluid flow channels in a fixed bed,as shown in Fig.5(a)and(b),the velocity magnitude of the fluid varies from 0 to 40 m·s?1.In addition to the solid region,there are zones with stagnant points at the wake of the particles and at the end of the fixed bed.Fig.6 calculates the volume fraction of fluid elements at different velocities in the fluid region at a Reynolds number of 3422.9.The results show that the packed bed of pentaphyllous particle has the larger volume fraction of the low velocity zone and the less volume fraction of the high velocity zone.This is because the packed bed packing with pentaphyllous particles has larger viodage flow channel.

One of the most important indicators for fixed beds is the bed pressure drop.Eisfeld and Schnitzlein[31]compared the pressure drop predictions in the literature with the relevant experimental data and found that the method of correcting the Ergun equation for the wall effect by Reichelt[32]is the most promising.The pressure drop correlation[31]is

where K1=154,k1=1.15 and k2=0.87 for the spherical particle packed beds,and for other shapes of particles K1=155,k1=1.42 and k2=0.83.Fig.7(a)compares the pressure drop of spherical particle packed bed of different tube to particle diameter ratios to Eq.(12).It shows that the simulation results using this method agree well with Eq.(12)over a wide range of Reynolds numbers(50-15000).The pressure drop per unit length increases with the increase of the tube to particle diameter ratios,which is due to the weakening of the wall limitation and the decrease of the voidage.Fig.7(b)shows the relationship between the pressure drop and the Reynolds number in a packed bed of particles of different shapes.Compared with the results of Eisfeld and Schnitzlein[31],the maximum error of spherical particles is 14%,while the maximum error of pentaphyllous particles is 42%.It shows that in Eq.(12),uniform parameters are used to predict the pressure drop of packed beds of differently shaped particles,and the accuracy of the results will be insufficient.Although the 4#packed bed packed with pentaphyllous particles has a higher filling height,the pressure drop is smaller than that of the 3# packed bed because of its larger voidage.

3.2.Heat transfer

The performance of a fixed bed reactor is affected by the mass/heat transfer process between the catalyst and the fluid.By the similarity of mass transfer and heat transfer processes and the assumption that the surface of the sphere is a constant Dirichlet boundary,heat and mass transfer do not need to be studied separately.The heat transfer process of a single homothermal spherical particle in a square tube with a=6d under different Reynolds numbers is simulated.With these simulation results,the heat transfer coefficient(h)could be determined:

Fig.7.Pressure drops across packed beds over a large Reynolds number range:(a)tubes 1#,2#and 3#packed with spherical particles,(b)pressure drop comparison of different packed beds with differently shaped particles.

Fig.8.Comparison of single particle heat transfer Nusselt numbers over a large Reynolds number range.

Fig.9.Heat transfer of packed beds 1#,2#and 3#:(a)average Nusselt number of the packed bed under different Reynolds numbers;(b)the outlet temperature of packed beds.

where Tp,T∞,Tinand Toutare the particle surface temperature,the infinity fluid temperature,the fluid inlet temperature and the fluid outlet temperature,respectively.Fig.8 shows the comparison of the Nusselt number of a single spherical particle with literatures[33-35]at different Reynolds numbers.It shows good agreement between the simulations using our method and literatures[33-35]in a large range of Reynolds numbers(6-14000).

Similarly,the heat transfer of packed bed was also studied in this work.The particle temperature of the packed bed is defined to be constant at 400 K and inlet temperature of the fluid is 300 K.Different from single particle heat transfer,the average heat transfer coefficient(have)of a packed bed can be calculated by the following equation:

where ΔTlnis the logarithmic mean temperature difference.As recommended by Whitaker[36],the expression of the average Nusselt number of the packed bed is

As shown in Fig.9(a),the Nuavesimulated by this work at different Reynolds numbers is compared with that by Whitaker[36].It can be seen that the change of Nuavewith Repis consistent with that by Whitaker[36],but the Nuaveobtained in this work is smaller.This is because Whitaker's fitting data do not contain low N cases.In the case of lower N,the influence of the tube wall cannot be ignored,which will affect the heat exchange between the fluid and the particles.Similar phenomena can also be seen in the work of Romkes[37].Fig.9(b)shows that the outlet temperature of the packed bed decreases as the Repand N increase.It is because the lower N value has a larger fill height and a longer residence time.

Fig.10 shows the fixed bed heat transfer of differently shaped particles in tubes of the same inner diameter.Fig.10(a)shows that the average Nusselt number Nuaveof the packed bed gradually increases with increasing Rep.The Nuaveof the spherical particles is larger than that of the pentaphyllous particles at the same Reynolds number.However,as shown in Fig.10(b),at the same Reynolds number,the outlet temperature of the packed bed of pentaphyllous particles is higher.It is because the packed bed of pentaphyllous particles has a larger surface area.The sub-graph(I)in Fig.10(a)shows the temperature distributions in the packed beds of 3#and 4#,respectively,when the Reynolds number is 684.59.The larger outlet temperature means that pentaphyllous particles have a better heat/mass transfer performance between the fluid and the catalyst particles.Compared with spherical particle,the packed bed packing with pentaphyllous particles has lower pressure drop and better heat/mass transfer performance.This shows that the method can be used to select a more suitable particle shape for the packed bed.In addition,this high-fidelity mesh model can easily define the mass transfer and reaction in the particle region,which can be used to simulate the adsorption and reaction processes in a fixed bed.

4.Conclusions

In the present work,a novel high-fidelity mesh model is proposed for the study of packed beds.The DEM is used to obtain the coordinates and orientations of the particles,and the particle regions are defined according to a user-defined subroutine.In contrast to known methods,this high-fidelity mesh model avoids the particle-particle and particle-wall“contact points”problems.The application of IBM and adaptive meshing method ensures that the total number of cells is not too large while refines the cell near the particle surface.The results obtained by this method are in good agreement with the predicted viodage distribution and pressure drop of a packed bed packed with different shapes of particles.The heat transfer results of single spherical particles are also verified with literatures[33-35].Because the effect of the tube wall cannot be ignored under low N,the Nuaveof the packed bed is smaller than that of Whitaker[36].Compared with spherical particle,the packed bed packing with pentaphyllous particles has lower pressure drop and better heat/mass transfer performance between the fluid and the particles,and it shows that this method can be used to select a more suitable particle shape for the packed bed.In addition,this high-fidelity mesh model can easily define the mass transfer and reaction in the particle region,and in the future work we will simulate the adsorption and reaction processes in a packed bed by this method.

Fig.10.Heat transfer of packed beds 3#and 4#with differently shaped particles:(a)average Nusselt number of the packed bed under different Reynolds numbers;(b)the outlet temperature of packed beds,sub-graph(I)is the temperature distribution of the packed bed at Rep=684.6(u0=1 m·s?1).

Nomenclature

a side length of the square pipe,m

apthe specific surface area particle,m?1

D diameter of the cylindrical tube,m

d diameter of the spherical particle,m

H length of the cylindrical tube,m

h length of the pentaphyllous particle,m

hTheat transfer coefficient,W·(m2·K)?1

k thermal conductivity of fluid,W·(m·K)?1

l distance from the center(central axis)of a particle,m

N tube to particle diameter ratio

Npnumber of particles

Nu Nusselt number

NuaveThe average Nusselt number of packed bed

Repparticle Reynolds number

r radius of a spherical particle,m

T temperature,K

u0superficial velocity,m·s?1

v velocity,m·s?1

ε viodage

μ viscosity,Pa·s

ρ density,kg·m?3