Stochastic DSMC method for dense bubbly flows Methodology

(cid:1) Extension of the Direct Simulation Monte Carlo (DSMC) method for dense bubbly ﬂows. (cid:1) Veriﬁcation and validation with DPM/DBM and experiments from literature. (cid:1) Increased computational performance up to two orders of magnitude. A stochastic Direct Simulation Monte Carlo (DSMC) method has been extended for handling bubble- bubble and bubble-wall collisions. Bubbly ﬂows are generally characterized by highly correlated velocities due to presence of the surrounding liquid. The DSMC method has been improved to account for these kind of correlated collisions along with a treatment allowing the method to be used also at relatively high volume fractions. The method is ﬁrst veriﬁed with the deterministic Discrete Particle/Bubble Model (DPM/DBM) using two problem cases: (a) dry granular ﬂow of particles through two impinging nozzles and (b) 3D periodic bubble rise for mono-disperse and poly-disperse systems. The veriﬁcation parameters are the total number of prevailing collisions within the system, the collision frequencies and the time-averaged liquid velocity proﬁles (only for the 3D-periodic bubble rise). Subsequently the method is applied to a lab-scale bubble column and validated with the experimental data of Deen et al. (2001). A computational performance comparison with the DBM is reported for the 3D periodic bubble rise case with varying overall gas fractions. The DSMC is approximately two orders of magnitude faster than the deterministic approach for the studied dense bubbly ﬂow cases without adverse effects on the quality of the computational results.


Introduction
Bubbly flow is one of the most widely used methods for gasliquid contacting operations in the process industries. Processes like absorption, fermentation, Fischer-Tropsch synthesis, waste water treatment and bio-reaction based processes are typically operated in bubble columns. These operations are generally mass transfer limited and to obtain reasonable yields, certain criteria have to be satisfied. These include a high gas-liquid inter-facial area and fast mixing to enhance the reactor performance. Both criteria are met in bubble columns with the additional advantage of no moving mechanical parts. Industrial columns possess volumes in the range of 100-300 m 3 . Even larger columns are employed for bio-processes like fermentation (3000 m 3 ) and waste water treatment (20,000 m 3 ) (Deen et al., 2012).
The bubble dynamics inside the column dictate the (induced) liquid flow where coalescence, break-up and mass transfer rates determine the bubble size distribution. Local mass transfer rates are highly influenced by the mixing within the liquid phase which in turn depends on the previously mentioned phenomena in the column. Due to this strong coupling between the parameters, modeling of such bubble columns is very challenging and inevitably empirical input information (i.e. correlations) is required. A review of such correlations and the scale-up of these columns can be found in the works of Chaudhari and Ramachandran (1980), Shaikh and Al-Dahhan (2013). Due to recent advances in computational power and techniques, considerable progress has been made in simulation of the flow structures in such columns which occur at multiple time and length scales.
The accurate simulation of bubble columns defines a complex multi-scale problem due to the wide range of length scales that needs to be covered. Methods such as Front-Tracking (Unverdi and Tryggvason, 1992) and Volume of Fluid (van Sint Annaland et al., 2005) methods have been applied to simulate the behavior of a single bubble or a swarm of interacting bubbles rising through the liquid phase. Closures obtained from these simulations have been developed and implemented into methods used at larger length and time scales, such as the Discrete Bubble Model (DBM) and the Multi Fluid Models (MFM). The DBM constitutes a suitable framework to understand the influence of bubble-bubble interactions and bubble wall encounters on large scale flow structures. There are several recent works which have used this framework to simulate bubbly flows Jain et al., 2014) using the hard-sphere model for bubble-bubble and bubble-wall interactions. Several recent studies have also incorporated a soft sphere model which is much more parallelizable for bubblebubble and bubble-wall collisions (Heitkam et al., 2017;Xue et al., 2017a,b;Lau et al., 2014). However, the method becomes computationally very expensive when simulating dense bubbly flows because every bubble is tracked explicitly and the bubblebubble encounters are resolved in time. Therefore, at the industrial scale, Euler-Euler approaches like the MFMs are preferred to save time (Deen et al., 2004).
Euler-Euler approaches such as the Multi Fluid Model represent the bubble phase with different size classes. Coupled momentum and mass conservation equations are solved for each size class with suitable closures. The interactions among the different considered phases need to be defined explicitly. Bubble-bubble interactions are not explicitly accounted for nor the naturally occurring phenomena of coalescence and break-up. Dedicated external models, such as population balance models, have to be integrated with the Multi-Fluid Model to account for coalescence and break-up. Some of the recent works in this field can be found in Zhang (2007), Fletcher et al. (2017), Giorgio et al. (2017), Buwa and Ranade (2002), Bhole et al. (2008), Yang and Xiao (2017). In the work of Zhang, 2007 the range of bubble sizes is assumed to be 2-4 mm while the other parameters such as drag, lift etc. are assumed to not vary significantly. Therefore a two-fluid approach is employed. Similarly, Krishna et al. (2000) have considered three phases, one for large bubbles, one for small bubbles and one for the liquid phase.
Stochastic Euler-Lagrangian methods provide interesting compromise with respect to consideration of both, the computational efficiency and the accurate representation of the discrete phase. There are several examples of stochastic methods being used for particulate flows in the literature (O'Rourke, 1981;Du et al., 2011;Pawar et al., 2014;Sungkorn et al., 2011). These methods have also been coupled with Eulerian methods to simulate multiphase flows. Pawar et al. (2015) have simulated droplet sprays using the Direct Simulation Monte Carlo (DSMC) method for droplet-droplet collisions. Sommerfeld (2001) has presented a fictitious ghost particle approach for the inter-particle collisions. This approach considers the velocity fluctuations in the particle phase due to the turbulence in the fluid phase. It has also been used by Witz et al. (2016) for the execution of bubble-bubble interactions in an industrial scale bio-reactor. Huilin et al. (2006) have simulated a 2D fluidized bed using the DSMC approach also incorporating the DSMC parcels; however the porosity calculation for the coupling with the fluid does not seem to be consistent since the real particle positions in a parcel are not known.
The main objective of this work is to develop a computationally inexpensive Euler-Lagrange model that can be used for large scale dense bubbly flows. Bubble-bubble and bubble-wall interactions are critical bottlenecks in the computational speed of deterministic simulations. To resolve this issue, a stochastic DSMC approach is used to execute the collisions in a more efficient manner. Pawar et al. (2014) and Du et al. (2011) have shown that the method is nearly two orders of magnitude faster than the DBM, but limitations to the approach were encountered with respect to the volume fraction of the discrete phase. This work also aims to extend the DSMC approach to alleviate the volume fraction limitation observed by Pawar et al. (2014) and improve the DSMC algorithm for particles/bubbles with highly correlated velocities.
A detailed methodology of the DSMC method extended to dense bubbly flow application is presented in this paper. Section 2 describes the method used, the modifications and the algorithm. Verification of the method is performed for two systems: a) Dry granular flow b) 3D periodic bubble rise and will be reported in Section 4. The verification model chosen is the Discrete Bubble Model . Furthermore, a validation with experimental data reported in the literature is performed through the simulation of a lab-scale bubble column (Deen et al., 2001). Lastly the computational performance of the method is investi-

Nomenclature
Roman symbols F force (N) g radial distribution function at contact (-) r position vector of the discrete phase (particles/bubbles) P ij probability of collision between particle ids (-) Stk Stokes number ( gated on a single core and compared with the performance of the Discrete Bubble Model.

Discrete phase
Bubble sizes depend on the regime in which a bubble column is operated as well as the characteristics of the gas distributor. The bubble size is typically in the range of a few millimeters to a few centimeters (Deen et al., 2012). The bubble motion is governed by the Newtonian equations of motion (see Eq. (1)) where the considered forces include here are gravitational (F G ), pressure ðF P Þ, drag ðF D Þ, lift ðF L Þ, virtual mass ðF VM Þ and wall lubrication ðF W Þ force (see Eq. (2)). The relevant closures used for these forces can be found in Table 1.
The closures used for drag, lift and wall forces used here are from Tomiyama et al. (1995Tomiyama et al. ( , 2002. The drag correction for dense bubble swarms that is included here can be found in Roghair et al. (2011). A more detailed explanation about the forces and closures can be found in works of Darmana et al. (2006), Jain et al. (2014) and Delnoij et al. (1997). The DSMC methods generally employ a parcel approach where a group of discrete particles is represented by one simulated particle. The particles in a parcel are considered to be of the same size and to move with the same velocity as the simulated particle. Velocity change only occurs due to collisions with other parcels. This approach is adopted to save computational costs and is based on the assumption that the parcels grow isotropically purely due to diffusion in time. The systems for which this method was initially developed and applied to are molecular systems that have finite/ large Knudsen number (Bird, 1994). In this work a parcel size of 1 is chosen to keep things simple and relatively more accurate. For a parcel size of 1, Pawar et al. (2014) have shown a major (60 to 70 times) speed boost compared to the deterministic Discrete Particle approach.
Extension to larger parcel sizes requires careful consideration of different forces for the bubbles as a swarm and the mapping functions for calculation of the porosity. Moreover, the trajectory of the bubbles is non-linear due to coupling with the liquid and the parcel needs to grow anisotropically in the direction of buoyancy with time. This will be considered in our future work.

Liquid phase hydrodynamics
The liquid phase hydrodynamics is obtained by solving the volume-averaged Navier-Stokes equations coupled with the continuity equation (see Eqs. (4) and (3)). @ðe l q l Þ @t þ r Á e l q l u ¼ 0 ð3Þ @ðe l q l uÞ @t þ r Á e l q l uu ¼ Àe l rP À r Á e l s l þ e l q l g þ U where e l is the local liquid fraction. U represents the local momentum source term for the force exerted by the bubbles on the liquid.
The shear stress term s l is given by: Here, l L;l is the dynamic viscosity and l T;l is the turbulent viscosity. The contribution with l T;l originates from the convection term when the Navier-Stokes equations are volume-averaged and the Boussinesq eddy viscosity assumption is applied. A large volume of literature exists comparing different turbulence models with experimental data in a bubble column (Tabib et al., 2008;Deen et al., 2001). They concluded that the large eddy simulations outperform the k À e model. LES was chosen for the explicit calculation of the eddy viscosity (l T;l ) in the liquid phase. The l T;l is closed with the sub-grid scale eddy viscosity expression given by Vreman (2004).
where c ¼ 2:5C 2 s ; C s is the Smagorinsky constant, 23 . D m is the size of the grid in the m direction. A similar approach has been used by Darmana et al. (2006) and also by Jain et al. (2014).

Time marching
The numerical solution of the governing equations for our Euler-Lagrange model involves time-marching where the equations for the discrete and continuous phase are solved sequentially Table 1 Closures for the different types of forces acting on bubbles.

Force
Closure Reference  Tomiyama et al. (1995Tomiyama et al. ( , 2002, Roghair et al. (2011) F L ¼ ÀC L q l V b ðv À uÞ Â ðr Â uÞ  Tomiyama et al. (1995Tomiyama et al. ( , 2002 Tomiyama et al. (1995,2002) while accounting for the coupling. The flow time step (dt flow ) is divided into a fixed number of bubble time steps (dt bub ). The forces on each bubble are calculated based on the Eulerian flow-field around them for every bubble time-step dt bub . dt bub may be divided into several more time steps based on the collision frequency calculated for every bubble. More details about this procedure are given in Section 3.3. The major difference in time marching from the deterministic approach employed by Darmana et al. (2006) and the stochastic approach used in this work is the resolution of the lowest time scale. In the deterministic approach, the collision times are explicitly calculated and a list is maintained to sequentially do the collisions in the ascending order of their times. In the current approach, the collision times are not determined. The collisions take place based on a probabilistic approach.

Coupling
Inter-phase coupling consists of porosity mapping and the computation of the volumetric momentum exchange rate between the Euler and Lagrangian frameworks. Kitagawa et al. (2001) have systematically laid down the rules for a Lagrangian template function for mapping the porosity based on the discrete phase positions. Along the same lines, Deen et al. (2004) have proposed a fourth order polynomial filter for mapping. The filter width is a userdefined variable. Darmana et al. (2006) have used a variable filter width based on the bubble size for their Discrete Bubble Model which makes it independent of the Eulerian grid size. This mapping function is used in this work for the inter-phase coupling. Lau et al. (2011) have reported that the effect of changing the filter width on the predicted time averaged axial liquid velocities and velocity fluctuations in a bubble column is negligible when a drag correction based on local porosity is used. In this work, a filter width size of 2R bub is used for a one-to-one comparison between both the collision models.

Discrete phase dynamics and collision algorithm
The bubble phase equation of motion is solved using a first order explicit scheme for Eq. (1): The forces are mapped to the discrete bubble locations from the Eulerian grid cells using the polynomial filter described in Section 3.2. The collision sequence is then initiated according to the algorithm shown in Fig. 1 and Eq. (8). This also includes the update of bubble positions. Forces on the bubbles are recalculated every bubble time-step (dt bub ) and are also collected in the volumetric momentum source term U with an opposite sign. This is repeated until the discrete phase has moved for a full dt flow . The forces collected in the volumetric source term are time-averaged over the flow time-step dt flow , since they are calculated multiple times in one flow time-step. Weights for the averaging are calculated based on the ratio dt bub dt flow . With the new bubble positions, the volume fraction in each Eulerian cell is determined using the polynomial mapping function. The collisions that occur during a bubble time-step can cause a big computational overhead for a dense, large scale system with bubble/particle numbers in the order of millions. The algorithm of a deterministic model such as the Discrete Bubble Model  is inherently serial even with a domain decomposition as the collisions need to occur in the order of their collision times. In a parallel environment this time needs to be communicated every time the encounter list is updated. These col-lision times are maintained in an encounter list from which, the minimum time needs to be searched every time a collision occurs due to changes in the bubble positions and velocities. This is computationally expensive for dense systems, even with the implementation of neighbour lists.
A speed-up with the DSMC algorithm compared to the DBM is obtained when the collision partner is chosen using probabilistic rules; thereby avoiding the necessity of an encounter list. There are several methods used for choosing the right collision pairs in DSMC such as the time counter method, Nanbu method, modified Nanbu method and many more (Lutišan, 1995). Pawar et al. (2014) and Du et al. (2011) have used the modified Nanbu method with their own modifications for treatment of the searching scope and collision conditions.
A detailed description and motivation behind the modifications to the algorithm in this work will be given after the description of the algorithm. The algorithm is described as follows: 1. Choose a bubble/particle id i from a randomly generated list of id's. 2. Calculate the collision frequency f i for i based on Eq. (9).
where v ij is the relative velocity between particles/bubbles i and j; d is the diameter, n j is the parcel size, R s;i is the searching scope size for particle/bubble i and g ij is the radial distribution function at contact for discrete entity i with particle type j. During this phase, the nearest neighbour to i based on surface to surface distance is determined, as well as a local list of N i neighbour ids of particle/bubbles within the range of R s;i . If the number of neighbours is less than N min , the searching scope is enlarged and step 2 is repeated. 3. The discrete phase time-step Dt p;i is found from the mean free path (k i ) and the magnitude of the velocity of the particle/bubble (v i ). If the remnant time (dt bub À Dt compl ) is smaller than the discrete phase time step, then the remnant time is taken into account.
4. Calculate new searching scope for the next discrete phase timestep of particle/bubble i (see Eq. (12)).
is chosen between 0 and 1. The particle j is chosen based on Eq. (13). The collision probability of i with j is calculated using Eq. (14).
6. The following conditions (see Eqs. (15) and (16)) are to be satisfied if a collision has to occur. Eq. (16) makes sure that the discrete phase entities are moving towards each other which leads to more realistic collisions (Pawar et al., 2014). 7. This step is executed only for the cases with fluid coupling. If a collision in step 6 is rejected then the bubble Stokes number (see Eq. (17)) is checked as per Algorithm 1. If the condition is satisfied then it is checked if the nearest neighbour of i will overlap with i. If this is true, then i is collided with its nearest neighbour.
Algorithm 1. Explicit check for bubbly flow systems.
if Stk bub < 1 then check for collision with nearest neighbour of i else no collisions; move to next step end if where s bub is approximated for a spherical bubble and the nonlinear part of the drag is assumed to be a constant for the current time-step. The liquid phase relaxation time is estimated from the turbulent dissipation calculated from the implemented LES model.
where m l is the liquid kinematic viscosity and is the dissipation rate of turbulent kinetic energy. ¼< 2m e jjS ij jj 2 > where m e is the eddy viscosity calculated from the sub-grid scale turbulence model and S ij is the rate of deformation tensor. 8. Check for wall collisions. If a particle/bubble overlaps with a wall, then it is bounced back. 9. Update position (r i ), searching scope (R s;i ) of i and time completed (Dt compl ). Check for time step completion It should be noted that the collision frequency (f i ) and probability term (P ij ) include the radial distribution function (g ij ). As the flow becomes denser, a structural factor becomes important to determine the correct collision frequency due to a non-uniform spatial distribution of the particles. In a packed or a highly clustered system, the local relative velocities are low which leads to low collision frequencies and probabilities based on the equation proposed by Bird (1994). In these cases the clusters travel together with continuous collisions occurring among them. DSMC can predict the right collision frequency as long as the velocities of the particles/bubbles remain uncorrelated and the non-uniform spatial distribution is accounted for via g ij (Sundaram and Collins, 1997).

Radial distribution function (g ij )
The expression for g ij is closed from Ma and Ahmadi (1986) for a mono-disperse system (see Eq. (19)). The extension to polydisperse systems is based on the work of Santos et al. (1999). The radial distribution function for a multi-component hardsphere mixture can be approximated from the compressibility of a single component hard-sphere system (see Eq. (20)).
where d represents the sauter mean diameter of the particles within the searching scope whereas e p is the average solids/bubble volume fraction in the neighborhood of particle i.

Interstitial liquid based velocity correlations
Clustering and layering in bubbly flow is a common phenomenon observed in DNS  as well as experiments (Van Wijngaarden, 1993). Bubble dynamics is strongly coupled with the bubble-induced turbulence especially for a bub- Fig. 2. Schematic of the systems in P1 and P2.

Table 2
Parameter space for the mono-disperse system of particle size 2 mm applied to P1 (impinging particle streams).

Case
Mean   Table 2, DPM vs DSMC.  Table 2 applied to P1 for a mono-disperse system with particle size 2 mm, DBM vs DSMC.
ble in the wake of another bubble or a bubble swarm. As a consequence, the bubble velocities in and near the same region in a liquid are correlated (or biased). This can be quantified in terms of the bubble Stokes number. The seminal paper of Sundaram and Collins (1997), reported the behavior of particles with different Stokes numbers in the dilute limit in a turbulent suspension. As the particle Stokes number increases, the inertia reduces the degree of correlation in velocities with the fluid. This leads to a more complex behavior with partly correlated and partly uncorrelated velocities of particles/bubbles. There exists no theoretical basis in predicting collision frequency for this intermediate regime. Therefore an explicit treatment with the nearest neighbour collisions is added to the algorithm when there is fluid coupling and a collision is rejected due to highly correlated velocities (or low Stokes number). As the slip velocity of bubbles increase, the bubble Stokes number increases therefore the condition is based on an approximation of the calculated bubble Stokes number which in turn is based on the experienced fluid drag.
DSMC simulations are typically initiated with a dynamic system and hence collisions occur because of its dynamic nature. With the addition of the above Stokes condition the discrete phase can be started from rest, as from rest most of the collisions occur because of the fluid-bubble/particle interactions. This is also shown in Section 4.

Liquid phase numerical scheme
The continuity equation and the momentum equations are solved in a coupled manner using SIMPLE (Patankar, 1980). The unsteady part is discretized using the first order Euler scheme. The convective fluxes in the continuity equation are treated implic-itly. The convective, diffusive and the source terms are treated explicitly in the momentum equations.

Verification
The Discrete Bubble Model from Darmana et al. (2006) is taken as verification model for the DSMC. This is done for two distinctly different systems: 1. Impinging particle streams (P1): For the verification of the DSMC method for larger particle/bubble sizes with a nonuniform spatial distribution. 2. 3D -periodic bubble rise from rest (P2): For the verification in a bulk bubbly flow system with 4 way coupling.
where P1 corresponds to the case without the presence of any interstitial fluid and P2 to the bubbly flow case (with interstitial liquid). P1 consists of two nozzles facing each other at an angle (see Fig. 2a). The solid particles enter through these nozzles, collide at the impact region and spread across the domain based on the outcome of the collisions. The velocity distribution purely depends on the collision statistics within the system. Absence of a fluid phase also means there is no correlation in particle velocities. The mass flow-rate is set for each nozzle. The parameters that are varied, are the mean particle velocities and the standard deviation of the velocities. A Gaussian distribution is used to generate velocities in the x and z direction with given mean velocity and standard deviation. The mean velocity in the y direction is set to 0 and perturbations are imposed based on the same Gaussian distribution.
P2 represents a 3D periodic system which starts as an array of regularly spaced bubbles (see Fig. 2b). The spacing between bub-  Table 2 applied to P1 for a mono-disperse system with particle size 2 mm, DBM vs DSMC.

Table 3
Parameter space for the mono-disperse system of different bubble sizes applied to P2 (3D periodic bubble rise).

Case
Bubble size ( bles across the periodic boundaries is a little larger than in the central bulk of the domain. Thus, a concentration bias exists in the bulk of the domain. The bubbles start from rest and their motion is coupled with the liquid phase motion. The physical properties used in this system correspond to that of the air-water system. The following sections are classified based on particle/bubble system properties. The measurement parameters are the total number of collisions and the average collision frequency as a function of time. 4.1. Mono-disperse systems 4.1.1. Impinging particle streams (P1) The probability ðP ij Þ depends on the relative velocity distribution. Therefore it is necessary to test the DSMC model at different velocity regimes and different solids fractions at the impact region. The parametric space is defined based on the mean velocity of the particle streams and the standard deviation of the inlet velocity distribution (see Table 2).
It can be clearly observed from Figs. 3 and 5 that the DSMC can handle collisions even at high dispersed phase fractions. High solids volume fractions occur in cases 1 and 3 due to a low standard deviation in the inlet velocity distribution. Equivalently the total number of collisions for cases 2 and 4 have a maximum error (based on the difference between DSMC and DPM) of 3 %. The initial evolution of the collision frequency within the first second after the particles enter the system is also captured reasonably well by DSMC in all cases. Visually, the DSMC controlled collisions lead to a similar particle concentration throughout the domain compared to their DPM counterparts (see Fig. 3). Quantitatively, the average collision frequency at steady state for different cases is compared for both methods in Fig. 4.
It is known for a fact that in particulate systems larger velocity deviations typically lead to more collisions which seems counterintuitive from the obtained results. However, for this to be true, the systems must be loaded at a comparable solid volume fraction. In the tests shown for different cases here, the smaller standard deviation in velocities leads to two very focused particle streams which create a small impact region with high solids fraction. This in turn leads to a higher number of collisions and collision frequencies compared to the situation with high standard deviation in the input velocity where the particle streams spread out from the onset to give a larger impact region with lower solids fraction. Similar explanation is applicable for test cases shown in Section 4.2.1. Fig. 6 shows the particle velocity probability density functions for the DBM and DSMC. The particle velocities are sampled at a frequency equal to the inverse mean residence time, calculated based on the average velocity of the particles. The particles taken into consideration during sampling are below the impact region to account for post collisional velocities. It can be seen that on an average there is no difference in the collision outcomes from DPM and DSMC.

3D-periodic bubble rise (P2)
For the 3D-periodic bubble rise (P2) two systems are simulated with different bubble diameters (see Table 3). The velocity of the system changes with time as the bubbles rise. They act as small momentum sources meaning that the periodic system is evolving from lower to higher bubble and liquid velocity fields.
The bubbles start from rest and are arranged in a periodic array. As the bubbles rise, the bubbles in the wake of other bubbles accelerate and the initial structure disappears. Fig. 7 shows the time evolving structure of bubbles as they rise through the periodic sys-  Table 3 applied to P2 for a monodisperse system with bubble size 2 mm, DBM vs DSMC. tem. The system starts with layers. As the momentum of the system increases these layers reduce to small clusters due to increased (pseudo) turbulence in the liquid. Finally, these clusters disintegrate to a completely homogeneous system at high slip velocities. These features are equivalently captured by the DSMC algorithm. The transitioning of the bubble rise velocity field with time can be observed from the velocity scales in Fig. 7.
Evolution for the 4 mm system is not shown here but the qualitative comparison between both methods reveals that both methods perform equally well. The bubbles reach high velocities much faster and the system is also relatively more dense due to the use of larger bubble sizes with the same number of bubbles in a similarly sized domain.
The total number of collisions and the collision frequency vs time for the 2 mm and 4 mm (Case 1 and Case 2) systems can be compared on basis of Figs. 8 and 9. It can be observed that the col-  Table 3 applied to P2 for a monodisperse system with bubble size 4 mm, DBM vs DSMC.  Table 3 applied to P2, DBM vs DSMC.

Table 4
Parameter space for the poly-disperse system of particle size range 1-6 mm applied to P1 (impinging particle streams).

Case
Mean  Fig. 11. Steady state domain snapshots for poly-disperse dry granular streams, DPM vs DSMC. Fig. 12. Steady state collision frequency snapshots in the xz plane for a cross section at y YMAX ¼ 0:5 for poly-disperse system cases in Table 4, DPM vs DSMC.  Table 4 applied to P1 for a poly-disperse system with particle size range of 1-6 mm, DBM vs DSMC.
lision frequencies in both 2 mm and 4 mm cases match reasonably well between DBM and DSMC. At high velocity regimes (in the order of several m/s) in the 2 mm case the collision frequency is slightly over-predicted. At the same time, it can be observed from Fig. 10 that the liquid velocities are not over-predicted by the modified DSMC method. This shows that the current method is executing the correct collisions along with some extra ones in the direction lateral to the flow. Pure DSMC method without the currently made modifications can predict the right collision frequency at very high bubble velocities but cannot generate the right bubble velocities. This is because the bubble velocities are not completely uncorrelated as in the case of P1. Some correlation with the liquid velocity field always exists, which means the explicit treatment is required for prediction of correct bubble velocities. It is also to be noted that such bubble velocities in the order of several meters per second are not seen in bubbly flows. This over-prediction of collisions is not seen in the 4 mm case even at higher velocity regimes because this is a denser system, where collisions in the direction lateral to the mean flow direction occur much more frequently due to lack of void space. Therefore the collision frequency is accurately predicted by the proposed method.
Time averaged liquid velocity profiles in the rise direction are taken along a line on the x axis, at y ¼ 0:5Y MAX and z ¼ 0:56Z MAX where Y MAX and Z MAX are max domain lengths in the respective directions (see Fig. 10). The maximum difference in each case is less than 6% relative to the equivalent DBM cases. The extrema in the velocity profiles in 2 mm and 4 mm cases are also captured by the modified DSMC method. In general, it can be observed that the liquid velocities in the DSMC cases are slightly underpredicted; this can be attributed to the fact that DSMC collisions are not exact. In bubbly flow, the collisions are mostly biased in the rise direction. Because of this, the bubbles on the top can attain even higher velocities whereas DSMC homogenizes these collisions isotropically which in turn averages the momenta in the searching scope around the bubble.

Poly-disperse systems
The newly developed algorithm is also verified for a set of cases in P1 and P2 for different poly-disperse systems. The particle size range in P1 is taken 1-6 mm. A similar size range is taken for P2 specified in Table 5.

Impinging particle streams (P1)
The parametric space is defined similar to the case in Section 4.1.1 (see Table 4). The particles are generated with a random size distribution in the range of 1-6 mm. The generated (or injected) particle is kept equivalent in both nozzles at a given time to preserve the same flowrate from both nozzles.
It can be visually observed from Fig. 11 that the DSMC method executes collisions in such a way that the average outcome of the simulation is similar to the DPM. This can be said from the corresponding particle spread in the system for both the methods for all cases defined in the parametric space. This is also verified quantitatively in Fig. 12 where the steady state collision frequency is compared between both DPM and DSMC for all the cases.
The velocity probability density functions measured in time also indicate that the average velocity distribution is nearly identical for both the methods in all cases (see Fig. 14).
The total number of collisions and collision frequency vs time can be seen in Fig. 13. The variation of the collision frequency as well as the total collisions are within a difference of 3%.

3D-periodic bubble rise (P2)
The parametric space in this section (see Table 5) is based on the number of bubbles in the box and the overall gas fraction. It should be noted that the overall gas fraction is calculated with respect to the whole simulation domain where as the local gas fractions reach up to 37% during the simulation.
The time based evolution of the bubble positions and flow structures are compared in Fig. 15. The bubbles are initially arranged sequentially in the following order of their sizes 1,2,3 and 4 mm (see in Fig. 15 at time t ¼ 0:018 s). The smaller bubbles have lower terminal rise velocity. Therefore, due to the initial alternate arrangement of the bubbles, the collision frequency is high until the system homogenizes. Due to displacement, the 1 mm bubbles occupy the interstitial space between the larger bubbles (see Fig. 15 at time t ¼ 0:237 s). Unlike the mono-disperse case, the bubbles undergo two successive stages of velocity core formation in the central region of the domain at times t ¼ 1 s and t ¼ 3:6 s. At t ¼ 1 s, a high velocity core is formed in the central region of  Table 4 applied to P1 for a poly-disperse system with particle size range of 1-6 mm, DBM vs DSMC. the domain, at t ¼ 2:7 s the velocities homogenize in the system. With increasing time a low velocity core is formed in the same region at t ¼ 3:6 s. Then the system homogenizes at t ¼ 4:2 s and the system momentum increases uniformly with random patches  Table 5) applied to P2 (Note: The bubbles are colored based on their velocity magnitudes in m/s; scale given on the right). of high and low velocity regions. This whole process is identically captured by the DSMC algorithm.
It can be seen clearly from Figs. 16 and 17 that the major features of the evolution of the collisions are captured very well by the DSMC method. The initial spike in collision frequency is due to the small bubbles getting displaced down which collide with larger bubbles on their way. Once the bubbles homogenize, it is followed by the successive high and low velocity core formations which mark the increase and decrease of collision frequency at around t ¼ 1 s and t ¼ 3:5 s. The DSMC algorithm over-predicts the total collisions in the initial regimes for both the cases. This is expected as the DSMC performs collisions with random neighbours. Since the relative velocity is low during the initial stages (system starts from rest), the collision direction is biased and not isotropic. This leads to small overlaps between the bubbles which are then treated by the explicit nearest neighbour collisions introduced in the algorithm. These overlaps lead to extra collisions. The collision frequency stabilizes to the same value as the DBM as the liquid momentum within the system increases. This happens  Table 5 applied to P2 for a polydisperse system, DBM vs DSMC.
because of increased bubble induced turbulence which corresponds to a higher bubble Stokes number. The evolution of the velocity profiles in case 1 is nearly identical for both the methods (see Fig. 18). The slight under-prediction by the DSMC in case 2 is for the same reason as explained in Section 4.1.2 for the 4 mm mono-disperse case.

Validation
Validation of the DSMC method applied to bubbly flows is performed by simulating a lab-scale bubble column. The simulation results are compared with experimental data from Deen et al. (2001). The column and simulation settings are specified in Table 6. The bubble size is assumed to be 4 mm (Deen et al., 2001;Darmana, 2006). A 3D schematic drawing of the column is shown in Fig. 19 with reference directions.
The boundary conditions implemented are the same as in Darmana et al. (2006). A no slip condition is implemented at x ¼ 0; x ¼ X MAX ; y ¼ 0; y ¼ Y MAX and z ¼ 0 for the liquid phase. At z ¼ Z MAX , a wall boundary with free slip is prescribed for the liquid phase. This is an approximation for the free surface. On each wall in x and y directions, one cell below the top wall, small windows of 3 cells each are prescribed with an inlet/outlet pressure prescribed boundary condition. Since the bubbles are allowed to enter and leave the column, the volume addition and loss has to be accounted for in the liquid phase.    Table 5 applied to P2, DBM vs DSMC.

Table 6
Simulation settings for the bubble column from Deen et al. (2001).

Parameter
Value Units  Fig. 20. The DSMC does not produce the exact time history in comparison to the deterministic DBM. This is expected since the exact collisions are not executed in the DSMC algorithm. Nevertheless the plume fluctuations are captured well by the DSMC collision model.
The measured parameters are time-averaged axial velocity, axial velocity fluctuations and the lateral (x-direction) velocity fluctuations (see . The z-velocity and z-velocity fluctuations match the experimental data better in the DSMC cases. DSMC smears out the exchange of momenta among bubbles by executing random collisions within the searching scope. This is also evident from the average velocity profiles shown in Figs. 10 and 18 where DSMC simulations under-predict the velocities. We believe that the smearing out of momentum exchange occurring in DSMC approximately mimics the smearing out of momentum due to sub-grid scale hydrodynamic interactions, causing nearby bubbles to equilibrate their velocities in the experiments. Of course the physical origin is different, but the effect is approximately the same, which we think is the reason for the better match. The overprediction in the center line z-velocity has been reported by Darmana et al. (2006) for their DBM. The x-velocity fluctuations are captured very similar to the DBM.

Single-core performance
The purpose of development of our model is to ultimately simulate dense bubbly flows in an industrial scale slurry bubble column for which the discrete bubble model is computationally too expensive. This section is dedicated to assess the computational performance of the algorithm in comparison with the DBM. The system used for the performance check is the 3D-periodic bubble rise discussed in the previous sections. The performance study is based on the time taken by the algorithm to execute a collision in a given time step. This is done to keep the comparison fair since  both algorithms may not produce exactly the same number of collisions for a given situation. The total execution time is therefore divided by the total number of collisions giving time consumed per collision (see Fig. 24, left). The total execution times per discrete phase time-step is also reported (Fig. 24, right).
It can be noted that for gas fractions less than 5 percent, the DSMC algorithm performs as well as the DBM model. As the gas fraction increases the DSMC becomes more efficient in handling the collisions. This can straight away be attributed to the overhead added by the expensive encounter list that needs to be searched and updated every time a collision occurs.
The parallelization potential of DSMC is much better than that of the DBM since the loops are individual particle based and not collision time based. Once the domain is decomposed, it can remain nearly independent of the remaining sections of the domain (except for the ghost cells). In the case of the DBM, the time frame in consideration for the movement of the particles needs to be communicated and this feature results in more communication overhead if we consider MPI parallelization.

Conclusion
A stochastic DSMC model has been developed to simulate dense bubbly flows. The algorithm is now capable of handling relatively high discrete phase volume fractions due to the addition of the radial distribution function gðrÞ. The model is then extended for bubble-bubble collisions in the presence of surrounding liquid via the addition of an explicit check for nearest neighbour collisions. This addition accounts for collisions in the environment of highly correlated bubble velocities and also alleviates the volume fraction limitation occurring due to overlapping bubbles in case of random collisions.
The model is first verified for a dry granular case and then verified in a 3D periodic bubbly flow system. This is done for both mono-disperse and poly-disperse systems. The verification model used here in both cases is the Discrete Particle/Bubble Model. It is then further validated with the experimental PIV data by simulating a lab-scale bubble column from Deen et al. (2001). Performance analysis shows that the DSMC shows speed-ups of up to 100 times at high gas fractions. The current study was limited to parcel size 1. Its extension to higher parcel sizes will be the subject of future research. DBM: total bubbles = 10000 DSMC: total bubbles = 10000 DBM: total bubbles = 25000 DSMC: total bubbles = 25000 DBM: total bubbles = 50000 DSMC: total bubbles = 50000