Numerical study on the effect of stacking on the apparent reactivity in fixed bed reactors


                      Extrudates catalyst particles



Encadré par : CLIMENT Eric et ROLLAND Mathieu  (INP et IFPEN)

Most catalytic and petrochemical reactions are operated with fixed bed reactors. In these reactors, catalyst pellets are generally randomly stacked in a large cylindrical vessel and the reactants, usually gas and liquid, are flowing through the bed to react inside the catalyst pellets. The aim of this project is to perform numerical simulations to study the effect of stacking on the mass transfer and the apparent reactivity in those pellets .The used tool is COMSOL-Multiphysics 3.5a.




Context: .

Most catalytic and petrochemical reactions are operated with fixed bed reactors. In these reactors, catalyst pellets are generally randomly stacked in a large cylindrical vessel and the reactants, usually gas and liquid, are flowing through the bed to react inside the catalyst pellets. Catalyst pellets are made of a porous support onto which are deposited active materials. They are typically from 2 to 5mm in size with spherical, cylindrical or more complex shapes (tri-lobic or quadri-lobic extrudates).  More efficient processes can be designed based on a better comprehension of the physics ruling reactor performances in term of chemical response. A finer representation of the catalytic process allows smaller over-designs which will have fine economical results.


The objectives of this project are to:

■ Validate the Ranz-Marshall correlation

■ Find new correlation for cylindre

■ Study the effect of different stacked configurations on the outlet concentration

Numerical Method:

The used code is Comsol multiphysics version 3.5a through which we solve the  incompressible Navier-Stokes equation coupled with the convection diffusion equation , the calculation is Stationary .


The Incompreisble N-S equation:

$$ \frac{ \partial u}{\partial t}+u\cdot{\nabla{u}}-\nu{\nabla^{2}{u}}=-\nabla{w}+g $$

The convection diffusion equation:

$$ \frac{ \partial c}{\partial t}=D \nabla^{2}{c} -\vec{v}\cdot{\nabla}{c} $$

For stationary calculation :                        

$$\frac{\partial{c}}{\partial t}=0$$

The equation becomes:

$$ 0=\nabla{\cdot}({D\nabla{c}}) -\nabla{\cdot}(\vec{v}{c})+R $$

The dimensionless numbers controlling the problem are :

 Reynolds number:



Schmidt number: 

                               $$Sc=\frac{\nu}{D}=\frac{\mu}{\rho{D}} $$    

Sherwood number:


Péclet number :



$ \bullet   \nu$   is the kienamtic viscosity

$ \bullet   D $   is the Diffusion coefficient

$ \bullet   \mu $  is the Dymanic viscosity

$ \bullet   \rho $   is the Density of the fluid

In order to validate our numerical method we have started by 3D simulations of a spherical particle in a flow to study the mass transfer and validate the Ranz-Marshall corellation:

$$ Sh=2+0.552Re^{\frac{1}{2}}Sc^{\frac{1}{3}} $$

This relation permits to calculate the mass transfer coefficient by integrating the flux over the particle surface and evaluating Sherwood number .

For the N-S equation:The boundary conditions are:


►No slip at the particle surface and the container walls

►Velocity inlet at the intlet

►Pressure outlet at the outlet


For the convection diffusion equation:The boundary conditions are:


►C=1 at the particle surface 

►Insulation symertry at the container walls

►Intlet concentration zero

►outlet :Convective flux



Ranz-Marshall :

Validation in 3D:

In order to validate Ranz-Marshall correlation: $$Sh=2+0.552Re^{\frac{1}{2}}Sc^{\frac{1}{3}}$$

Numerous simulations were carried out for solving the Navier-Stockes equation coupled with the convection-diffusion equation around a spherical particle placed in a box, the sphere has radius equal to $2.{10^{-4}}$ and the box dimensions are 15r x 15r where r is radius of the sphere. The boundary conditions are similar to those explained in the introduction.

Geometry and Mesh:



                                 Geometry                                                                      Mesh

The taken geometry is a sphere of radius $2.{10^{-4}}$ m inside a cubic box of side 15r .The inlet is from down and the outlet is up .






                          Slice velocity profile                                 Slice concentration profile         

Validation in 2D:

Due to memory constrains we were not able to increase Sherwood number beyond ~6 in 3D as we need to refine the mesh when we increase Péclet number because the thickness of the boundary layer $ \delta$ is inversely proportional to Peclet number:


                                                                 $\delta $ ~ $\frac{1}{Pe^{\frac{3}{4}}}$

For this reason we also used 2D axisymetric model as the mesh size stays reasonable for higher values of Pe.

Geometry and Mesh:


                                   Geometry                                                                                 Mesh



                        Velocity profile                                                    Concentration profile


The velocity at the surface is always zero (No slip) but we have  seen that the boundary layer at the particle is getting thinner as we increase $Re$, it’s known that the thickness $\delta$ is inversly proportional to $Re$ and so when Re increases the thickness of the boundary layer decreases, the same applies for the diffusion, the layer thicknes decreases as we increase Péclet number.





In the above figure we compare the Ranz-Marshall correlation to the results obtained by 3D and 2D-axisymetric simulation. We see clearly that the 2D-axisymetric is closer to the correlation, this is because we were able to increase the domain size surrounding the particle in 2D as the mesh size stays finite whereas this was not possible in 3D because the increment of the mesh cells needs higher memory which was limited .


                   The variation of Sherwood number as function of Péclet number

We see that 2D simulations are better fitting the Ranz-Marshall  correlation and the range of Péclet is higher , this is because in 2D the mesh can be refined and can be solved by the available memory.

Cylindre Correlation :


Cylinder  Correlation:        

In this section we study the mass transfer around cylinder, for this reason we put a cylinder  of hight $H$ and diameter $D$ such that $H=D$ and we perform simulations for constant Schimit number $Sc=100$ and Reynold number $Re \in [0,1]$ . We vary the angle $\alpha$ of the cylindre as explained in the figure below and we calculate the flux for $\alpha \in [0,90]$ in order to see the effect of the angle on the mass transfer and in order to find a correlation similar to Ranz-Marshall correlation for the cylinder.



                                                     The angle $\alpha$ between the cylinder axis and the vertical




                           Geometry                                                                                          Mesh

As proposed by some papers , the correlation in convection-diffusion phenomena can be written according to the power law:

$$ Sh=A+BRe^nSc^m$$

Where the constants can be detrmined from the fitting of the data . From pure diffusion A can be determined because at velocity equal to zero there is no convection and $Re=0$ and so the second term vanishes , from pure diffusion we have obtained  A=2.48 . The other constants will be detrmined from the fitting .


                           Velocity profile-Re=0.4

                                                        Velocity profile-Re=0.4


                                                           Concentration profile-Re=0.4,Pe=40


                            The calculated Sherwood number for diffirent angles with the basic fitting  .

The above figure shows the Sherwood number variation for different angels $\alpha$ in a laminar flow , the flux decreases when $\alpha$ increases untill it reaches a minimum at $\alpha=90$. At $\alpha=0$ the calculated values have a basic fitting $Sh=Sh_0+Re^{0.4}Sc^{0.46}$ where $Sh_0=2.48$ and has been determined from pure diffusion.





Stacking Effect

In this section we  study the effect of stacking on the outlet concentration.

By evaluating the integral :

   $$ Sh=-\frac{1}{2{\pi{r}{\Delta{c}}}} \int{( \frac{\partial{c}}{\partial{x}}n_x+\frac{\partial{c}}{\partial{y}}n_y+\frac{\partial{c}}{\partial{z}}n_z} )dS$$

over the surface of the particle we calcuate Sherwood number and we use it to calculate the mass transfere coefficient from the fact that $ Sh=\frac{KL}{D}$ .

Studied geometries:

Due to memory constrains our study was limited by some specific geometries through which Comsol Multiphsics 3.5a can solve. We start with a very simple case through which we show the variation of the Sherwood number .



     (A)$\varepsilon=0.7$            (B)$\varepsilon=0.6$               (C)$\varepsilon=0.55$


In our study spheres where inserted inside the box where we consider that these spheres act as a catalyst in the packed bed reactor. We studied three cases where void fration has been taken as 0.7 in case(A), 0.6 in case (B), 0.55 in case(C) respectively. Free Mesh has been taken in all the three different geometry.

Concentration profile has been analyzed for all the different cases.  



                                                                                                         Configuration(A)-Concentration profile




                                               Configuration(B)-Concentration profile



                                                 Configuration(C)-Concentration profile

In the above figure we find that for the configuration (A) which has highest void fraction, concentration at the outlet of the box is maximum, which means that less catalysts are present inside the box and less reactants are consumed. For the configuration (C) which has the least void fraction, which means more catalysts are present inside the box, more reactants are consumed and less  outlet concentration is observed. For the Configuration (B) which has the value of  void fraction in between Case A and Case B, Value of the Concentration found at the outlet for configuration (B) is in between the Value of concentration for configuration (A) and configuration (B).

                 Outlet concentration as function of k for different geometries

Above figure represents the variation of the outlet concentration with the different reactivity constant for system having different void fraction. One having higher void fraction has more outet concentration and one with lower void fraction has lesser outer concentration.


The objective of our project was to study the effect of stacking on the hydrodynamic and the mass transfer inside a fixed bed reactor.

Ranz-marshall correlation has been validated for a single sphere inside the box for 2D as well as 3D. Outlet concentration profile has been studied in the case of stacking with different void fraction and the effect of stacking on the concentration has been studied.

Hydrodynamic and mass transfer with cylinder inside the box has been studied by changing the orientation of cylinder. Correlation of Sherwood with Reynolds and Schmidt for a single cylinder placed inside finite box has been obtained.

When the diffusivity of fluid is decreased peclet number increases and the boundary layer thickness $\delta$ decreases. For $D_e$  less than $10^{-6}$, we were not able to refine the mesh because when the number of cells increase the memory required by the simulation increases beyond the available limit and Comsol was giving "out of memory error" so we were not able to carry out our study for wide range of parameters as well as for complex geometry.The advance version of the comsol or may be some different code would be better for studying complex geometry where we can use wide range of physical parameters.