Yanni Chi,Rui Zhang,Xianghai Meng,Jian Xu,Wei Du,Haiyan Liu,Zhichang Liu

State Key Laboratory of Heavy Oil Processing,China University of Petroleum-Beijing,Beijing 102249,China

Keywords: Liquid-liquid cyclone reactor Ionic liquid,separation Droplet deformation Droplet breakage Computational fluid dynamics

ABSTRACT A liquid-liquid cyclone reactor (LLCR) was designed to achieve mixing-reaction-separation integration during isobutane alkylation catalyzed by ionic liquids.However,studies of the droplets deformation and breakage in the kind of reactors are lacking.In this work,the research studied the velocity distribution,pressure field,and turbulent field to investigate the flow pattern and the main energy loss location in the LLCR through the computational fluid dynamics(CFD)method.The simulation results were verified by experiemnts to prove the correctness of the model.Then the deformation and breakage process of droplets,and the influencing factors of droplets breakage were studied by remodeling which was based on the tangential velocity distribution result of the three dimensional model.The three dimensional simulation results clearly showed that the pressure of the LLCR was mainly concentrated in the cone section and fluid turbulent motion was the most intense near the lateral wall.The reconstruct the results of the two dimensional model clearly showed that the deformation and breakage location of droplets were mainly occurred in the velocity boundary layer,while it was difficult to break in the mainstream region.In addition,low surface tension and high Weber number had a positive effect on droplet breakage.

1.Introduction

Alkylation of isobutane with butenes catalyzed by ionic liquids(ILA) occurs on the liquid-liquid heterogeneous surface or in the thin layer near the surface [1–3].The reaction rate is very fast and there are many side reactions including polymerization,disproportionation,and degradation reaction.Therefore,the process requires good liquid-liquid mixing to ensure the selectivity of the main reaction.On the other hand,rapid separation of the products from the ionic liquid is required to improve the quality of alkylate thanks to the low degree of side reactions [4].Thus,ILA requires both fast mixing and quick separation of ionic liquid and hydrocarbon phases.Liuet al.[5]and Zhouet al.[6]applied the Strator reactor and a novel liquid-liquid loop reactor for ILA and studied the flow behavior,mixing performance,and separation efficiency under different operating conditions to optimize the reactor geometry.However,the residence time of those reactors is too long for ILA.Therefore,it is necessary to study and develop a new reactor suitable for ILA.

Aiming at the rapid mixing and quick separation,our research group developed a liquid-liquid cyclone reactor (LLCR) [7,8]for ILA,as is shown in Fig.1.The isobutane and butene feed was injected through symmetric tangential inlet 2 and inlet 3,and the ionic liquid flowed axially from the heavy phase inlet 1 into the reactor.Four guided vanes (24 mm length,11° outlet angle)lay below the inlet 1.The ionic liquid gradually entered the highvelocity rotation owing to the action of the guided vane.Then,the centrifugal force caused the hydrocarbon phase to move toward the reactor center,which was mixed and separated with the heavy phases[9].Finally,the light phase exited from the overflow and the heavy phase from the underflow.Therefore,the LLCR had the advantages of mixing-reaction-separation integration,short contact time between the hydrocarbon phase and ionic liquid.Zhanget al.[10]concluded that the velocity of the dispersed phase and the relative velocity between the two phases were the main factors affecting the mixing and separation in the LLCR.Duanet al.[11]combined the electrical resistance tomography and numerical simulation to conclude that the increase of the recirculation eddy size was beneficial to the increase of the mixing performance and the decrease of the separation efficiency,while the tangential velocity had the opposite effect.The above research showed that the mixing and separation performance of the LLCR was affected by many factors.However,there was a lack of research on the energy and turbulence field in the LLCR.

In recent years,many scholars have done a lot of researches on the velocity field,the concentration field,and the pressure field of hydrocyclone [12–14].But reports on the phenomenon of droplet breakage are scarce.At present,there are many reports about droplet breakage in other fields.For example,scholars [15–17]used the theoretical analysis and experimental research to conduct a detailed study on droplet size distribution in a stirred tank,which continuously improved the frequency function and the criterion of droplet breakage.In the field of spray atomization,the researchers[18,19]combined theoretical exploration,experimental research,and simulation to carefully analyze the droplet size and velocity distribution.Some workers[20–23]studied the mechanism of droplet breakage.Their researches indicate that the droplet breakage is mainly caused by the shear force and viscous shear force.The shear force is generated by the turbulent flow.The viscous shear force is produced by the time-averaged velocity gradient between the dispersed and the continuous phases.However,the size of the droplets studied by researchers is generally larger at the millimeter level.In practical applications,the droplet diameter in the LLCR can be as small as the micrometer level[24].Besides,the nature of the two-phase system in LLCR is different from that of previous studies.Recently,Zhanget al.[25]studied the droplet size distribution and the Sauter average diameter near the wall under different operating parameters in the LLCR,and the results showed that the droplet diameter was 70–110 μm under the optimal conditions.However,the study only measured the local droplet size distribution.Therefore,the phenomenon of droplet breakage needs further research.

In this work,the two-phase pressure field and turbulence field in the LLCR were simulated to investigate the main energy loss location and the turbulence situation.After that,according to the research results of the two-phase flow,the droplet breakage behaviors at the wall and the factors affecting droplet breakage were studied by re-modeled to provide data support for the design and optimization of LLCR.

2.Numerical Methods

Fig.1.Schematic diagram of the LLCR structure.

Computational fluid dynamics (CFD) has been widely used by researchers and process engineers to analyze the flow performance in equipment.Two types of CFD models are often employed in the multiphase flow:the Eulerian-Lagrangian model and the Eulerian-Eulerian model.The method depends on the physical properties and flow patterns of the two phases.The Eulerian-Lagrangian approach can be used to establish the separation model between primary and dispersed phases,but a volume fraction of the disperse phase is required for less than 10% [26–28]in the current case.The Eulerian-Eulerian approach is used to simulate with high loading (more than 10%) of the dispersed phase.There are three submodels in the Eulerian-Eulerian model,the Euler model,the hybrid model,and the VOF (Volume of Fluid) model.The VOF model is a surface tracking method that is under a fixed Euler grid.Two or three incompatible fluids are simulated by solving a separate momentum equation and processing the volume fraction of each fluid passing through the region.VOF uses transportable scalars to implicitly track interface information and does not need additional techniques to deal with interface topological changes.There were two reasons for choosing the VOF model.One reason was that both phases in the LLCR were mutually incompatible liquids.The other reason was that ILA reactions occurred in or near the interface region between the ionic liquid and hydrocarbon phases and the phase interface required to be traced.

2.1.Turbulence model

The turbulence closure model is an important part in the fluid dynamics calculation in the LLCR.Previous numerical simulation shows that the Reynolds stress model(RSM)is suitable for swirling flow in hydrocyclones,since it can predict anisotropic turbulence[29,30].The RSM model strictly considers the rapid changes of streamline curvature,vortex,rotation,and tension,which solves the Reynolds pressure problem and dissipation rate and has a high precision for the simulation of the complex flow field.The large eddy simulation (LES) also provides a good prediction of turbulence encountered in the hydrocyclones [31,32],however,it requires higher computational configuration than RSM.Therefore,the RSM model was chosen to simulate the momentum transport induced by turbulence.

The interphase forces mainly include lift force,virtual mass force,and drag force.Drewet al.[33]simulated the oil-water emulsion flow,and found that the drag force had the most important influence on the flow distribution.The main interaction force between two phases in the LLCR was the drag force and the Schiller-Naumann model was generally applicable to liquidliquid multiphase flow.The drag force between the continuous and dispersed phases can be expressed as Eq.(1)andReis the relative Reynolds number.

2.2.Geometry and meshes

The geometry and sizes of the LLCR are displayed in Fig.1 and Table 1.The cylindrical section of the LLCR was equipped with a guide vane with a spiral structure to enhance the turbulence flow of the heavy phase.Fig.1 demonstrates the structure of the guide vane.

Table 1Size distribution in the cyclone reactor

ICEM CFD 14.5 was used for the creation of geometry and meshing as the pre-processing tool.In consideration of the complex structure of the intermediate vanes,the geometry was divided into three parts.Structured grids were used for the upper cylinder,lower sections,and the cones,while unstructured grids were used for the intermediate vanes.The LLCR grid model is shown in Fig.2.

2.3.Solution methodology

The governing equations were discretized with the QUICK scheme.The coupling between velocity and pressure was dealt with SIMPLE algorithm.The solution was considered acceptable when the scaled residual of each term was below 1×10-4.The time step was set to 1×10-4s.The calculation reached stability when the time reached 5 s.

2.4.Boundary conditions

Glycerin-water mixture and kerosene were selected to replace the heavy phase (ionic liquid) and the light phase (hydrocarbons)since they had similar density and viscosity,as shown in Table 2.Also,they were safer and less toxic as experimental materials.Both heavy and light phase inlets were set as velocity inlets.The total flow rate was 1.2 m3·h-1and the volume ratio of the heavy to the light phase was 2.According to the structural parameters of the LLCR,the inlet boundary conditions of the light and heavy phases are listed in Table 3.Since the length of the overflow and the underflow in the LLCR were more than 10 times the diameter of the pipe,it could be considered that the flow at the outlet was fully developed.Therefore,both the overflow and the underflow were set to outflow.

Table 2Physical properties of relevant fluids

Table 3Calculation parameter values of two phases in the LLCR

3.Results and Discussion

3.1.Grid dependence

Fig.3(a)and(b)shows the velocity at the cone section and Fig.3(c)shows the oil content of the overflow with different calculation times.There was a little difference in the velocity between 240,000 and 440,000 grids,which was slightly higher than 800,000 grids.In addition,440,000 grids had the highest oil content of the overflow and were relatively stable,because every small grid of 440,000 grids was cubing,which geometry was neat and the size was suitable.However,every small grid of 240,000 grids was big,which led to flow fluctuation.800,000 grids were fine,which meant the big deformation of a small grid,leading to computation instability.Therefore,the 440,000 meshes were chosen for the subsequent study to ensure reliability.

3.2.Model validation

It is necessary to validate the model before the analysis of the numerical simulation.The separation efficiency (η) was defined as the ratio of light phase flow rate at the overflow to that of the light phase inlet.The values of η from both experiments and simulations were compared to prove that the model was reasonable.The η is calculated byqol/qil.

Whereqolandqilare the mass flow rate of light phase at the overflow and the inlet,respectively.Fig.4 illustrates the η values from both experiment and simulation,and the maximum and average relative error between experiment and simulation was 6.53%and 4.58%,indicating that the results of the calculation model were acceptable.

Fig.2.Mesh for the simulation of the LLCR.

3.3.Flow characteristics

3.3.1.Turbulent flow distribution

Fig.5 shows turbulent kinetic energy for the axial position of the LLCR.The turbulent kinetic energy of the cross-section in the LLCR exhibited axis symmetry and it was intense near the wall,which indicated that droplet breakage occurred mainly near the wall.The turbulent kinetic energy decreased rapidly with the decrease of the radius.As the fluid flowed downward,the turbulent energy near the sidewall in the mixing and tail sections decreased,while that of the cone section showed an increasing tendency.It meant that the energy loss increased in the cone section and decreased in the mixing and the tail sections.The simulation had similar rules comparing with the research of Penget al.[34]on the turbulent field of cyclone separator.

Fig.6 shows the turbulent intensity for the axial position of the LLCR.The distribution of turbulence intensity in the three parts was similar to the turbulent kinetic energy.

Fig.7 shows the turbulent viscosity for the axial position of the LLCR.The turbulent viscosity near the wall and in the axial center was relatively large in the mixing section,while it was small in the cone section.Besides,the tail section showed in the opposite direction from the mixing section.Comparing the three parts,the turbulent viscosity was the largest in the cone.The turbulent viscosity was proportional to the turbulent Reynolds number,therefore,the turbulent phenomenon was more serious in the cone section.

3.3.2.Pressure field distribution

Fig.3.Velocity distributions under different grid quantities at(a)z=30 mm and(b)z=-5 mm,(c) Variation of overflow oil content with calculation time under different grid quantity.

The concentration of the light phase in the upper part of the tail section was much higher than that of in the lower part,as shown in Fig.8(a).Therefore,the tail section was divided into two parts.Finally,the LLCR was divided into four regions in the axial direction: I,II,III,and IV,which were the mixing,separation,circular pipe,and tailpipe section,respectively.Fig.8(b) shows the static pressure distribution for the axial position of the LLCR.The static pressure was relatively large in region I because the sudden expansion of the section caused great resistance when the heavy phase entered the LLCR cylinder from the entrance.The light and heavy phases were mixed and the descending height of liquid was relatively small in region II.The static pressure in regions III and IV decreased linearly owing to stable flow patterns and regular geometry.

Fig.4.Separation efficiency of the experiments and simulations.

Fig.8(c)shows the dynamical pressure distribution for the axial position of the LLCR.The dynamic pressure was proportional to the square of the fluid velocity,so the large dynamic pressure led to the large velocity.The dynamic pressure was drastically reduced in region I,because the velocity of the two phases began to reduce drastically due to the collision and the fluid motion resistance.It indicated that the heavy phase accelerated by the guide vane was fastly mixed with the light phase.The dynamic pressure decreased slowly at first and then increased slowly in region II.At this time,the heavy phase and the light phase were gradually separated.The light phase converged to the center to form the inner swirling upward,while the heavy phase gathered near the sidewall to form the outer swirling downward.The dynamic pressures in regions III and IV remained almost the same,indicating that the light and heavy phases were almost completely separated.Only the heavy phase flowed in the tail section and its velocity remained constant.

Fig.8(d)shows the total pressure distribution for the axial position of the LLCR.The total pressure was equal to the sum of the static and dynamic pressures.The pressure was the largest in region I owing to the intensification of secondary flow in the guide vane.The change of pressure was relatively stable in region II,indicating that the forces separating the heavy and light phases were mainly provided by centrifugal forces.The pressure in regions III and IV decreased linearly because the steady flow pattern of the light and heavy two-phase flow showed a straight downward flow.

The pressure drop was represented by the resistance coefficient to confirm the energy loss of each part of the LLCR.

where ξ,Δp,ρa(bǔ)ndv are the resistance coefficient,the pressure drop,density of mixing phases and flow rate,respectively.

Table 4 shows the resistance coefficient in each zone of the LLCR.The resistance coefficient was the largest in region II,which meant the pressure loss was mainly located in the separation section.According to the analysis of the turbulence field,the turbulence was the strongest in the cone section,which indicated that the strong swirling flow was one of the main causes of pressure drop.

Fig.5.Turbulent kinetic energy in the LLCR,(a)mixing section,(b)cone section,(c)tail section.

3.4.Change of droplet morphology

The droplet breakage in the LLCR generally consists of processes such as droplets deformation,droplets breakage,and small droplets coalescence.To simulate the droplets deformation and fragmentation in the swirling field,the following simplifications and assumptions were made.(1) The effect of the tangential velocity on the two-dimensional droplets was considered only.(2)Droplets were considered to be broken if the concentration of droplets in the heavy phase was less than 1%.(3) The effect of gravity on the droplets was ignored.Based on the above assumptions,the geometrical model was simplified to a two-dimensional structure of 180° curved,as shown in Fig.9(a).

Fig.6.Turbulent intensity in the LLCR,(a) mixing section,(b) cone section,(c) tail section.

To make the inlet velocity of the two-dimensional structure consistent with the tangential velocity distribution in the LLCR,the tangential velocity distribution along the radial direction from 9 to 15 mm atz=30 mm was selected.The inlet velocity was implemented through a user-defined UDF.Then three droplets A,B,and C with a diameter of 100 μm were placed in the computational domain filled with the glycerin-water mixture,as shown in Fig.9(b).

The broken time and position of three droplets are listed in Table 5.Comparing with droplets A and B,the broken time the droplet C was the shortest.Based on the analysis of the turbulence field,the turbulent kinetic energy near the wall surface was relatively large to make the time and space scale of droplet deformation faster,which was conducive to the occurrence of breakage.

Table 4Pressure drop,flow rate,density and resistance coefficient in each area

Table 5Broken time and broken position of three droplet

Fig.7.Turbulent viscosity in the LLCR,(a) mixing section,(b) cone section,(c) tail section.

The morphological changes of droplet C at different times are shown in Fig.10.At first,the droplet was round,then the stretching length and the spreading radius increased with the prolonged time.As a result,a long droplet with two cores was formed.The droplet was broken into two small droplets when the droplet length reached a critical value.Then the two small droplets continued to move,stretched and broken into three or more droplets.The two droplets in the middle were coalesced due to the small distance,then separated,continuously deformed and broken.Finally,droplet C was broken into 8 small droplets.Comparing with the droplet breakage in the jet field in the literature by Raoet al.[35],the deformation and breakage time was shorter owing to the shear force produced by the tangential velocity gradient.However,the droplet breakage also occurred at a violent position in the turbulence field.

3.5.Analysis of factors affecting droplet breaking

3.5.1.Effect of the droplet size

Droplets with diameters of 40,60,80,100,120,140,160,and 180 μm were examined to study the effect of droplet size on breakage.Fig.11 shows the relationship between diameter and broken angle of droplets.Droplet C with a size ranging from 40 to 180 μm was broken.The broken angle first decreased then changed a little with the increase of size,which indicated that the time from the inlet boundary to the droplet breakage became shorter.Droplet breakage required the small surface energy per unit area,which led to the decrease of the droplet stability.Droplet B was not broken when the size was less than 100 μm.Because small droplets followed the heavy phase and flowed out of the outlet quickly.The broken angle first decreased then increased under the size above 100 μm,since droplets became more and more stable in the middle flow.

3.5.2.Effect of the droplet slip velocity

The initial velocity of droplets with a diameter of 40 μm was set as 0,1,1.5,2,and 2.5 m·s-1to investigate the effect of relative slip velocity on breakage.The broken angle decreased with the increase of slip velocity,as shown in Fig.12,since the droplet got more energy than the surface energy.The effect of the heavy phase on the inertial force of droplet led to the time required for breakage short.

3.5.3.Effect of the surface tension coefficient

The function of surface tension was to maintain droplet morphology.Fig.13 shows the relationship between the surface tension coefficient and a broken angle.The droplet was broken when the surface tension coefficient was not above 0.016 N·m-1.It indicated that the critical surface force per unit length to break the droplet was 0.016 N·m-1.The sum of shear stress and flow resistance of the droplet surface that the heavy phase acted on was greater than the surface tension.Owing to surface tension preventing droplet fragmentation,the droplets were directly spread into a film on the arc-shaped wall surface when the surface tension coefficient reached a critical value.This rule was similar to the conclusions of Galinatet al.[36].

3.5.4.Effect of Weber number

The broken angle decreased with the increase of the Weber number,as shown in Fig.14.Weber number represented the ratio of inertial force to surface tension,so the larger Weber number meant the larger inertial force.Inertia force and surface tension were competing for relations to determine the droplets breakage.The variation of inertial force greatly exceeded surface tension,when other factors remained constant,so it was inferred that the inertial force played a dominant role in determining the droplet breakage.

4.Conclusions

(1) The static pressure distribution conformed with the distribution rules of the eddy field.The pressure decreased from top to bottom along the axis of the LLCR.The pressure loss in each part of the LLCR was mainly concentrated in the cone section (the separation area).

(2) The turbulent flow near the wall of the LLCR was large.The fluid movement near the sidewall was complicated and the vortex motion was very intense.Therefore,the turbulent flow of the LLCR was related to the motion state and droplet breakage.

(3) The deformation and breakage location of droplets were mainly occurred in the velocity boundary layer,while it was difficult to break droplet in the mainstream region.

(4) The broken angle in the boundary layer decreased with the increase of the droplet diameter,the relative slip velocity,the Weber number,and the decrease of the surface tension coefficient,indicating the large surface tension prevented droplet breakage,but large Weber number had a positive effect.

Fig.8.(a) Concentration distribution of light phase and pressure distribution in the LLCR,(b) static pressure,(c) dynamical pressure,(d) total pressure.

Fig.9.(a) Calculation model,(b) droplet injection coordinates.

Fig.10.Morphological changes of droplet C at different time,(a)1×10-4 s,(b)5×10-4 s,(c)10×10-4 s,(d)15×10-4 s,(e)20×10-4 s,(f)25×10-4 s,(g)30×10-4 s,(h)35×10-4 s,(i) 40×10-4 s,(j) 45×10-4 s,(k) 50×10-4 s,(l) 60×10-4 s,(m) 70×10-4 s,(n) 80×10-4 s,(o) 90×10-4 s,(p) final morphology.

Fig.11.Relationship of droplet diameter and broken angle,(a) droplet C,(b) droplet B.

Fig.12.Relationship of slip velocity and broken angle.

Fig.13.Relationship of surface tension coefficient and broken angle.

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 gratefully acknowledge the financial support of the Natural Science Foundation of China (21890763),and National Science and Technology Innovation Leading Talents Project of National Ten Thousand People Plan.

Fig.14.Relationship of Weber number and broken angle.

Nomenclature

CDdrag coefficient

Ddiameter of LLCR,mm

Dodiameter of overflow,mm

Dudiameter of downflow,mm

Foverflow split ratio,%

f drag force,N

L1length of guide vane,mm

L2length of light phase inlet,mm

L3length of column,mm

L4length of cone,mm

L5length of downflow,mm

qilflow rate of light phase at the inlet,mm

qolflow rate of light phase at the overflow,mm

ReReynolds number

rradial radius of LLCR,mm

Δppressure drop,kPa

vflow rate of mixing phases,m·s-1

η separation efficiency,%

ξ resistance coefficient

ρ density of mixing phases,kg·m-3

ω broken angel,(°)

Subscripts

D drag

il light phase at the inlet

o overflow

ol light phase at the overflow

u downflow

1 guide vane

2 light phase inlet

3 column

4 cone

5 downflow