- Email: [email protected]
Direct Numerical Simulations of gas–liquid–solid three phase flows
Chemical Engineering Science 100 (2013) 293–299
Contents lists available at SciVerse ScienceDirect
Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces
Direct Numerical Simulations of gas–liquid–solid three phase flows M.W. Baltussen, L.J.H. Seelen, J.A.M. Kuipers, N.G. Deen n Multiphase Reactors Group, Department of Chemical Engineering and Chemistry, Eindhoven University of Technology, P.O. Box 513, 5612 AZ Eindhoven, The Netherlands
H I G H L I G H T S c
c
c
c
G R A P H I C A L
A B S T R A C T
A novel combined Front Tracking Immersed Boundary method is presented. The combined method is verified and validated. The combined method is used to study dense three phase flows. The drag coefficient of the bubbles increased with increasing solids fraction.
a r t i c l e i n f o
abstract
Article history: Received 3 September 2012 Received in revised form 12 February 2013 Accepted 22 February 2013 Available online 1 March 2013
Gas–liquid–solid three phase flows are encountered in for example the Fischer–Tropsch process for the production of synthetic fuels in bubble slurry columns. To predict the hydrodynamics in large slurry bubble columns a multi-scale modeling approach can be used, which accounts for the large variation in time and length scales. In this paper, the smallest scale model has been developed using the Front Tracking model of Dijkhuizen et al. (2010b) and the Immersed Boundary model of Kriebitzsch (2011). In the Front Tracking model, each bubble is tracked separately. Furthermore, the Immersed Boundary method introduces the particle–fluid and the particle–particle (via a hard sphere model) interactions in the model. The resulting hybrid Front Tracking Immersed Boundary model is able to simulate dense three-phase flows and accounts for swarm effects in a fundamental manner. From our simulations we found that the relative drag coefficient of bubbles in three-phase flows seems to increase with increasing solids volume fraction. However, longer averaging periods are needed to derive a fully predictive correlation for the relative drag coefficient with respect to the solid volume fraction. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Multiphase flow Front Tracking model Immersed Boundary method Bubble columns Slurries Hydrodynamics
1. Introduction Fischer–Tropsch synthesis, methanol synthesis and other gasto-liquid processes become more feasible, due to the increase in the international oil prices. These processes are mostly performed in slurry bubble columns, which are considerably larger than the currently operated slurry bubble columns. Understanding of the hydrodynamics and the complex phase interactions of threephase systems is important to arrive at improved scale-up procedures of slurry bubble columns (Kantarci et al., 2005; Wang et al., 2007; Yang et al., 2007). In recent years Multiphase
Computational Fluid Dynamics (MCFD) has become an important tool to simulate complex multiphase flows and enhance our understanding of the complex phase interactions. Industrial slurry bubble columns typically are tens of meters tall. At this scale, it is not possible to resolve the flow phenomena at the level of individual dispersed elements, consequently a multiscale modeling approach is adopted. In such approaches, the larger scale models use closure equations obtained from the smaller scale models (Deen et al., 2010; Roghair et al., 2011a; Yang et al., 2007).
1.1. Multi-scale modeling n
Corresponding author. Tel.: þ31 40 247 3681; fax: þ31 40 247 5833. E-mail address: [email protected] (N.G. Deen).
0009-2509/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2013.02.052
In the multi-scale modeling approach, there are three basic types of models: Euler–Euler models, Euler–Lagrangian models
294
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
and Direct Numerical Simulations (DNS). These three types are shown in Fig. 1 for gas–liquid flow. The Euler–Euler models (Fig. 1C) are able to simulate industrial scale bubble columns. In these models, both the liquid phase and the dispersed phases (the bubbles and the particles) are treated as inter-penetrating continua. However, these models need closures for the bubble– bubble, the particle–particle, the bubble–particle, the bubble– liquid and the particle–liquid interactions (Roghair et al., 2011a). The interaction between the dispersed elements can partly be obtained from the Euler–Lagrangian models, see Fig. 1B. In these models, the dispersed phases are represented in a Lagrangian fashion, while the continuous phase is modeled using the continuum approach. However, these models still need closures for the bubble–liquid and the particle–liquid interactions (Deen et al., 2010; Roghair et al., 2011a; Yang et al., 2007). The bubble–liquid and the particle–liquid interactions can be determined using DNS, Fig. 1A. In these models, the Navier– Stokes equations are solved without assumptions. However, due to computational limitations only a small part of the column can actually be simulated (Roghair et al., 2011a; Yang et al., 2007).
allowed. However, it can be introduced via a sub-grid model. The main objective of this work is to combine the FT method with the Immersed Boundary (IB) method of Kriebitzsch (2011). The combined FT-IB model will be verified and validated. In Section 2, the numerical method will be explained. Section 3 gives the verification and validation results. Furthermore, this section also presents the results of the simulations for three-phase systems. Finally, in Section 4 the conclusions of the work will be presented.
1.2. Objectives
r
In this work, a novel DNS model will be presented for gas–liquid– solid three phase flows. The model Li et al. (2001) used was a combined Euler–Lagrangian and DNS approach. However, this model still needs assumptions for the particle–liquid interactions. Ge and Fan (2006) used a combined Level-Set Immersed Boundary method (LS-IBM) and Jain et al. (2012) used a combined Volume-Of-Fluid Immersed Boundary method (VOF-IBM). The disadvantage of the use of front capturing techniques, like Level-Set and Volume-Of-Fluid, is the unphysical merging of the bubbles. This is a disadvantage when swarm effects are studied (van Sint Annaland et al., 2005). Therefore, the Front Tracking (FT) of Dijkhuizen et al. (2010b) will be used in this work. This method tracks the interface of the bubble explicitly, using marker points. By default coalescence is not
2. Numerical model In this work, a short introduction to the combined FT-IB model is given. A more detailed explanation on the FT model and the IB model is given by Dijkhuizen et al. (2010b) and Kriebitzsch (2011), respectively. 2.1. Governing equations and discretization The combined FT-IB model solves the incompressible Navier– Stokes equation and the continuity equation, as follows: @u ¼ rprr ðuuÞ þ r s þ rg þ F s þ F IB @t
ru ¼0
ð1Þ ð2Þ
Compared to the Navier–Stokes equation two additional source terms appear: the force density due to the surface tension, F s , and the force density due to the presence of the immersed boundaries, F IB . F s is calculated directly from the triangular markers that represent the gas–liquid interface, as shown in Fig. 2. Each marker exerts a tensile force on the neighboring markers. Summing the contributions of all the neighboring markers, the surface force on the central marker is calculated. The resulting force is mapped to the surrounding cells using a mass-weighing function. F IB is introduced to enforce the no-slip boundary condition at the surface of the particles. To this end, a number of force points
Fig. 1. The multi-scale modeling approach for gas–liquid flow (Roghair et al., 2011a). (A) shows a simulation using Direct Numerical Simulations, (B) a simulation with a Euler–Lagrangian model and (C) a simulation with the Euler–Euler model.
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
295
particular particle n consists of a translational part, wn,p , and a rotational part, on,p , as shown in Eq. (3). In this equation, the vector between a force point at the surface and the center of mass of the particle is given by ðr m r i,j,k Þ v p,n ¼ w p,n þ o p,n ðr m r i,j,k Þ
Fig. 2. The Lagrangian presentation of a single bubble. The figure on the right shows an enlargement of the surface mesh.
ð3Þ
Eq. (4) gives the translational part of the velocity of the particle. In this equation, F l-p,n is the hydrodynamic force and F k-n are the forces between particles k and n ðrp,n rl ÞV p,n
X @w n ¼ ðrp,n rl ÞV p,n g þ F l-p,n þ F k-n @t kan
ð4Þ
The rotational part of the velocity can be determined using Eq. (5), in which T k-n is the torque between particles k and n and Ip,n the moment of inertia, which is given by Eq. (6) for spherical particles ðrp,n rl ÞV p,n
Ip,n ¼
X Ip,n @o n ¼ ðr m r i,j,k Þ F l-p,n þ T k-n rp,n @t kan
1 2 mp,n dp,n 10
ð5Þ
ð6Þ
The equations for the advection of the particles are solved using an explicit Euler scheme. After the new velocities of all particles have been obtained, they are moved to their new positions. Inter-particle collisions, which take place during the movement, are handled with a hard sphere model (Hoomans et al., 1996). Fig. 3. The Lagrangian presentation of a particle. The blue sphere shows the particle, and the red points are the force points used for enforcing the no-slip boundary condition. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
are distributed over the surface of the particle as shown in Fig. 3. The forcing term is calculated using an iterative procedure. The velocity at the discrete force points is first estimated, using Eq. (1) without F IB . Using the difference between the velocity at the force point and the particle velocity, F IB is estimated. Using this estimated immersed boundary force, the velocity at the force point is again calculated using Eq. (1). The immersed boundary force can then be adjusted using this new velocity at the force point. This procedure will be repeated until the velocity at the force point is equal to the particle velocity. The Navier–Stokes equations are solved on a staggered grid. To solve Eqs. (1) and (2), a projection-correction method is used. First, the velocity is determined using Eq. (1). The convective term is discretized using a second order flux-delimited Barton scheme. The diffusion term is discretized using a second order central scheme. All the contributions in Eq. (1) are treated explicitly, except for the diffusion term, which is treated semi-explicit. The implicit part of the diffusion term is chosen such that the velocities in the different directions can be solved separately. Subsequently, the solution for the velocity is corrected using Eq. (2). Both these steps are calculated using a ICCG matrix solver. When the velocity field is known, the position of the particles and the bubbles can be adjusted, which will be explained in detail in respectively Sections 2.2 and 2.3. Using the markers of the bubbles, the density can be determined via volume-weighing averaging. While the viscosity can be determined from the phase fraction using harmonic averaging. 2.2. Particles Knowing the fluid phase velocity field, the particle positions need to be updated. First of all, the velocity of the particles should be determined using the second law of Newton. The velocity of a
2.3. Bubbles After updating the position of the particles, the position of the bubbles can be updated. Due to the one-field approximation, the velocity of the marker points can be determined from a piecewise cubic spline interpolation, using the Eulerian grid cells around the marker points. The markers are displaced using a fourth order Runge–Kutta time stepping. By advecting each marker point separately, the bubble changes its position and its shape. However, due to this, the distance between the marker points changes, which leads to a lower quality of the surface mesh. To overcome these problems, a remeshing procedure is used. This procedure consists of four elementary remeshing operations: edge splitting, edge collapsing, edge swapping and smoothing (Roghair, 2012). The Front Tracking method is not volume conservative. Due to the separate advection of each marker point and the remeshing procedures, the bubble volume will change during the simulation. Due to the large number of time steps used in the simulation, the small volume changes per time step will become significant in a simulation. To conserve the volume in the bubbles, the method of Kuprat et al. (2001) is implemented in the smoothing procedure.
3. Results 3.1. Verification To verify the combined model, several tests have been performed. First of all, the implementation of the immersed boundary force was verified using Stokes flow around a sphere. In this test, a single sphere is subjected to fluid flow in the Stokes regime where the drag force acting on the sphere is calculated. The simulation settings are given in Table 1. The simulations were performed using cubic domains with several solid volume fractions (by changing the width of the domain). The results of these simulations are similar to the analytical results derived by
296
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
Table 1 The simulation settings for a single sphere subjected to fluid flow in the Stokes regime. The computational grid is a cubic domain. Solid fractions Computational grid Grid size (m)
0.01y0.002 74y128
Time step (s)
6:6 104
Particle diameter (m)
1:0 103 1.0
Gas density (kg m3 ) Gas viscosity (Pa s) Inlet velocity (m s1 )
5:0 105
1:0 105 1:0 103
Table 2 The simulation settings for a settling sphere. Computational grid (W D H) Particle positions Grid size (m)
160 160 180 (80,80,40)
Time step (s)
1:0 105
Particle diameter (m)
1:0 103 1.0
Fluid density (kg m3 ) Fluid viscosity (Pa s) Inlet velocity (m s1 )
5:0 105
1:0 105 1.0
Hasimoto (1959), for a 2D domain. The deviation is maximally 2% at low solid fractions. At higher solid fractions ( Z 0:8%) the difference between the results increases to 4%, which is caused by the fact that the comparison with the drag obtained by Hasimoto (1959) only holds for low solid fractions. Subsequently, the implementation of the surface tension model was tested using a stationary bubble test of van Sint Annaland et al. (2005). In this test a bubble is placed in a zerogravity field. The excess pressure in the bubble, the Laplace pressure, will in this case exactly counter balance the surface forces. The simulation settings used are similar to the settings of van Sint Annaland et al. (2005). The resulting pressure jump is equal to 200.02 Pa, which only deviates 0.01% from the Laplace pressure. To test the advection of the surface of the bubbles, a standard advection test is performed. In this test, a bubble, with a diameter of 1 1 3 1 4 the domain, is initialized at (2, 4, 2). In the domain, there is a constant rotational velocity. Halfway the simulation, the direction of the velocity field is reversed. It is expected that the bubble returns to its original position. Again the same simulation settings of van Sint Annaland et al. (2005) are used in this simulation. However, the total simulation time is decreased to 1 s. The computational error of the position of the bubble is 0.5%, which is similar to the results of Rider and Kothe (1998) and Dijkhuizen et al. (2010b) and better than the results of van Sint Annaland et al. (2005). Furthermore, the volume conservative remeshing procedure decreased the error in the volume of the bubble to 5 1011 %. Finally, the implementation of the flow solver was tested using an oscillating bubble test. In this test, a bubble is placed in a zerogravity field. The bubble is a sphere which has a slightly larger diameter in one direction. During the test, the bubble diameter will start to oscillate. Lamb (1932) derived an analytical solution for the frequency of these oscillations. The simulation settings are equal to the settings of Dijkhuizen et al. (2010b). The difference between the simulation results and the analytical solution is 0.7% for the oscillation period and 2.5% for the amplitude of the oscillation, which is similar to the results reported by Dijkhuizen (2008).
On the basis of these tests, we conclude that we have correctly implemented the combined FT-IB model. 3.2. Validation Besides verification of the model, the combined FT-IB model should also be validated. First of all, a single sphere is subjected to a flow with a Reynolds number of 100. This will ensure that the model is also able to predict flows with high Reynolds numbers (compared to Stokes flow). The simulation settings used in this paper are shown in Table 2. The resulting drag coefficient is 1.269 (using the forcing diameter), which is equal to the results of Beetstra et al. (2007) and Deen et al. (2009). In a bubble slurry column, the solid fraction is much higher than the performed simulations. Therefore, a test is performed with a cubic array of particles. In the test the simulation settings are equal to the settings of Deen et al. (2009) with a domain of 1603. The dimensionless drag force is determined for the central particle in the array. The resulting dimensionless drag force is 19.2. To compare the results with the results of Zick and Homsy (1982), Deen et al. (2009) made a fit of the data of Zick and Homsy (1982). The obtained drag force is within the 3% accuracy of the fit. Besides the comparison with other codes, the FT-IBM code, can also be compared to experiments. In these simulations, the terminal settling velocity of three particles is determined. The parameters used in these simulations are given in Table 3. The table also shows the experimental and the simulation results. The simulation results give similar results as the experiments. The table shows that the terminal velocity determined in the simulations is similar to the experimental results (Mordant and Pinton, 2000) and the terminal velocity calculated with the correlation of Clift et al. (1978). In the simulations with high Reynolds numbers, the results are influenced by the confinement of the wake in the simulation domain. Due to this confinement, the wake could not fully develop. This might explain the difference between the simulation results and the experimental and correlation results. To validate the FT part of the combined code, five single bubble simulations have been performed in different regimes. The used ¨ os ¨ and Morton numbers are given in Table 4. The resulting Eotv Reynolds numbers are also given in the table. The results of the simulations are similar to the results of Dijkhuizen et al. (2010a). The small differences are caused by slightly different positioning of the bubble in the domain (Dijkhuizen et al., 2010a places his bubbles at (0.5 0.5 0.6) domain size while in the simulations presented here the bubble is placed at (0.5 0.5 0.67) domain size). Furthermore, the code used by Dijkhuizen et al. (2010a) did not include volume conservative remeshing. Therefore, the bubbles decreased in size during the simulation. The difference between Tomiyama (1998) and the presented results can be caused by the slightly contaminated liquids used by Tomiyama (Dijkhuizen et al., 2010a). From these validation results, we conclude that the combined FT-IB model is able to simulate particles and bubbles correctly. 3.3. Three phase systems Following the verification and the validation of the FT-IB method, the method is used for the simulation of three phase systems using periodic domains. In total, 10 test cases were simulated using the combined model. The settings used in the simulations are shown in Table 5. The number of bubbles is fixed to 16. Roghair et al. (2011a) established that 16 bubbles are enough to simulate a bubble swarm using periodic boundary conditions. The size of the cubic domain was adjusted based on the gas volume fraction. Subsequently, the number of particles used in the simulations was determined from the size of the domain and the imposed particle volume fraction.
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
297
Table 3 The parameters used for settling sphere cases. The experimental results for the terminal rise velocity, vt, are taken from Mordant and Pinton (2000). Furthermore, the terminal velocity of Clift et al. (1978) is calculated with the correlation proposed by Clift et al. (1978). Case
1
2
3
Computational grid (W D H) Particle positions Grid size (m)
140 140 180 (70,70,40)
120 120 220 (60,60,40)
120 120 260 (60,60,40)
2:5 105
4:0 105
5:0 105
Time step (s)
1:0 105
1:0 105
1:0 105
Particle diameter (m)
5:0 104 (0,0, 9.81) 1000
8:0 104 (0,0, 9.81) 1000
1:0 103 (0,0, 9.81) 1000
Particle density (kg m3 )
9:0 104 2560
9:0 104 7710
9:1 104 7850
Experimental vt (m s1 ) vt of Clift et al. (1978) (m s1 ) Computed vt (m s1 )
0.074 0.076 0.076
0.316 0.325 0.312
0.383 0.393 0.372
g (m s2 ) Liquid density (kg m3 ) Liquid viscosity (Pa s)
Table 4 The parameters used in the simulation of single rising bubbles. The resulting Reynolds numbers are also shown in this table. For comparison, the experimental results of Tomiyama (1998) and the simulation results of Dijkhuizen et al. (2010a) are also shown. Case
1
2
3
4
5
Mo
1:0 103 1.0
1:0 103 10
1.0
Eo
100
1:0 103 100
1:0 1012 2.0
Re of Tomiyama (1998)
2.11
26.5
22.8
2.10
2:06 103
Re of Dijkhuizen et al. (2010a)
2.21
23.6
18.4
2.07
2:33 103
Computed Re
2.33
24.2
19.7
1.99
2:36 103
Table 5 The simulation settings for the three phase systems. The grid size of the cubic grid and the number of particle used in the simulations are determined by the volume fraction of gas and the volume fraction of particles used. Bubble diameter (m)
2:0 103
Particle diameter (m)
1:0 103
Grid size (m)
5:0 105
Time step (s)
1:0 105
Fluid density (kg m3 )
1:0 103
Gas density (kg m3 )
1:0 102
Particle density (kg m3 ) Gas viscosity (Pa s)
2:0 103
Fluid viscosity (Pa s) Number of bubbles Number of particles Bubble volume fraction (–) Particle volume fraction (–)
1:8 105 1:0 103 16 8y64 0.20, 0.30, 0.40 0.025, 0.05, 0.075, 0.10
force and F W the forces due to the walls mb
dvb ¼ F G þ F P þ F D þ F L þ F VM þ F W dt
Assuming a pseudo-steady state, the time-averaged force balance reduces to Eq. (8). In this equation, the virtual mass force is zero due to the steady state assumption and the lift force is zero when the velocity gradients are calculated over a length scale equal to the computational domain size. Furthermore, the walls are not taken into account in the DNS simulations. Therefore, the F W will be equal to zero /F D S ¼ /F G S þ/F P S
ð8Þ
1 2 C D rl pdb 9vb u9ðvuÞ ¼ V b rg g þ rpV b 8
ð9Þ
1 2 C D rl pdb 9vb u9ðvuÞ ¼ V b g ðrg þ ðarg þ frs þ ð1afÞrl ÞÞ 8 Fig. 4 shows two typical snapshots of the three phase simulations. Both configurations clearly show a horizontal alignment of the bubbles, see Fig. 4c and f. A similar type of clustering was also observed for bubble swarms by Roghair et al. (2011b). The horizontal alignment of the bubbles is probably caused by the lack of large scale circulation phenomena and the absence of walls in the simulations. Furthermore, due to the alignment of the bubbles, the particles also show preference for horizontal and vertical alignment (Fig. 4b and e). The results of the DNS simulations can provide a closure for the drag force in the Discrete Bubble Model (DBM) (i.e. the unresolved model using Lagrangian tracking of the dispersed elements). In the DBM, the velocity of the bubbles is calculated using Newton’s second law, Eq. (7). In this equation v b is the bubble velocity, F G the gravity force, F P the force due to pressure gradients, F D the drag force, F L the lift force, F VM the virtual mass
ð7Þ
ð10Þ
Insertion of the expressions for the definition of the drag force and the pressure force in Eq. (8) leads to Eq. (10). In this equation db is the bubble diameter, Vb the volume of the bubble, f the solids volume fraction and a the gas volume fraction. Rearranging Eq. (10), the relative drag coefficient can be determined, according to Eq. (11) (Roghair et al., 2011a). C D,1 and v b,1 in Eq. (11) represent the drag coefficient and the velocity of a single bubble in an infinite, initially quiescent liquid, respectively, as defined in Roghair et al. (2011a) C D,rel ¼
CD
rs rl rl rg
C D,1 ð1aÞ þ f
!¼
/v b,1 S2 ð/vuSÞ2
ð11Þ
The relative drag coefficient is determined for all 10 simulations. The results are presented in Fig. 5. The figure shows an increase in
298
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
Fig. 4. Two snapshots (a) and (d) of a three phase simulation containing 30% bubbles and 7.5% particles. (b) and (e) show the configuration of the particles in the simulation, while (c) and (f) the configuration of the bubbles.
20
2.5
CD,rel three phases/CD, rel two phases
α = 20 % α = 30 % α = 40 %
CD, rel (−)
15
10
5
0
2
4
6
8
10
φ (%)
α = 20 % α = 30 % α = 40 %
2
1.5
1
0
2
4
6
8
10
φ (%)
Fig. 5. The effect of the solids fraction (f) on the relative drag coefficient. The line around the point a ¼ 40% and f ¼ 2:5% shows the standard deviation of the determined relative drag coefficient.
Fig. 6. The effect of the solids fraction (f) on the relative drag coefficient over the relative drag coefficient of a bubble swarm with the same void fraction. The dashed line is the trend line.
the drag coefficient with an increase in the solids volume fraction. This can be explained by the particles that obstruct the bubbles during their rise. Fig. 5 also shows the standard deviation of the simulation with an a of 40% and a f of 2.5%. The value for the standard deviation is typical for all the simulations. The relative large standard deviation can be explained by the short averaging periods used in these simulations (0.1–0.35 s), however, the trend is
very clear. The development of a fully predictive correlation will be the subject of a future contribution. Therefore, the relative drag coefficient is divided by the relative drag coefficient of a bubble swarm (which can be calculated using Eq. (11) with a f of zero). The results are shown in Fig. 6. All the data points seem to coincide with the trend line shown in this figure.
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
299
4. Conclusions
Acknowledgments
In this paper, a new model is presented for the Direct Numerical Simulation of three phase flows. The model combines the FT model of Dijkhuizen et al. (2010b) and the IB model of Kriebitzsch (2011). The combined model is verified and validated for both the FT and the IB part of the model. The FT-IB model was used to simulate dense threephase flows. The relative drag coefficient of the bubbles in the simulation of the three-phase flow seems to increase with an increasing solids fraction. However, longer averaging periods are needed to derive a fully predictive correlation for the relative drag coefficient with respect to the solid volume fraction.
The authors would like to thank the European Research Council for its financial support, under its Starting Investigator Grant scheme, contract number 259521 (CuttingBubbles).
Nomenclature Roman symbols C d F g I m p r T t u v w V
Coefficient Diameter, m Force, N Gravity constant, m s 2 Moment of inertia, kg m2 Mass, kg Pressure, Pa Position, m Torque, m N Time, s Liquid, fluid velocity, m s 1 Particle velocity, m s 1 Translation velocity, m s 1 Volume, m3
Greek symbols
a m f
r s o Abbreviations and subscripts b D DBM DNS FT g G IB i, j, k l L m Mo Eo p P Re s t VM W
s 1
Gas volume fraction Viscosity, Pa s Solid volume fraction Density, kg m 3 Stress tensor Rotational velocity, s 1 Bubble Drag Discrete Bubble Model Direct Numerical Simulations Front Tracking Gas phase Gravity Immersed boundary Center of mass Liquid phase Lift Marker Morton number ¨ os ¨ number Eotv Particle Pressure Reynolds number Solid Terminal Virtual mass Wall Surface tension Single bubble infinite liquid
References Beetstra, R., Van Der Hoef, M.A., Kuipers, J.A.M., 2007. Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres. AIChE J. 53 (2), 489–501. Clift, R., Grace, J., Weber, M., 1978. Bubbles, Drops and Particles. Academic Press, New York. Deen, N.G., Mudde, R., Kuipers, J.A.M., Zehner, P., Kraume, M., 2010. Bubble columns. In: Ullmann’s Encyclopedia of Industrial Chemistry. Wiley-VCH Verlag GmbH & Co, KGaA, Weinheim. Deen, N.G., van Sint Annaland, M., Kuipers, J.A.M., 2009. Direct numerical simulation of complex multi-fluid flows using a combined front tracking and immersed boundary method. Chem. Eng. Sci. 64 (9), 2186–2201. Dijkhuizen, W., 2008. Derivation of Closures for Bubbly Flows Using Direct Numerical Simulations. Ph.D. Thesis, University of Twente. Dijkhuizen, W., Roghair, I., Van Sint Annaland, M., Kuipers, J.A.M., 2010a. DNS of gas bubbles behaviour using an improved 3D front tracking model—drag force on isolated bubbles and comparison with experiments. Chem. Eng. Sci. 65 (4), 1415–1426. Dijkhuizen, W., Roghair, I., Van-Sint-Annaland, M., Kuipers, J.A.M., 2010b. DNS of gas bubbles behaviour using an improved 3D front tracking model—model development. Chem. Eng. Sci. 65 (4), 1427–1437. Ge, Y., Fan, L.-S., 2006. 3-d direct numerical simulation of gas–liquid and gas–liquid– solid flow systems using the level-set and immersed-boundary methods. Adv. Chem. Eng. 31, 1–63. Hasimoto, H., 1959. On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres. J. Fluid Mech. 5 (02), 317–328. Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J., Van Swaaij, W.P.M., 1996. Discrete particle simulation of bubble and slug formation in a two-dimensional gasfluidised bed: a hard-sphere approach. Chem. Eng. Sci. 51 (1), 99–118. Jain, D., Deen, N.G., Kuipers, J.A.M., Antonyuk, S., Heinrich, S., 2012. Direct numerical simulation of particle impact on thin liquid films using a combined volume of fluid and immersed boundary method. Chem. Eng. Sci. 69 (1), 530–540. Kantarci, N., Borak, F., Ulgen, K., 2005. Bubble column reactors. Process Biochem. 40 (7), 2263–2283. Kriebitzsch, S., 2011. Direct Numerical Simulation of Dense Gas–Solid Flows. Ph.D. Thesis, Eindhoven, University of Technology. Kuprat, A., Khamayseh, A., George, D., Larkey, L., 2001. Volume conserving smoothing for piecewise linear curves, surfaces, and triple lines. J. Comput. Phys. 172, 99–118. Lamb, H., 1932. Hydrodynamics, 6th ed. Cambridge University Press, London. Li, Y., Yang, G., Zhang, J., Fan, L., 2001. Numerical studies of bubble formation dynamics in gas–liquid–solid fluidization at high pressures. Powder Technol. 116 (2–3), 246–260. Mordant, N., Pinton, J.-F., 2000. Velocity measurement of a settling sphere. Eur. Phys. J. B 18 (2), 343–352. Rider, W.J., Kothe, D.B., 1998. Reconstructing volume tracking. J. Comput. Phys. 141, 112. Roghair, I., 2012. Direct Numerical Simulations of Hydrodynamics and Mass Transfer in Dense Bubbly Flows. Ph.D. Thesis, Eindhoven, University of Technology. Roghair, I., Lau, Y., Deen, N.G., Slagter, H., Baltussen, M., Annaland, M.V.S., Kuipers, J.A.M., 2011a. On the drag force of bubbles in bubble swarms at intermediate and high Reynolds numbers. Chem. Eng. Sci. 66, 3204–3211. Roghair, I., Martı´nez Mercado, J., Annaland, M.V.S., Kuipers, H., Sun, C., Lohse, D., 2011b. Energy spectra and bubble velocity distributions in pseudo-turbulence: numerical simulations vs. experiments. Int. J. Multiphase Flow 37 (9), 1093–1098. Tomiyama, A., 1998. Struggle with computational bubble dynamics. In: Third International Conference on Multiphase Flow, pp. 369–405. van Sint Annaland, M., Deen, N.G., Kuipers, J.A.M., 2005. Numerical simulation of gas bubbles behaviour using a three-dimensional volume of fluid method. Chem. Eng. Sci. 60 (11), 2999–3011. Wang, T., Wang, J., Jin, Y., 2007. Slurry reactors for gas-to-liquid processes: a review. Ind. Eng. Chem. Res. 46 (18), 5824–5847. Yang, G., Du, B., Fan, L., 2007. Bubble formation and dynamics in gas–liquid–solid fluidization—a review. Chem. Eng. Sci. 62 (1–2), 2–27. Zick, A.A., Homsy, G.M., 1982. Stokes flow through periodic arrays of spheres. J. Fluid Mech. 115, 13–26.
Contents lists available at SciVerse ScienceDirect
Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces
Direct Numerical Simulations of gas–liquid–solid three phase flows M.W. Baltussen, L.J.H. Seelen, J.A.M. Kuipers, N.G. Deen n Multiphase Reactors Group, Department of Chemical Engineering and Chemistry, Eindhoven University of Technology, P.O. Box 513, 5612 AZ Eindhoven, The Netherlands
H I G H L I G H T S c
c
c
c
G R A P H I C A L
A B S T R A C T
A novel combined Front Tracking Immersed Boundary method is presented. The combined method is verified and validated. The combined method is used to study dense three phase flows. The drag coefficient of the bubbles increased with increasing solids fraction.
a r t i c l e i n f o
abstract
Article history: Received 3 September 2012 Received in revised form 12 February 2013 Accepted 22 February 2013 Available online 1 March 2013
Gas–liquid–solid three phase flows are encountered in for example the Fischer–Tropsch process for the production of synthetic fuels in bubble slurry columns. To predict the hydrodynamics in large slurry bubble columns a multi-scale modeling approach can be used, which accounts for the large variation in time and length scales. In this paper, the smallest scale model has been developed using the Front Tracking model of Dijkhuizen et al. (2010b) and the Immersed Boundary model of Kriebitzsch (2011). In the Front Tracking model, each bubble is tracked separately. Furthermore, the Immersed Boundary method introduces the particle–fluid and the particle–particle (via a hard sphere model) interactions in the model. The resulting hybrid Front Tracking Immersed Boundary model is able to simulate dense three-phase flows and accounts for swarm effects in a fundamental manner. From our simulations we found that the relative drag coefficient of bubbles in three-phase flows seems to increase with increasing solids volume fraction. However, longer averaging periods are needed to derive a fully predictive correlation for the relative drag coefficient with respect to the solid volume fraction. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Multiphase flow Front Tracking model Immersed Boundary method Bubble columns Slurries Hydrodynamics
1. Introduction Fischer–Tropsch synthesis, methanol synthesis and other gasto-liquid processes become more feasible, due to the increase in the international oil prices. These processes are mostly performed in slurry bubble columns, which are considerably larger than the currently operated slurry bubble columns. Understanding of the hydrodynamics and the complex phase interactions of threephase systems is important to arrive at improved scale-up procedures of slurry bubble columns (Kantarci et al., 2005; Wang et al., 2007; Yang et al., 2007). In recent years Multiphase
Computational Fluid Dynamics (MCFD) has become an important tool to simulate complex multiphase flows and enhance our understanding of the complex phase interactions. Industrial slurry bubble columns typically are tens of meters tall. At this scale, it is not possible to resolve the flow phenomena at the level of individual dispersed elements, consequently a multiscale modeling approach is adopted. In such approaches, the larger scale models use closure equations obtained from the smaller scale models (Deen et al., 2010; Roghair et al., 2011a; Yang et al., 2007).
1.1. Multi-scale modeling n
Corresponding author. Tel.: þ31 40 247 3681; fax: þ31 40 247 5833. E-mail address: [email protected] (N.G. Deen).
0009-2509/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2013.02.052
In the multi-scale modeling approach, there are three basic types of models: Euler–Euler models, Euler–Lagrangian models
294
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
and Direct Numerical Simulations (DNS). These three types are shown in Fig. 1 for gas–liquid flow. The Euler–Euler models (Fig. 1C) are able to simulate industrial scale bubble columns. In these models, both the liquid phase and the dispersed phases (the bubbles and the particles) are treated as inter-penetrating continua. However, these models need closures for the bubble– bubble, the particle–particle, the bubble–particle, the bubble– liquid and the particle–liquid interactions (Roghair et al., 2011a). The interaction between the dispersed elements can partly be obtained from the Euler–Lagrangian models, see Fig. 1B. In these models, the dispersed phases are represented in a Lagrangian fashion, while the continuous phase is modeled using the continuum approach. However, these models still need closures for the bubble–liquid and the particle–liquid interactions (Deen et al., 2010; Roghair et al., 2011a; Yang et al., 2007). The bubble–liquid and the particle–liquid interactions can be determined using DNS, Fig. 1A. In these models, the Navier– Stokes equations are solved without assumptions. However, due to computational limitations only a small part of the column can actually be simulated (Roghair et al., 2011a; Yang et al., 2007).
allowed. However, it can be introduced via a sub-grid model. The main objective of this work is to combine the FT method with the Immersed Boundary (IB) method of Kriebitzsch (2011). The combined FT-IB model will be verified and validated. In Section 2, the numerical method will be explained. Section 3 gives the verification and validation results. Furthermore, this section also presents the results of the simulations for three-phase systems. Finally, in Section 4 the conclusions of the work will be presented.
1.2. Objectives
r
In this work, a novel DNS model will be presented for gas–liquid– solid three phase flows. The model Li et al. (2001) used was a combined Euler–Lagrangian and DNS approach. However, this model still needs assumptions for the particle–liquid interactions. Ge and Fan (2006) used a combined Level-Set Immersed Boundary method (LS-IBM) and Jain et al. (2012) used a combined Volume-Of-Fluid Immersed Boundary method (VOF-IBM). The disadvantage of the use of front capturing techniques, like Level-Set and Volume-Of-Fluid, is the unphysical merging of the bubbles. This is a disadvantage when swarm effects are studied (van Sint Annaland et al., 2005). Therefore, the Front Tracking (FT) of Dijkhuizen et al. (2010b) will be used in this work. This method tracks the interface of the bubble explicitly, using marker points. By default coalescence is not
2. Numerical model In this work, a short introduction to the combined FT-IB model is given. A more detailed explanation on the FT model and the IB model is given by Dijkhuizen et al. (2010b) and Kriebitzsch (2011), respectively. 2.1. Governing equations and discretization The combined FT-IB model solves the incompressible Navier– Stokes equation and the continuity equation, as follows: @u ¼ rprr ðuuÞ þ r s þ rg þ F s þ F IB @t
ru ¼0
ð1Þ ð2Þ
Compared to the Navier–Stokes equation two additional source terms appear: the force density due to the surface tension, F s , and the force density due to the presence of the immersed boundaries, F IB . F s is calculated directly from the triangular markers that represent the gas–liquid interface, as shown in Fig. 2. Each marker exerts a tensile force on the neighboring markers. Summing the contributions of all the neighboring markers, the surface force on the central marker is calculated. The resulting force is mapped to the surrounding cells using a mass-weighing function. F IB is introduced to enforce the no-slip boundary condition at the surface of the particles. To this end, a number of force points
Fig. 1. The multi-scale modeling approach for gas–liquid flow (Roghair et al., 2011a). (A) shows a simulation using Direct Numerical Simulations, (B) a simulation with a Euler–Lagrangian model and (C) a simulation with the Euler–Euler model.
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
295
particular particle n consists of a translational part, wn,p , and a rotational part, on,p , as shown in Eq. (3). In this equation, the vector between a force point at the surface and the center of mass of the particle is given by ðr m r i,j,k Þ v p,n ¼ w p,n þ o p,n ðr m r i,j,k Þ
Fig. 2. The Lagrangian presentation of a single bubble. The figure on the right shows an enlargement of the surface mesh.
ð3Þ
Eq. (4) gives the translational part of the velocity of the particle. In this equation, F l-p,n is the hydrodynamic force and F k-n are the forces between particles k and n ðrp,n rl ÞV p,n
X @w n ¼ ðrp,n rl ÞV p,n g þ F l-p,n þ F k-n @t kan
ð4Þ
The rotational part of the velocity can be determined using Eq. (5), in which T k-n is the torque between particles k and n and Ip,n the moment of inertia, which is given by Eq. (6) for spherical particles ðrp,n rl ÞV p,n
Ip,n ¼
X Ip,n @o n ¼ ðr m r i,j,k Þ F l-p,n þ T k-n rp,n @t kan
1 2 mp,n dp,n 10
ð5Þ
ð6Þ
The equations for the advection of the particles are solved using an explicit Euler scheme. After the new velocities of all particles have been obtained, they are moved to their new positions. Inter-particle collisions, which take place during the movement, are handled with a hard sphere model (Hoomans et al., 1996). Fig. 3. The Lagrangian presentation of a particle. The blue sphere shows the particle, and the red points are the force points used for enforcing the no-slip boundary condition. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
are distributed over the surface of the particle as shown in Fig. 3. The forcing term is calculated using an iterative procedure. The velocity at the discrete force points is first estimated, using Eq. (1) without F IB . Using the difference between the velocity at the force point and the particle velocity, F IB is estimated. Using this estimated immersed boundary force, the velocity at the force point is again calculated using Eq. (1). The immersed boundary force can then be adjusted using this new velocity at the force point. This procedure will be repeated until the velocity at the force point is equal to the particle velocity. The Navier–Stokes equations are solved on a staggered grid. To solve Eqs. (1) and (2), a projection-correction method is used. First, the velocity is determined using Eq. (1). The convective term is discretized using a second order flux-delimited Barton scheme. The diffusion term is discretized using a second order central scheme. All the contributions in Eq. (1) are treated explicitly, except for the diffusion term, which is treated semi-explicit. The implicit part of the diffusion term is chosen such that the velocities in the different directions can be solved separately. Subsequently, the solution for the velocity is corrected using Eq. (2). Both these steps are calculated using a ICCG matrix solver. When the velocity field is known, the position of the particles and the bubbles can be adjusted, which will be explained in detail in respectively Sections 2.2 and 2.3. Using the markers of the bubbles, the density can be determined via volume-weighing averaging. While the viscosity can be determined from the phase fraction using harmonic averaging. 2.2. Particles Knowing the fluid phase velocity field, the particle positions need to be updated. First of all, the velocity of the particles should be determined using the second law of Newton. The velocity of a
2.3. Bubbles After updating the position of the particles, the position of the bubbles can be updated. Due to the one-field approximation, the velocity of the marker points can be determined from a piecewise cubic spline interpolation, using the Eulerian grid cells around the marker points. The markers are displaced using a fourth order Runge–Kutta time stepping. By advecting each marker point separately, the bubble changes its position and its shape. However, due to this, the distance between the marker points changes, which leads to a lower quality of the surface mesh. To overcome these problems, a remeshing procedure is used. This procedure consists of four elementary remeshing operations: edge splitting, edge collapsing, edge swapping and smoothing (Roghair, 2012). The Front Tracking method is not volume conservative. Due to the separate advection of each marker point and the remeshing procedures, the bubble volume will change during the simulation. Due to the large number of time steps used in the simulation, the small volume changes per time step will become significant in a simulation. To conserve the volume in the bubbles, the method of Kuprat et al. (2001) is implemented in the smoothing procedure.
3. Results 3.1. Verification To verify the combined model, several tests have been performed. First of all, the implementation of the immersed boundary force was verified using Stokes flow around a sphere. In this test, a single sphere is subjected to fluid flow in the Stokes regime where the drag force acting on the sphere is calculated. The simulation settings are given in Table 1. The simulations were performed using cubic domains with several solid volume fractions (by changing the width of the domain). The results of these simulations are similar to the analytical results derived by
296
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
Table 1 The simulation settings for a single sphere subjected to fluid flow in the Stokes regime. The computational grid is a cubic domain. Solid fractions Computational grid Grid size (m)
0.01y0.002 74y128
Time step (s)
6:6 104
Particle diameter (m)
1:0 103 1.0
Gas density (kg m3 ) Gas viscosity (Pa s) Inlet velocity (m s1 )
5:0 105
1:0 105 1:0 103
Table 2 The simulation settings for a settling sphere. Computational grid (W D H) Particle positions Grid size (m)
160 160 180 (80,80,40)
Time step (s)
1:0 105
Particle diameter (m)
1:0 103 1.0
Fluid density (kg m3 ) Fluid viscosity (Pa s) Inlet velocity (m s1 )
5:0 105
1:0 105 1.0
Hasimoto (1959), for a 2D domain. The deviation is maximally 2% at low solid fractions. At higher solid fractions ( Z 0:8%) the difference between the results increases to 4%, which is caused by the fact that the comparison with the drag obtained by Hasimoto (1959) only holds for low solid fractions. Subsequently, the implementation of the surface tension model was tested using a stationary bubble test of van Sint Annaland et al. (2005). In this test a bubble is placed in a zerogravity field. The excess pressure in the bubble, the Laplace pressure, will in this case exactly counter balance the surface forces. The simulation settings used are similar to the settings of van Sint Annaland et al. (2005). The resulting pressure jump is equal to 200.02 Pa, which only deviates 0.01% from the Laplace pressure. To test the advection of the surface of the bubbles, a standard advection test is performed. In this test, a bubble, with a diameter of 1 1 3 1 4 the domain, is initialized at (2, 4, 2). In the domain, there is a constant rotational velocity. Halfway the simulation, the direction of the velocity field is reversed. It is expected that the bubble returns to its original position. Again the same simulation settings of van Sint Annaland et al. (2005) are used in this simulation. However, the total simulation time is decreased to 1 s. The computational error of the position of the bubble is 0.5%, which is similar to the results of Rider and Kothe (1998) and Dijkhuizen et al. (2010b) and better than the results of van Sint Annaland et al. (2005). Furthermore, the volume conservative remeshing procedure decreased the error in the volume of the bubble to 5 1011 %. Finally, the implementation of the flow solver was tested using an oscillating bubble test. In this test, a bubble is placed in a zerogravity field. The bubble is a sphere which has a slightly larger diameter in one direction. During the test, the bubble diameter will start to oscillate. Lamb (1932) derived an analytical solution for the frequency of these oscillations. The simulation settings are equal to the settings of Dijkhuizen et al. (2010b). The difference between the simulation results and the analytical solution is 0.7% for the oscillation period and 2.5% for the amplitude of the oscillation, which is similar to the results reported by Dijkhuizen (2008).
On the basis of these tests, we conclude that we have correctly implemented the combined FT-IB model. 3.2. Validation Besides verification of the model, the combined FT-IB model should also be validated. First of all, a single sphere is subjected to a flow with a Reynolds number of 100. This will ensure that the model is also able to predict flows with high Reynolds numbers (compared to Stokes flow). The simulation settings used in this paper are shown in Table 2. The resulting drag coefficient is 1.269 (using the forcing diameter), which is equal to the results of Beetstra et al. (2007) and Deen et al. (2009). In a bubble slurry column, the solid fraction is much higher than the performed simulations. Therefore, a test is performed with a cubic array of particles. In the test the simulation settings are equal to the settings of Deen et al. (2009) with a domain of 1603. The dimensionless drag force is determined for the central particle in the array. The resulting dimensionless drag force is 19.2. To compare the results with the results of Zick and Homsy (1982), Deen et al. (2009) made a fit of the data of Zick and Homsy (1982). The obtained drag force is within the 3% accuracy of the fit. Besides the comparison with other codes, the FT-IBM code, can also be compared to experiments. In these simulations, the terminal settling velocity of three particles is determined. The parameters used in these simulations are given in Table 3. The table also shows the experimental and the simulation results. The simulation results give similar results as the experiments. The table shows that the terminal velocity determined in the simulations is similar to the experimental results (Mordant and Pinton, 2000) and the terminal velocity calculated with the correlation of Clift et al. (1978). In the simulations with high Reynolds numbers, the results are influenced by the confinement of the wake in the simulation domain. Due to this confinement, the wake could not fully develop. This might explain the difference between the simulation results and the experimental and correlation results. To validate the FT part of the combined code, five single bubble simulations have been performed in different regimes. The used ¨ os ¨ and Morton numbers are given in Table 4. The resulting Eotv Reynolds numbers are also given in the table. The results of the simulations are similar to the results of Dijkhuizen et al. (2010a). The small differences are caused by slightly different positioning of the bubble in the domain (Dijkhuizen et al., 2010a places his bubbles at (0.5 0.5 0.6) domain size while in the simulations presented here the bubble is placed at (0.5 0.5 0.67) domain size). Furthermore, the code used by Dijkhuizen et al. (2010a) did not include volume conservative remeshing. Therefore, the bubbles decreased in size during the simulation. The difference between Tomiyama (1998) and the presented results can be caused by the slightly contaminated liquids used by Tomiyama (Dijkhuizen et al., 2010a). From these validation results, we conclude that the combined FT-IB model is able to simulate particles and bubbles correctly. 3.3. Three phase systems Following the verification and the validation of the FT-IB method, the method is used for the simulation of three phase systems using periodic domains. In total, 10 test cases were simulated using the combined model. The settings used in the simulations are shown in Table 5. The number of bubbles is fixed to 16. Roghair et al. (2011a) established that 16 bubbles are enough to simulate a bubble swarm using periodic boundary conditions. The size of the cubic domain was adjusted based on the gas volume fraction. Subsequently, the number of particles used in the simulations was determined from the size of the domain and the imposed particle volume fraction.
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
297
Table 3 The parameters used for settling sphere cases. The experimental results for the terminal rise velocity, vt, are taken from Mordant and Pinton (2000). Furthermore, the terminal velocity of Clift et al. (1978) is calculated with the correlation proposed by Clift et al. (1978). Case
1
2
3
Computational grid (W D H) Particle positions Grid size (m)
140 140 180 (70,70,40)
120 120 220 (60,60,40)
120 120 260 (60,60,40)
2:5 105
4:0 105
5:0 105
Time step (s)
1:0 105
1:0 105
1:0 105
Particle diameter (m)
5:0 104 (0,0, 9.81) 1000
8:0 104 (0,0, 9.81) 1000
1:0 103 (0,0, 9.81) 1000
Particle density (kg m3 )
9:0 104 2560
9:0 104 7710
9:1 104 7850
Experimental vt (m s1 ) vt of Clift et al. (1978) (m s1 ) Computed vt (m s1 )
0.074 0.076 0.076
0.316 0.325 0.312
0.383 0.393 0.372
g (m s2 ) Liquid density (kg m3 ) Liquid viscosity (Pa s)
Table 4 The parameters used in the simulation of single rising bubbles. The resulting Reynolds numbers are also shown in this table. For comparison, the experimental results of Tomiyama (1998) and the simulation results of Dijkhuizen et al. (2010a) are also shown. Case
1
2
3
4
5
Mo
1:0 103 1.0
1:0 103 10
1.0
Eo
100
1:0 103 100
1:0 1012 2.0
Re of Tomiyama (1998)
2.11
26.5
22.8
2.10
2:06 103
Re of Dijkhuizen et al. (2010a)
2.21
23.6
18.4
2.07
2:33 103
Computed Re
2.33
24.2
19.7
1.99
2:36 103
Table 5 The simulation settings for the three phase systems. The grid size of the cubic grid and the number of particle used in the simulations are determined by the volume fraction of gas and the volume fraction of particles used. Bubble diameter (m)
2:0 103
Particle diameter (m)
1:0 103
Grid size (m)
5:0 105
Time step (s)
1:0 105
Fluid density (kg m3 )
1:0 103
Gas density (kg m3 )
1:0 102
Particle density (kg m3 ) Gas viscosity (Pa s)
2:0 103
Fluid viscosity (Pa s) Number of bubbles Number of particles Bubble volume fraction (–) Particle volume fraction (–)
1:8 105 1:0 103 16 8y64 0.20, 0.30, 0.40 0.025, 0.05, 0.075, 0.10
force and F W the forces due to the walls mb
dvb ¼ F G þ F P þ F D þ F L þ F VM þ F W dt
Assuming a pseudo-steady state, the time-averaged force balance reduces to Eq. (8). In this equation, the virtual mass force is zero due to the steady state assumption and the lift force is zero when the velocity gradients are calculated over a length scale equal to the computational domain size. Furthermore, the walls are not taken into account in the DNS simulations. Therefore, the F W will be equal to zero /F D S ¼ /F G S þ/F P S
ð8Þ
1 2 C D rl pdb 9vb u9ðvuÞ ¼ V b rg g þ rpV b 8
ð9Þ
1 2 C D rl pdb 9vb u9ðvuÞ ¼ V b g ðrg þ ðarg þ frs þ ð1afÞrl ÞÞ 8 Fig. 4 shows two typical snapshots of the three phase simulations. Both configurations clearly show a horizontal alignment of the bubbles, see Fig. 4c and f. A similar type of clustering was also observed for bubble swarms by Roghair et al. (2011b). The horizontal alignment of the bubbles is probably caused by the lack of large scale circulation phenomena and the absence of walls in the simulations. Furthermore, due to the alignment of the bubbles, the particles also show preference for horizontal and vertical alignment (Fig. 4b and e). The results of the DNS simulations can provide a closure for the drag force in the Discrete Bubble Model (DBM) (i.e. the unresolved model using Lagrangian tracking of the dispersed elements). In the DBM, the velocity of the bubbles is calculated using Newton’s second law, Eq. (7). In this equation v b is the bubble velocity, F G the gravity force, F P the force due to pressure gradients, F D the drag force, F L the lift force, F VM the virtual mass
ð7Þ
ð10Þ
Insertion of the expressions for the definition of the drag force and the pressure force in Eq. (8) leads to Eq. (10). In this equation db is the bubble diameter, Vb the volume of the bubble, f the solids volume fraction and a the gas volume fraction. Rearranging Eq. (10), the relative drag coefficient can be determined, according to Eq. (11) (Roghair et al., 2011a). C D,1 and v b,1 in Eq. (11) represent the drag coefficient and the velocity of a single bubble in an infinite, initially quiescent liquid, respectively, as defined in Roghair et al. (2011a) C D,rel ¼
CD
rs rl rl rg
C D,1 ð1aÞ þ f
!¼
/v b,1 S2 ð/vuSÞ2
ð11Þ
The relative drag coefficient is determined for all 10 simulations. The results are presented in Fig. 5. The figure shows an increase in
298
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
Fig. 4. Two snapshots (a) and (d) of a three phase simulation containing 30% bubbles and 7.5% particles. (b) and (e) show the configuration of the particles in the simulation, while (c) and (f) the configuration of the bubbles.
20
2.5
CD,rel three phases/CD, rel two phases
α = 20 % α = 30 % α = 40 %
CD, rel (−)
15
10
5
0
2
4
6
8
10
φ (%)
α = 20 % α = 30 % α = 40 %
2
1.5
1
0
2
4
6
8
10
φ (%)
Fig. 5. The effect of the solids fraction (f) on the relative drag coefficient. The line around the point a ¼ 40% and f ¼ 2:5% shows the standard deviation of the determined relative drag coefficient.
Fig. 6. The effect of the solids fraction (f) on the relative drag coefficient over the relative drag coefficient of a bubble swarm with the same void fraction. The dashed line is the trend line.
the drag coefficient with an increase in the solids volume fraction. This can be explained by the particles that obstruct the bubbles during their rise. Fig. 5 also shows the standard deviation of the simulation with an a of 40% and a f of 2.5%. The value for the standard deviation is typical for all the simulations. The relative large standard deviation can be explained by the short averaging periods used in these simulations (0.1–0.35 s), however, the trend is
very clear. The development of a fully predictive correlation will be the subject of a future contribution. Therefore, the relative drag coefficient is divided by the relative drag coefficient of a bubble swarm (which can be calculated using Eq. (11) with a f of zero). The results are shown in Fig. 6. All the data points seem to coincide with the trend line shown in this figure.
M.W. Baltussen et al. / Chemical Engineering Science 100 (2013) 293–299
299
4. Conclusions
Acknowledgments
In this paper, a new model is presented for the Direct Numerical Simulation of three phase flows. The model combines the FT model of Dijkhuizen et al. (2010b) and the IB model of Kriebitzsch (2011). The combined model is verified and validated for both the FT and the IB part of the model. The FT-IB model was used to simulate dense threephase flows. The relative drag coefficient of the bubbles in the simulation of the three-phase flow seems to increase with an increasing solids fraction. However, longer averaging periods are needed to derive a fully predictive correlation for the relative drag coefficient with respect to the solid volume fraction.
The authors would like to thank the European Research Council for its financial support, under its Starting Investigator Grant scheme, contract number 259521 (CuttingBubbles).
Nomenclature Roman symbols C d F g I m p r T t u v w V
Coefficient Diameter, m Force, N Gravity constant, m s 2 Moment of inertia, kg m2 Mass, kg Pressure, Pa Position, m Torque, m N Time, s Liquid, fluid velocity, m s 1 Particle velocity, m s 1 Translation velocity, m s 1 Volume, m3
Greek symbols
a m f
r s o Abbreviations and subscripts b D DBM DNS FT g G IB i, j, k l L m Mo Eo p P Re s t VM W
s 1
Gas volume fraction Viscosity, Pa s Solid volume fraction Density, kg m 3 Stress tensor Rotational velocity, s 1 Bubble Drag Discrete Bubble Model Direct Numerical Simulations Front Tracking Gas phase Gravity Immersed boundary Center of mass Liquid phase Lift Marker Morton number ¨ os ¨ number Eotv Particle Pressure Reynolds number Solid Terminal Virtual mass Wall Surface tension Single bubble infinite liquid
References Beetstra, R., Van Der Hoef, M.A., Kuipers, J.A.M., 2007. Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres. AIChE J. 53 (2), 489–501. Clift, R., Grace, J., Weber, M., 1978. Bubbles, Drops and Particles. Academic Press, New York. Deen, N.G., Mudde, R., Kuipers, J.A.M., Zehner, P., Kraume, M., 2010. Bubble columns. In: Ullmann’s Encyclopedia of Industrial Chemistry. Wiley-VCH Verlag GmbH & Co, KGaA, Weinheim. Deen, N.G., van Sint Annaland, M., Kuipers, J.A.M., 2009. Direct numerical simulation of complex multi-fluid flows using a combined front tracking and immersed boundary method. Chem. Eng. Sci. 64 (9), 2186–2201. Dijkhuizen, W., 2008. Derivation of Closures for Bubbly Flows Using Direct Numerical Simulations. Ph.D. Thesis, University of Twente. Dijkhuizen, W., Roghair, I., Van Sint Annaland, M., Kuipers, J.A.M., 2010a. DNS of gas bubbles behaviour using an improved 3D front tracking model—drag force on isolated bubbles and comparison with experiments. Chem. Eng. Sci. 65 (4), 1415–1426. Dijkhuizen, W., Roghair, I., Van-Sint-Annaland, M., Kuipers, J.A.M., 2010b. DNS of gas bubbles behaviour using an improved 3D front tracking model—model development. Chem. Eng. Sci. 65 (4), 1427–1437. Ge, Y., Fan, L.-S., 2006. 3-d direct numerical simulation of gas–liquid and gas–liquid– solid flow systems using the level-set and immersed-boundary methods. Adv. Chem. Eng. 31, 1–63. Hasimoto, H., 1959. On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres. J. Fluid Mech. 5 (02), 317–328. Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J., Van Swaaij, W.P.M., 1996. Discrete particle simulation of bubble and slug formation in a two-dimensional gasfluidised bed: a hard-sphere approach. Chem. Eng. Sci. 51 (1), 99–118. Jain, D., Deen, N.G., Kuipers, J.A.M., Antonyuk, S., Heinrich, S., 2012. Direct numerical simulation of particle impact on thin liquid films using a combined volume of fluid and immersed boundary method. Chem. Eng. Sci. 69 (1), 530–540. Kantarci, N., Borak, F., Ulgen, K., 2005. Bubble column reactors. Process Biochem. 40 (7), 2263–2283. Kriebitzsch, S., 2011. Direct Numerical Simulation of Dense Gas–Solid Flows. Ph.D. Thesis, Eindhoven, University of Technology. Kuprat, A., Khamayseh, A., George, D., Larkey, L., 2001. Volume conserving smoothing for piecewise linear curves, surfaces, and triple lines. J. Comput. Phys. 172, 99–118. Lamb, H., 1932. Hydrodynamics, 6th ed. Cambridge University Press, London. Li, Y., Yang, G., Zhang, J., Fan, L., 2001. Numerical studies of bubble formation dynamics in gas–liquid–solid fluidization at high pressures. Powder Technol. 116 (2–3), 246–260. Mordant, N., Pinton, J.-F., 2000. Velocity measurement of a settling sphere. Eur. Phys. J. B 18 (2), 343–352. Rider, W.J., Kothe, D.B., 1998. Reconstructing volume tracking. J. Comput. Phys. 141, 112. Roghair, I., 2012. Direct Numerical Simulations of Hydrodynamics and Mass Transfer in Dense Bubbly Flows. Ph.D. Thesis, Eindhoven, University of Technology. Roghair, I., Lau, Y., Deen, N.G., Slagter, H., Baltussen, M., Annaland, M.V.S., Kuipers, J.A.M., 2011a. On the drag force of bubbles in bubble swarms at intermediate and high Reynolds numbers. Chem. Eng. Sci. 66, 3204–3211. Roghair, I., Martı´nez Mercado, J., Annaland, M.V.S., Kuipers, H., Sun, C., Lohse, D., 2011b. Energy spectra and bubble velocity distributions in pseudo-turbulence: numerical simulations vs. experiments. Int. J. Multiphase Flow 37 (9), 1093–1098. Tomiyama, A., 1998. Struggle with computational bubble dynamics. In: Third International Conference on Multiphase Flow, pp. 369–405. van Sint Annaland, M., Deen, N.G., Kuipers, J.A.M., 2005. Numerical simulation of gas bubbles behaviour using a three-dimensional volume of fluid method. Chem. Eng. Sci. 60 (11), 2999–3011. Wang, T., Wang, J., Jin, Y., 2007. Slurry reactors for gas-to-liquid processes: a review. Ind. Eng. Chem. Res. 46 (18), 5824–5847. Yang, G., Du, B., Fan, L., 2007. Bubble formation and dynamics in gas–liquid–solid fluidization—a review. Chem. Eng. Sci. 62 (1–2), 2–27. Zick, A.A., Homsy, G.M., 1982. Stokes flow through periodic arrays of spheres. J. Fluid Mech. 115, 13–26.