Comparison of Eulerian QBMM and classical Eulerian–Eulerian method for the simulation of polydisperse bubbly flows

Funding information Helmholtz-Zentrum Dresden-Rossendorf Open access funding enabled and organized by Projekt DEAL [Correction updated on 31 October 2020 after first online publication.] Abstract The spatial gas distribution of poly-disperse bubbly flows depends greatly on the bubble size. To reflect the resulting polycelerity, more than two momentum balance equations (typically for the gas and liquid phases) have to be considered, as done in the multifluid approach. The inhomogeneous multiple-size group model follows this approach, also combined with a population balance model. As an alternative, in a previous work, an Eulerian quadrature-based moments method (E-QBMM) was implemented in OpenFOAM; however, only the drag force was included. In this work, different nondrag forces (lift, wall lubrication, and turbulent dispersion) are added to enable more complex test cases to be simulated. Simulation results obtained using EQBMM are compared with the classical E–E method and validated against experimental data for different test cases. The results show that there is good agreement between E-QBMM and E–E methods for mono-disperse cases, but E-QBMM can better simulate the separation and segregation of small and large bubbles.

motion. It may change its sign depending on the bubble size, as the lift force experienced by small spherical bubbles differs from that experienced by large ellipsoidal ones. In upward bubble flows, the lift force pushes the small bubbles toward the wall and the large bubbles toward the center of the pipe. The wall lubrication force prevents bubbles from collecting at the wall. Moreover, the turbulent dispersion force considers the bubble dispersion due to the effect of turbulent eddies. It is proportional to the gas volume fraction gradient and flattens the corresponding profiles. It has been shown to improve the stability of the E-E method. 4 As shown in Figure 1, all these forces have an effect on the bubbles' migration, and as a result the profile of the phase fraction distribution is established as the outcome of multiple factors. In particular, when only the drag force is included, all the bubbles move upward with negligible radial migration. The predicted phase fraction is flattened somewhat, since no lateral forces are exerted. If the lift force is included, the small bubbles tend to move toward the wall forming a gas volume fraction profile with a maximum directly at the wall. If the wall lubrication force is also included, the predicted phase fraction profile exhibits a peak next to the wall, since the wall lubrication force pushes the bubbles away from the wall. When in addition the turbulence dispersion model is considered, the phase fraction profile is flattened, since it acts as a diffusion term. It should be noted here that the statement above is quite general. For example, whether the wall lubrication model developed by Hosokawa et al 5 or Antal et al 6 is used, the predicted phase fraction exhibits similar resulting features with a wall peak. Besides the phase fraction distribution, the gas holdup is also important. However, the consideration of the column aspect ratio and its influence on the holdup concerns global parameter. The investigation of critical values is not the objective of this work. Readers are referred to other works for more information. [7][8][9][10] Thanks to increasing computational power, computational fluid dynamics (CFD) methods have become more and more feasible as a means of examining such complex flows. The macroscopic Eulerian-Eulerian (E-E) method has been employed extensively to investigate bubbly flows in many studies. [11][12][13][14][15][16][17] However, one drawback of the classical two-phase E-E method is that it can only be employed for monodisperse systems. Here, mono-dispersity refers to situations in which all the bubbles have the same properties (e.g., bubbles with identical diameters). In contrast, poly-dispersity refers to situations in which the properties of the disperse-phase entities are different for each entity (e.g., bubbles with different diameters). In practice, mono-disperse bubbly flows are relatively rare, hence, it will be important to have a modeling framework that naturally accounts for poly-dispersity. In order to extend the method to include poly-disperse systems and to allow for the consideration of bubble coalescence and breakup, CFD was coupled with population balance models and employed to simulate poly-disperse systems. [18][19][20][21][22][23][24] In the CFD-PBM coupling procedure, the mean Sauter diameter is calculated by the PBM and fed into the momentum interfacial exchange terms. Table 1 shows the different PBM solving methods.
Readers interested in this subject are referred to another recent work. 25 These methods are valid only for small Stokes flows (e.g., liquid-liquid emulsions with relatively homogeneous droplet diameters) since the disperse phase is still convected by one velocity field for all bubble or droplet sizes. 26 However, many experimental studies on polydisperse bubbly flows have shown that, when polydispersity is significant, bubbles tend to migrate at different velocities depending on their size. These effects are caused by the fact that interfacial forces depend on the bubble size; this is particularly important for the lift force, where the inversion of sign of the lateral lift force in vertical upward bubbly flows even causes small and large bubbles move toward opposite directions. [27][28][29] To predict such polydisperse multiphase flows, CFD-PBM coupling was extended to include other multiphase E-E methods, in which bubbles within a range of sizes are handled as separate phases (the inhomogeneous multiple-size group [MUSIG] model). [34][35][36] Another solution is to systematically couple the CFD with a more gen- with "Eulerian" referring to the approach for the continuous phase and "QBMM" denoting the quadrature-based moments method. In fact, the transport equations for the moments of the NDF are closed using a quadrature approximation. To model the phase velocity of the disperse phase, different methods can be chosen, such as the conditional quadrature method of moments (CQMOM) 37 or the velocity polynomial approximation (VPA) model. [38][39][40] In CQMOM, the disperse phase velocity is treated as a separate "internal coordinate," and the quadrature approximation is used to overcome the closure problem.
The quadrature approximation consists of different abscissas or nodes that can be thought of as separate bubble classes or "phases." This approach was mainly used to predict particle trajectory crossing. In the VPA model, by contrast, the disperse phase velocity is not formally treated as an independent internal coordinate, and a polynomial relationship is assumed between the bubble velocity and bubble size.
In this work, we implemented the full set of interfacial momentum where α c is the phase fraction of the continuous phase, ρ c is its density, and U c is its average velocity, R c represents the stress tensor, p is the average pressure, g is the gravity acceleration vector, and A is the momentum interface exchange term, which will be discussed in the next section.

| GPBE/QBMM for the disperse phase
Considering the NDF n(t, x, d, V d ), the GPBE for the disperse phase can be written as 38 where V d is the bubble velocities, S is the possible source term, which is neglected in this work, A is the continuous rate of change of particle where d is the particle size and P i are the velocity polynomial coefficients (VPCs). Equation 4 can be written in the following where v x , v y , v z are the velocity components in the Cartesian coordinate system and p i,x , p i,y , and p i,z are the components of P i . Once the VPCs are determined, the conditional particle velocities V(d) can be updated from the particle size.
To solve Equation 3 numerically with the VPA, the operator splitting procedure is employed in this work. 38 In the first step of operator splitting, only convection is considered: This equation can be transformed into the mixed-moments transport equations by defining of the mixed moments of the NDF: When m j,k,l,i in Equation 7 are solved, m j + 1,k,l,i , m j,k + 1,l,i , and m j,k,l + 1,i need to be known, which implies that Equation 7 is not closed. In this work, it is closed using a quadrature approximation, as formulated in the original quadrature method of moments, 30 in which any higher order mixed moment can be calculated as where w β and d β are the weights and abscissas of the quadrature approximation calculated from the pure moments of the NDF with respect to bubble size, N is the number of nodes, v are the three components of the bubble velocity and an assumption is made on the dependence of bubble velocity versus bubble size. The N weights and abscissas can be calculated from the first 2 N moments using the moment inversion algorithm (e.g., the Wheeler algorithm ). 43 The calculation is done using only the first 2 N moments of the NDF with respect to bubble size m 0,0,0,i with i = 0,1,…2 N-1. Moreover, the velocity polynomial coefficients P i with i = 0, 1, 2 as defined in Equa- where α β is the phase fraction of the disperse phase with abscissas (particle size) d β , A β, buo denotes the buoyancy force, A β, drag denotes the drag force, A β, lift denotes the lift force, A β, wall denotes the wall lubrication force, and A β, bubblePres denotes the bubble pressure force.
The buoyancy force can be modelled as in which g is the gravity acceleration vector and ρ β is the density of the disperse phase with size d β . For the bubbly pipe flows investigated in this work, the buoyant force acting on the radial direction accelerates the bubbles' upward movement. The drag force can instead be modeled as In vertical upward bubbly flows, the drag force counterbalances the buoyancy force and this determines bubble's terminal velocity.
The lift force plays a critical role in the prediction of the lateral behavior of bubbly flows. Therefore, a correct description of the lift coefficient in bubble columns is crucial in order to model this transversal force correctly. It can be calculated by the following expression: where C L, β is the lift force coefficient, which can be calculated by different models. 27 The main feature of the lift force is that the model is capable of predicting the so-called cross-over point, at which bubble distortion causes a reversal in the sign of the lift force. The coefficient C L, β becomes negative for large bubbles (e.g., for air bubbles in pure water larger than 5.6 mm), largely affecting the dynamics of bubble radial and axial redistribution in horizontal pipes. It should be noted that the critical value at which the sign of C L, β changes may differ slightly. It depends mainly on the shape and dimension of the bubble.
Besagni and Inzoli 45 found that the bubble aspect ratio depends on the bubble size in bubbly flows. Spherical bubbles are relatively rare compared with ellipsoidal bubbles. Different models have been proposed to model the bubble aspect ratio. 45,46 Recently new measurements on the lift force for air-water flows were done, and the effect of the bubble shape was discussed. 47 The validity of the change of the sign of the lift force in dependence on the bubble size was also shown for ploy-disperse flows turbulent flows. 28,47 There might be some limits for highly turbulent flows. However, the proper correlation for the evaluation of the lift force is still controversial. Some authors claim that the intensity of turbulent fluctuations also plays an important role, but most researchers agree that bubbles of different sizes are subjected to a lift force acting in opposite directions.
Another aspect of the lift model is that the results obtained using Equation 13 show a gas bubble radial distribution that peaks at the wall, due to the high continuous phase velocity gradient. The predicted accumulation is abnormal and not reflected by experimental findings. In order to handle this problem, a lift force damping model is usually employed. In the lift damping model, the lift coefficient is multiplied by a limiter, λ β , which can be calculated using the following expression: where y w is the distance from the center of the cell (here, the model is formulated for a finite-volume discretization) to the nearest wall.
Once the distance is smaller than 1.5 times the bubble diameter, the limiter gradually decreases to zero.
Another way to handle this problem is to include the wall lubrication force, which tends to push the secondary phase away from the walls. In bubbly upward flows in a vertical pipe, this force results in the disperse phase concentrating in a region near, but not immediately adjacent, to the wall. It is calculated as follows: where n is the unit normal pointing away from the wall and C w is the wall lubrication coefficient, which can be calculated using different models. Readers should refer to the work of Hosokawa et al 5 and Antal et al 6 for more information.
The bubble pressure force acts as a driving force for bubbles to move from areas with higher phase fractions to areas of lower phase fraction. It arises due to the pressure variations in the continuous phase that are not resolved at the mesoscale. The bubble pressure force is defined by 48 where C bp is the bubble pressure coefficient.
Once the V t + Δt β have been updated, the mixed moments can be updated in terms of their definition as reported in Equation 8. Readers interested in details of the second step of the operator splitting procedure are invited to read our previous work. 40 At last, we finalize the discussion of the numerics of the E-QBMM and the E-E method by

| Coupling and numerical discretization
From the initial settings, the disperse phase fraction can be calculated (e.g., α d = π 6 m 0,0,0,3 ). Neglecting the shared pressure gradient force, lift force, wall lubrication force, bubble pressure force and the contribution of the disperse phase in the drag force, the semidiscretised form of the continuous phase momentum equation is where a P and a N are the matrix diagonal and nondiagonal coefficients, which are a function of U c ; U c, P, and U c, N are the unknown velocity of the continuous phase defined at and near the cell center. The solution of Equation 18 is the predicted continuous phase velocity, U * c,P , defined by The nondrag forces together with the contribution of the disperse phase in the drag force are defined as the momentum flux, ϕ c, force, f : where N dβ is the number of abscissas, S f is the surface-normal vector, and f denotes the variables defined at cell faces. Substituting these into Equation 1 leads to the semi-discretised form: Once the pressure is computed from Equation 21, the corrected continuous phase velocity can be computed by the flux reconstruction in the PISO procedure.

| Turbulence model
The Reynolds stress arises in the momentum equations as a result of the averaging process. Different turbulence models can be employed to calculate the Reynolds stress, such as the k -ε model, 1 k -ω model, 49 or the k -ω SST model. 50 It was shown that the k -ω turbulence model yielded a better qualitative prediction of the bubble plume than the k -ε model, due to the low Reynolds number treatment of the former model. 49 Some other works show that the k -ε model can still yield good results in bubbly flows. 51,52 In this study, a two-phase k -ε model was employed. The equations are omitted for brevity and readers are refereed to other works for more details. 53 It should be noted that in general the bubble-induced turbulence (BIT) plays an important role. The presence of bubbles modifies the structure of the liquid turbulence field and the production of shear-induced turbulence, which in turns modifies the bubble distribution and the break-up and coalescence processes. These bubbles act as a source of the BIT, also generating turbulence in flows that would otherwise be laminar. In general, the BIT model includes a source term in the turbulence transport equations to account for the turbulence generated by the bubbles, and different models have been developed. [54][55][56][57] However, because BIT may not play a major role in the flows considered here, it neglected in this work.

| Numerical configurations
The momentum flux instead of source terms. 58 The higher order fluxsplitting realizable scheme is implemented for the moment transport equations 44 with a MIN-MOD limiter. 59 The flux-corrected transport (FCT) scheme with a multidimensional universal limiter with explicit solution (MULES) is employed for the phase fraction equation to ensure the boundedness of the phase fraction. 60 The coupling of the shared pressure and the continuous phase velocity is solved by the PISO procedure. 61 The adiabatic solver reacting TwoPhaseEulerFoam without mass-transfer was also employed as a baseline solver in which the E-E method is implemented.
All the grids in the following test cases are generated by blockMesh. The mesh size is selected according to a grid indepen- was shown in our preliminary investigations that such meshes achieve mesh-independent solutions. Mesh with higher resolution did not improve the predictions. The initial moments at the inlet are calculated by assuming a log-normal distribution. 62 The discretization scheme, sparse matrix solver, momentum closure models, and other details of these different test cases are listed in Tables 3-5. Since we employ the PISO algorithm, the pressure is solved iteratively at each time step. The relevant tolerance and the final tolerance equal 0.01 and 1 × 10 -7 , respectively. Other variables (e.g., k and ε) are solved after the pressure and velocity iteration procedure, and the tolerance equals 1 × 10 -7 . Since the moment transported equations are solved using the explicitly realizable scheme, a diagonal solver can be constructed and a relatively low tolerance (1 × 10 -15 ) is used. In order to minimize the time discretization error and to ensure the moments realizability, the time step is adjusted to ensure that the Courant number is smaller than 0.05. Unless otherwise stated, all the other simulations in the following sections are performed with identical settings.
It should be noted that the bubble diameter, phase fraction, and the velocity at the inlet work critically. Generally, the inlet disperse phase fraction and the bubble velocity are unknown except for the test case studied by Banowski et al. 29 The inlet velocity can be calculated from the superficial velocity using the following expression: In this work, we employ another approach for all the test cases due to its simplicity. In this method, the inlet phase fraction is fixed at 0.5, and then the inlet bubble velocity satisfying the superficial velocity condition is calculated. It was found in our preliminary investigation that there is no difference between these methods.   Similar results can be also found in our previous work in which the one-way coupled E-QBMM was employed to simulate particle-size segregation. 40 The diffusive feature of the first-order spatial scheme can smash the wall peak phase fraction profile into vertical upward bubbly flows. The continuous phase variables (e.g., the liquid velocity  and the turbulent kinetic energy) predicted by the first-order and higher order two-way coupled E-QBMM show no substantial difference, which implies that the spatial scheme of the disperse phase has little effect on the continuous phase. Therefore, in the following sections, all the test cases predicted by the E-QBMM were launched by the higher order spatial schemes.

| Test case studied by Žun
Since the bubble size was assumed to be 4 mm, the wall peak was successfully predicted by the E-QBMM and the E-E method due to the lateral lift force and wall lubrication force, as shown in Figures 2   and 3. Moreover, the phase fraction in the vicinity of the wall decreases sharply because the lateral wall lubrication force pushes the bubbles away from it. As the bubble pressure force and the turbulent dispersion force are not included in the E-QBMM and the E-E method, the wall peaks predicted by both methods are rather strong, but still consistent with one other. The predictions of the turbulent kinetic energy and the turbulent energy dissipation rate also show the typical wall peak trend in bubbly flows, which is consistent with other works. 28,65 To further validate the algorithm, three points were selected to monitor the local liquid velocities predicted by the E-QBMM and the E-E method. It can be seen in Table 6 that the local liquid velocities predicted by these two different methods agree well with each other, which further implies that the implementation is correct.
To validate the turbulent dispersion force and the bubble pressure gradient force, the simulation results predicted by the E-QBMM and the E-E method were compared against measured data. The constant and uniform bubble size (abscissas value) was assumed to be 4.1 mm.
As the author did not report the superficial gas velocity, it was assumed to be 0.5 m/s after a fitting procedure. A similar procedure can also be found in other works. 65 The turbulent dispersion force model proposed by reference 66, and the coefficient of the bubble pressure gradient force were assumed to be 2.0. 39     Now it is interesting to employ the E-QBMM to simulate a polydisperse system by adjusting the bubble sizes. Figure 10 represents the phase fraction predicted by the E-QBMM for bubble sizes of 4.1/4.9/5.8 mm and 4.3/4.9/5.6 mm. It can be seen that bubbles of size smaller than 5.5 mm tend to move toward the wall, but the maximum value of the phase fraction is located at different positions. The smaller the bubbles are, the closer they accumulate to the wall. On the other hand, the large bubbles (i.e., larger than 5.5 mm) tend to move toward the pipe center, which is consistent with the physics.
Moreover, the overall phase fraction profile predicted by E-QBMM still represents a wall peak even in this polydisperse system since the mean Sauter diameter is smaller than 5.5 mm.

| Test case by Banowski et al
Finally, we complete the study with a relatively recent work by Banowski et al. 29 In that work, air/water two-phase co-current upward flows in a vertical pipe are investigated using X-ray tomography. The test section comprises a vertical titanium pipe with an inner diameter of 54.8 mm and a length of 6 m. Gas is injected into the The numerical configurations in the E-QBMM simulations are identical with previous test cases. The phase fraction predicted by the E-QBMM and the E-E method are reported in Figure 11. Again, it can be seen that the double peak was predicted by the E-QBMM due to large bubbles moving toward the pipe center and small bubbles moving toward the wall. The two-phase E-E method fails to predict the double peak since it is only applicable for monodisperse phase systems. A multiphase E-E method combined with the inhomogeneous MUSIG model should be able to handle the problem, in which multiple momentum equations are employed for the disperse phase. 34 However, in order to make a justified comparison with the E-QBMM, in which only one GPBE is employed for the disperse phase, only the two-phase E-E method is investigated in this work.
The contour plots of the phase fraction predicted by the E-QBMM and the E-E method at different horizontal sections are shown in Figure 12. It can be seen that the plots predicted by the E-QBMM also capture the double peak in vertical upward poly-disperse bubbly flows in upper sections. It should be noted that in our previous work, 40 we found that the flow field information (e.g., the phase fraction distribution) is highly dependent on the inlet conditions for a sudden-enlargement gas-liquid test case with a relatively low L/R value. However, even if a uniform phase fraction is given for the test

| CONCLUSIONS
In this work, the higher order fully coupled E-QBMM with a full set of momentum interfacial exchange terms (e.g., the drag/lift/wall lubrication/bubble pressure forces) was implemented in open-source CFD code OpenFOAM-5.x, in which the conditioned disperse velocity is modeled by the VPA. The solver is called twoWayGPBEFoam, and it was employed to simulate monodisperse and polydisperse bubbly flows in vertical upward channels. reactingTwoPhaseEulerFoam was also employed as a benchmark model comparison, in which the E-E method is implemented. Different experimental data from different works featuring wall peaks and double peaks were selected to verify the algorithm and the implementations.
We show that the results predicted by the higher order E-QBMM are identical to those predicted by the E-E method. For bubbly flows, when the bubbles are small, both methods can predict the wall peak.
For bubbly flows with large diameters, the predictions using both methods show the bubbles moving toward the channel/pipe center due to the negative lift force coefficient, which is consistent with experimental data. When large bubbles and small bubbles exist alongside one another, constituting a polydisperse system, the double peak can be successfully predicted by the E-QBMM, since the bubbles of different sizes are transported at different velocities (see Equation 4). Moreover, the bubble pressure force and turbulent dispersion force can smooth the lateral phase fraction distribution, and the overprediction of the wall lubrication force incurs an underestimation of the phase fraction in a small region in the vicinity of the wall. Improvements can be achieved by a more sophisticated combination of the nondrag forces and the breakage and coalescence model.

ACKNOWLEDGMENT
This work was partially supported by Helmholtz-Zentrum Dresden-Rossendorf. One author (Dongyue Li) wants to acknowledge the continued encouragement of the CFD-China community. Open access funding enabled and organized by Projekt DEAL.