Compressibility induced bubble size variation in bubble column reactors:Simulati

Dongyue Li,Zhipeng Li,Zhengming Gao*

State Key Laboratory of Chemical Resource Engineering,School of Chemical Engineering,Beijing University of Chemical Technology,Beijing 100029,China

Keywords:Bubble column reactors Computational fluid dynamics Population balance equation Compressibility

ABSTRACT Bubble column reactors can be simulated by the two fluid model(TFM)coupled with the population balance equation(PBE).For the large industrialbubble columns,the compressibility due to the pressure difference may introduce notable bubble size variation.In order to address the compressibility effect,the PBE should be reformulated and coupled with the compressible TFM.In this work,the PBE with a compressibility term was formulated from single bubble dynamics,the mean Sauter diameters predicted by the compressible TFM coupled with the PBE were compared with the analyticalsolutions obtained by the idealgas law.Itwas proven thatthe mesoscale formulations presented in this work were physically consistentwith the macroscale modeling.Itcan be used to simulate large industrial plants when the compressibility induced bubble size variation is important.

1.Introduction

Numerous chemical engineering systems involve discrete population of bubbles.One notable example is the bubble column reactors[1].The information of the flow fields in the bubble column can be predicted by the computational fluid dynamics(CFD)and it is common to employ the two fluid model(TFM)[2]to simulate multiphase flows with high phase fractions[3-11].However,the hydrodynamic description and the volume averaging approaches employed in the TFM lead to the problem that it is only applicable for monodisperse systems.For the realistic polydisperse system,the evolution of the bubble population in the bubble column can be modeled by means of the population balance equation(PBE)[12],whose unknown is the bubble size distribution(BSD)associated with the bubble population.The BSD can be further distinguished by the length-based BSD or volume-based BSD[13].By coupling the TFM and the PBE,the knowledge of the BSD in bubble columns can be obtained.

In the past,the TFM-PBE coupling was extensively employed to simulate the bubble columns[9,10,14,15].In the work of Liang et al.[15],the class method was applied to solve the PBE.The TFM-PBE coupling was achieved by the momentum exchange closure term.Three different drag models were employed to investigate the multiphase flow behaviors.The wake acceleration effect of large bubbles was firstly evaluated.In the work of Buffo et al.[14],the quadrature-based moment method was employed.The effect of the bubble swarm was investigated by a correction on the drag model.However,the compressibility induced bubble size variation has never been investigated.For the industrial plants where the operating pressure at the bottom may be much larger than that on the top,the high pressure difference may lead to significant variation of the bubble density and the bubble size.

In such cases the compressibility should be addressed in the TFM and PBE.In this work,the compressible TFM was coupled with the PBE to predict the BSD in conceptual one-dimensional bubble columns with different height.The compressibility induced bubble size variation was novelly introduced by a “growth”term in the PBE which was solved by the method of moments(MOM).The continuous BSD was reconstructed by the log-normal extended quadrature method of moments(EQMOM).Compared with the class method(CM),readers will find that the compressibility induced bubble size variation can be much easily captured in the MOM.Finally,the prediction of the mean Sauter diameter of the bubbles was verified with the analytical solutions.

2.Background Theory

2.1.Compressible TFM

The phase fraction equation of the compressible TFM is reported as the following:

where αdis the volume fraction of the dispersed phase,ρdis the density of the dispersed phase,and Udis the average velocity of the dispersed phase.The volume fraction ofthe continuous phase(αc)is generally calculated by knowing that volume fractions sum to unity(αd+ αc=1).The average velocities ofthe dispersed phase(Ud)and ofthe continuous phase(Uc)are calculated by solving the corresponding momentum balance equations:

where p is the pressure shared by the two phases,τdand τcare the viscous stress tensors,Rdand Rcare the Reynolds-stress tensors,and g is the gravitational acceleration vector.Mdis the interfacial force term,which describes the momentum exchange between the two phases and requires a certain degree of modeling.It is common to break it down into the axial drag force Mdrag;the lateral lift force Mlift,which acts perpendicular to the direction of the relative motion of the two phases;the virtual mass force Mvm,which is proportional to relative phase accelerations;the wall lubrication force Mwallwhich acts to drive the bubble away from the wall;and the turbulent dispersion force Mturb,which is the result of the turbulent fluctuations of the liquid velocity.Specifically,the drag force can be calculated as follows:

where CDis the drag force coefficient,which was calculated by the correlation of Schiller and Naumann[16];ddis the diameter of the dispersed phase,or the mean Sauter diameter in each computational cell if the PBE was coupled with the TFM.The liftforce can be calculated as follows:

where CLis the lift force coefficient.The virtual mass force can be calculated as follows:

where Cvmis the virtual mass force coefficient.Readers interested in the theory of the TFM are suggested to other works[17]and the numerical solving procedure was summarized in our previous work[18].

2.2.Compressible PBE

For bubble column simulations,the internal coordinate of the BSD is often characterized by the bubble diameter in order to investigate the size variation intuitively.The one-dimensional length-based BSD transport equation for n(L,x,t)can be written as the following(the dependence on L,x and t is omitted for simplicity):

where udis the advection velocity for the bubbles.Atthis stage,the possible source terms such as the growth,breakage and coalescence are neglected.To include the compressibility on Eq.(7),one may include the density on each term to address the density variation,and Eq.(7)can be reformulated to where ρdis the bubble density which is not spatially homogeneous for the compressible fluids.Unfortunately,it should be stressed here Eq.(8)does not obey the statistic theory from which the PBE is developed and it leads to wrong results.

The bubble size variation due to the compressibility is reported in Fig.1.

Fig.1.Left:single bubble growth due to the compressibility from the bottomto the top in a bubble column.Right:the “drift”shape of the BSD due to the compressibility.

It can be seen that the bubble size variation is continuous,which implies that it is not physically reasonable to address the compressibility effect as first-order or higher-order point processes[13].It is therefore reasonable to qualify the bubble size variation due to the compressibility by adding a zero-order point process(e.g.,a “growth”term)in Eq.(7).Under such assumption,the BSD transport equation with the growth term can be written as

where G is the diameter-based growth coefficient,which can be calculated by

It represents the continuous rate of change of bubble length L.The unknown growth coefficient introduced by the compressibility should be derived from the single bubble dynamics.The relationship between the pressure and the volume can be represented by the ideal gas law as follows:

where V is the volume,p is the pressure,N is the mole number,R is the universal gas constant,and T is the temperature.Using this assumption,we obtain

Combining Eqs.(12)and(13)yields

Substituting Eq.(14)into Eq.(9)leads to the BSD transportequation with effect of the compressibility due to the pressure difference:

In this work,Eq.(15)was numerically solved by the method of moments(MOM)and the continuous BSD was reconstructed by the extended quadrature method of moments(EQMOM)[19].The corresponding moment transport equations are defined by

The details of the MOM and EQMOM are neglected here.Reader interested in this subject is suggested to other latest works[13].It should be stressed here that Eq.(8)can be also transformed to the moment transport equations as follows:

Itcan be seen thatthe R.H.S.of Eqs.(16)and(17)is different when k≠3.For example,the R.H.S.of Eqs.(16)and(17)for m2isand,respectively.If Eq.(17)was employed to predict themean Sauter diameter(m3/m2),it overpredicts the value compared with the analyticalsolution due to the largersink source term.In the following testcase,the correctness ofEq.(16)was proved by comparisons between numerical predictions with analytical solutions.

Table 1 Numerical schemes and boundary conditions used in the test cases

3.Numerical Test Case

3.1.1-meter-high bubble column

A one-dimensional conceptual bubble column was employed in this work to investigate the effect of the compressibility on the bubble size variation.The grids include 75 cells in height(1 m),as this results in grid-independent predictions.Further re finement of the mesh(e.g.,150 cells)did not improve the predicted results.The computational domain was initially fed by mixture ofwaterand gas.The inletboundary of the gas phase fraction was 9.549%.The inlet gas velocity was 0.1 m·s-1as the Dirichlet boundary condition,and the flow fields were assumed to be laminar.The operating pressure at the top of the column was set to 101325 Pa.The boundary conditions are reported in Table 1.The information of the flow fields was calculated by the compressible TFM solver in OpenFOAM-5.x.

It can be seen in Fig.2 that the plots of the predicted phase fraction by Eq.(1)and thecalculated from Eq.(16)overlap.It is because that when k=3,Eq.(16)can be written as

which is mathematically equivalent to the compressible phase fraction equation[18]:

The overlapped plots verify the correctness of Eq.(16),whereas Eq.(17)is incorrect.Moreover,itcan be seen thatthe gas phase fraction inside the bubble column(around 2.5%)ismuch smallerthan thatatthe inlet(9.549%)due to the upward buoyant force.

Letus now analyze the bubble size variation due to the compressibility.For the testcase investigated in this work,the pressure can be calculated analytically by p=p0+ρcgh.From Eq.(11),the analytical solution for the bubble size can be written as

where L0and p0are the reference diameter and pressure,respectively.The reference pressure can be calculated from the hydrodynamics,and it equals 111115 Pa at the bottom of the one-meter-high column.The corresponding mean Sauter diameter is 3.136 mm.The comparison of the predicted mean Sauter diameter(d32=m3/m2)with the analytical solution calculated by Eq.(20)was reported in Fig.2.It can be seen that the predicted mean Sauter diameter predicted by Eq.(16)agrees well with the analytical solution.

Fig.2.Left:plots ofthe predicted phase fractionαd by Eq.(1)(black line)and by Eq.(16)(red line).Right:plots ofthe mean Sauterdiameter predicted by Eq.(16)(black line)and the analytical solution(red line).

The reconstructed BSD predicted by the log-normalEQMOMwas reported in Fig.3 in order to verify the drift phenomenon as reported in Fig.1.It can be seen that the shape of the BSD shifts horizontally from left to the right,which implies the bubble literally grows due to the compressibility.Such drift phenomenon of the reconstructed BSD predicted by the EQMOM is consistent with Fig.1.

Fig.3.Plots of the reconstructed BSD predicted by the log-normal EQMOM for x=0 m(black line),x=0.5 m(green line)and x=1 m(red line).

Fig.4.Plots of the predicted phase fractionby PBE with(red line)and without(black line)the compressibility induced term.

3.2.9-meter-high bubble column

In this test case,the bubble column was heightened to 9 m,which is a common operating height for industrial gas-liquid flows[20].The plots of the predicted phase fraction are reported in Fig.4.It can be seen that the predicted phase fraction(π6m3)by the PBE with the compressibility induced term captures the bubble volume growth due to the pressure difference in the bubble column.However,the predicted bubble volume by PBE without the compressibility induced term keeps constant due to the absence of the source term in Eq.(18),which is not physically meaningful.The plots of the predicted mean Sauter diameters are reported in Fig.5.It can be seen that the predicted mean Sauter diameter predicted by Eq.(16)agrees well with the analytical solution.It was further proven by that the relative value of the predicted mean Sauter diameter with analytical solution is approximately 1.

Atlast,it is interesting to discuss the relative importance ofthe compressibility induced bubble size variation compared with the breakage/coalescence induced bubble size variation.From the 9-meter-high bubble column test case,the bubble diameter increases from 3.14 mm to 3.83 mm due to the compressibility effect.The variation is approximately 20%which should not be overlooked.Moreover,when the PBE was coupled with the TFM,the predictedand αdshould be identical.If the compressibility term was neglected or not correct(e.g.,Eq.(17)),the predictedby the PBE deviates with the predicted αdby the compressible TFM.Last but not least,the compressibility term should be included together with the breakage/coalescence source term in the PBE,since it scales up/down the bubble size on the basis of the bubble size variation due to the breakage/coalescence.

4.Conclusions

In this work,the bubble size variation due to the compressibility was investigated by the compressible TFM coupled with the PBE.The predicted mean Sauter diameters were compared with the analytical solutions by virtue to the ideal gas law.The drift phenomenon of the BSD was successfully predicted by the log-normal EQMOM.It was proven that the PBE with the compressibility effect derived from single bubble dynamics is physically reasonable compared with Eq.(8).For the 9-meter high bubble column test case,the compressibility induced bubble size variation is approximately 20%which cannotbe overlooked.We suggest employing Eq.(15)in the CFD-PBE simulations for industrial bubble column reactors.

Acknowledgment

Dongyue Liwould like to thank Ronald OertelofHelmholtz-Zentrum Dresden-Rossendorf and Antonio buffo of Politecnico di Torino for the discussion of the bubble dynamics.

Fig.5.Left:plots ofthe mean Sauterdiameterpredicted by Eq.(16)(black line)and the analyticalsolution(red line).Right:plotofthe relative value ofthe predicted mean Sauterdiameter with analytical solution.