Study on bubble-induced turbulence in pipes and containers with Reynolds-stress models

Bubbly flow still represents a challenge for large-scale numerical simulation. Among many others, the understanding and modelling of bubble-induced turbulence (BIT) are far from being satisfactory even though continuous efforts have been made. In particular, the buoyancy of the bubbles generally introduces turbulence anisotropy in the flow, which cannot be captured by the standard eddy viscosity models with specific source terms representing BIT. Recently, on the basis of bubble-resolving direct numerical simulation data, a new Reynolds-stress model considering BIT was developed by Ma et al. (J Fluid Mech, 883: A9 (2020)) within the Euler—Euler framework. The objective of the present work is to assess this model and compare its performance with other standard Reynolds-stress models using a systematic test strategy. We select the experimental data in the BIT-dominated range and find that the new model leads to major improvements in the prediction of full Reynolds-stress components.


Introduction
Bubble-laden turbulent flows occur in a large number of processes in the energy generation, chemical industry, and nature (Lohse, 2018). In large-scale CFD simulations of these flows, the Euler-Euler (EE) method coupled with Reynolds-Averaged Navier-Stokes (RANS) modelling is the only viable approach. In this EE-RANS method, only continuous statistical quantities are calculated, and hence, all two-phase phenomena related to the gas-liquid phase boundaries need to be modelled (Ishii and Hibiki, 2006). In particular, one of the most challenging and complicated topics is turbulence modelling (Balachandar and Eaton, 2010) since the bubbles can strongly alter the carrier phase turbulence. One can readily cite several related mechanisms: (i) turbulence production due to the bubble wakes (Riboux et al., 2010); (ii) enhanced turbulence dissipation rate in the vicinity of the bubble surface (Ilić, 2006;; (iii) modulation of the liquid mean velocity profile due to interphase momentum transfer, resulting in a change of shear-induced turbulence (SIT) (Lu and Tryggvason, 2013;du Cluzeau et al., 2019;Bragg et al., 2021); (iv) increased turbulence intermittency in the dissipation range . These mechanisms can significantly affect the interfacial processes, such as bubble coalescence and breakup (Liao et al., 2014;Liao, 2020;Qi et al., 2020) and interphase momentum transfer (Lucas et al., 2016;Michaelides et al., 2016;Sommerfeld, 2017). Hence, providing an improved representation of some of these mechanisms is the goal of the turbulence modelling development for bubbly flows (Wörner and Erdogan, 2013;Lucas et al., 2015).
Over the last three decades, remarkable works have been accomplished in developing two-equation RANS models with source terms capturing the so-called bubble-induced turbulence (BIT) (Lopez de Bertodano et al., 1994;Morel, 1997;Troshko and Hassan, 2001;Rzehak and Krepper, 2013;Colombo and Fairweather, 2015;Ma, 2017;Vaidheeswaran and Hibiki, 2017;Liu et al., 2018;Liao et al., 2019;Magolan et al., 2019). These linear eddy-viscosity models take the bubble influence into account by adding source terms in the turbulence transport equations for both the turbulent kinetic energy (TKE, k) and a turbulence length-scale related parameter, e.g., the turbulence dissipation rate ε. As a result, this modulate as well the eddy-viscosity, by taking the expression 2 t / μ ν C k ε = for a k-ε type model for instance ( 0 09 μ C = . , typically). However, despite the fact that these Vol. 4, No. 2, 2022, 121-132 Experimental andComputational Multiphase Flow https://doi.org/10.1007/s42757-021-0128-0 Y. Liao, T. Ma 122 two-equation RANS models may be able to compute k and ε sufficiently, this type of model suffers from some issues concerning the standard eddy-viscosity concept: when applied to bubble-laden turbulent flows, as discussed in Ma et al. (2020aMa et al. ( , 2020b. Here, δ ij is the Kronecker delta, u u u ¢ ºis the fluctuating liquid velocity, and ⋅⋅⋅ denotes a Reynolds average (estimated using appropriate space and time average).
For most single-phase flows, the turbulence production highly relates to the mean velocity gradients doing work on the Reynolds stresses, i j u u ¢ ¢ , in the fluid, so that there is some rational for the eddy-viscosity concept. In BIT-dominated flows, however, the turbulence production mainly depends upon the interfacial energy transfer between liquid and bubbles (Kataoka and Serizawa, 1989). As a result, closures for i j u u ¢ ¢ based on the mean velocity gradient may not be meaningful for such a flow-at the basic level, the deficiency remains for the representation of i j u u ¢ ¢ in Eq. (1) that is intrinsic to this general hierarchy of EE k-ε type model, hence, cannot be overcome by more complex BIT terms in the corresponding turbulence transport equations. In view of these issues, several groups have sought to develop differential Reynolds-stress models (RSMs) based on the EE method that are more appropriate for BIT-dominated flows. We try to category their different Reynolds-stress closures into two classes. The first class refers to the works simply adapting the single-phase RSM in EE simulations of bubble-laden flows without considering the BIT effect. One of those examples can be taken in the study of Masood et al. (2014), who applied different single-phase RSMs in EE simulations of a bubble plume of Deen et al. (2001) and showed good agreement in the liquid single-point statistics with the experimental data. However, the success in this special case cannot be generalized to all kinds of bubbly flows. In contrast, the experiment of Deen et al. (2001) has also been simulated very accurately in terms of averaged liquid statistics by using different kinds of scale-resolving simulations without considering any BIT effect (Deen et al., 2001;Ničeno et al., 2008;Ma et al., 2015a;Ullrich et al., 2021). In these cited references, they show that compared to BIT the undulatory modulation of bubble plume is the dominant effect in such a flow, generating the most velocity fluctuations. Hence, considering the BIT effect or not influences the simulation only to a relatively small extent.
The BIT effect is considered in the second class of RSM applied to bubbly flows. Early attempts in this category, e.g., Lopez de Bertodano et al. (1990) combined the RSM of Launder et al. (1975) (LRR for short) with the BIT expression by Biesheuvel and van Wijngaarden (1984) as a source tensor to capture the turbulence generated by bubbles. This algebraic BIT expression is derived by using the assumption of a potential flow and considers the liquid velocity fluctuations influenced by the moving bubbles primarily from the displacement of the liquid. However, such an assumption cannot be applied for real situations of finite-size bubbles (say several millimetres), where the wake effect is dominant (Riboux et al., 2010). Further EE-RSMs for bubbly flows were performed by Cokljat et al. (2006), Colombo and Fairweather (2015), and recently by Ullrich et al. (2021). While Cokljat et al. (2006) and Ullrich et al. (2021) employed an isotropic source tensor to capture BIT, Colombo and Fairweather (2015) used a larger BIT source term in the direction of slip velocity between bubble and fluid. In all the above studies, the additional BIT tensor in the RSM resulted from various phenomenological modelling arguments.
Recently, Ma et al. (2020b) developed an RSM based on a term-by-term modelling strategy, having the interface resolving direct numerical simulation (DNS) data Fröhlich, 2015, 2016) as a reference for the particular term in the modelled Reynolds-stress transport equation. The focus of the present work is on the application of this model in different BIT dominant flows to assess its prediction on all the Reynolds-stress components. After introducing some technical details of this RSM in Section 2, a brief summary of the validation cases and the simulation setup are provided in Section 3. The results obtained by the RSM and its comparison with the experimental data are described in Section 4. A conclusion on the model performance is given in Section 5.

Reynolds-stress model with BIT source tensor
In this section, the complete formulation of the Reynoldsstress model of Ma et al. (2020b) is presented. First, the transport equations for the Reynolds stresses are given: Study on bubble-induced turbulence in pipes and containers with Reynolds-stress models 123 where ,

R ij
S and k S denote the BIT sources for Reynolds stresses and turbulent kinetic energy, respectively. Here, k S is adopted from Ma et al. (2017), which has been validated for a number of bubbly flow cases (Liao et al., 2019).
The modelled pressure-strain term ij  in Eq. (2) is identical to the SSG model (Speziale et al., 1991), and its two-phase version is given in Appendix A1. α l is the liquid volume fraction, v l is the liquid molecular kinematic viscosity, and s 1 63 c = . is an empirical constant. Moreover, the bubble Reynolds number p p r / is based on bubble diameter d p , u r is the mean relative velocity between the liquid and bubble (with r r u   º u ), and the drag force The rate of turbulence dissipation ε in Eq. (2) can be obtained from its own transport equation: with all single-phase related terms having the standard form (Hanjalic and Launder, 2011), and the empirical constants are 1 1 45 ε C = . , 1 36 ε σ = . , and 2 1 83 ε C = . . The interfacial term for the ε-equation, ε S , was modelled in Ma et al. (2017) and is adopted here without modification, i.e., / . In this expression, k S is given by Eq. (4), C D is the drag coefficient, and is the time scale characterizing BIT. The major advancement in this model is that not only the source terms k S and ε S , but also the expression of ij b * in Eq.
(3) are derived on the basis of the DNS data from Fröhlich (2015, 2016)

Investigated cases and numerical settings
The Reynolds-stress model presented above is applied in simulating the vertical pipe flows of Hosokawa and Tomiyama (2013), where the unladen bulk flow is laminar in the experiment and a homogeneous bubble column case of Akbar et al. (2012). These cases are chosen in order to minimize the effect of SIT as well as its interaction with BIT, which is desirable for a pure validation of BIT closures. Furthermore, the dataset provides full Reynolds-stress components of the liquid phase, which is rare in the published experimental data. Besides the comprehensive turbulence information in the carrier phase, the dispersed phase in the presently considered cases exhibits constant bubble diameters, leading to another benefit for validating BIT model, namely, we can focus on checking the RSM without being influenced by the uncertainty caused by modelling the polydispersity in the momentum and turbulence equations.
The experimental conditions of the three pipe flow cases are briefly summarized in Table 1. The bulk Reynolds number Re b = 900 is defined according to the bulk liquid velocity J l and pipe diameter D, which is the same for all the cases considered here. Water and air were fed into a vertical pipe of 2 m in height at room temperature and atmospheric pressure (see Fig. 1). The gravitational force acts in the negative x-direction. Gas volume fraction, mean liquid/gas velocity, and Reynolds-stresses were measured radially at an axial distance of 1.7 m from the bottom, where the flow is regarded fully developed. The complete experimental setup and the applied measurement techniques are described in Hosokawa and Tomiyama (2013), and we, therefore, refer the reader to that paper for more details.
The homogeneous bubble column case of Akbar et al. (2012) is used as another validation case. This case features a rectangular water/air bubble column, with a gas superficial velocity of 3 mm/s. Its height (x), width (y), and depth (z) are 800, 240, and 72 mm, respectively. Here, the coordinate system of the previous cases is used for improved readability, which is different from the one employed in the original paper. The global gas void fraction is 1.285%, and the bubble  (Akbar et al., 2012). The data for comparison are taken along the measurement line from the wall to the centre (half the column width) at a height of 500 mm in the midplane (z = 36 mm), see Fig. 2. More details are provided in the cited reference.
The simulations in the present study are performed with the ANSYS CFX 19.5 software. The computational domain and setup are similar to those adopted in our previous work , where the same pipe flow cases were investigated using OpenFOAM 6 with an SST k-ε turbulence model including the BIT term of Ma et al. (2017). Interphase momentum transfer resulting from both drag and non-drag forces, i.e., shear and wall lift, turbulence dispersion, and virtual mass, is taken into account. The corresponding interfacial transfer models employed here are summarized in Table 2. A complete description and further validation studies about all interfacial transfer models can be found in the so-called baseline model (Liao et al., , 2018Rzehak et al., 2017;Ziegenhein et al. 2017;Colombo et al., 2021) It is important to stress that our previous study , however, could not well reproduce the radial void fraction distribution in the pipe flow cases (shown in Section 4). The reason is related to the Tomiyama lift model,  which would cause a wall-peak or intermediate profile of void fraction based on the relatively small bubble sizes in the considered cases. In contrast, core-peak profiles were measured in the experiment for these three cases. Recall that the major objective of the current work is to better understand BIT and validate the Reynolds-stress closure model. If, however, the simulated void fraction does not fit the experimental value, the BIT model cannot be assessed. Hence, a prerequisite is that the uncertainty created by the other submodels is removed as much as possible.
To accomplish this, adjustment is performed on the lift coefficient C L in order to obtain a satisfactory agreement between the simulated and measured void fraction profiles. More specifically, an iterative process was carried out by running EE-RSM simulations simultaneously optimizing the lift coefficient C L , while the other interfacial force models were employed as they are in Table 2. The adjusted values of C L for each case are given in Appendix A2.
Considering that the pipe flow is uniform and axisymmetric, a two-degree wedge instead of the whole pipe is simulated, and a quasi-two-dimensional grid with one layer of cells in the circumferential (z) direction is applied. The mesh resolution in the radial and axial directions is presented in Table 3. The maximum cell size is determined based on the study of Liao et al. (2020), which showed a resolution finer than 0.5 mm and 10 mm in the radial and axial directions, respectively, is sufficient to provide mesh-independent results (see a study of grid refinement in Appendix A3). For the bubble column case, the full water region (x < 700 mm) is considered with a 175 × 240 × 36 grid in the x, y, and z directions, respectively.
Towards boundary conditions, velocity inlet with measured gas/liquid mass flow rates, pressure outlet for pipe flow, and degassing for bubble column is employed. On the walls of all the considered cases, non-slip and free-slip conditions are assumed for the liquid and gas phases, respectively, and the scalable wall function is selected for modelling flow near the wall. A maximum domain imbalance of mass and momentum of 0.1% and residual of 10 −5 are set as convergence criteria.

Vertical pipe flow
The calculated radial profiles of void fraction, liquid/gas velocity, turbulent kinetic energy, and Reynolds stresses using the present RSM have now been validated against the measurements from Hosokawa and Tomiyama (2013). For comparison reason, additional simulations are carried out using other standard models from the literature under exactly the same conditions, namely employing the same overall optimisation procedure described in Section 3. These models are from Cokljat et al. (2006) and Colombo and Fairweather (2015). The former proposed an isotropic BIT source term, i.e., just the same as the present RSM. Figure 3 shows the distribution of gas volume fraction. As mentioned above, the profiles are optimized by adjusting the lift coefficient in order to put the focus only on the BIT issue. Consequently, the simulation results, which are illustrated with dotted (for isotropic BIT source), dashed (for "1:0.5:0.5"), and solid (for the present RSM) lines, coincide with each other and agree well with the experimental data in all the three cases. As additional information, the results of our previous work  employed with the lift coefficient of Tomiyama et al. (2002) are included (see the dash-dot lines). One can notice that Cases H1, H2, and H3 all exhibit a core-peak profile in the experiment, and the peak increases with the average void fraction and bubble diameter (see Table 1). Based on the lift correlation of Tomiyama et al. (2002), all the cases have a positive lift coefficient C L = 0.288 for the considered bubble parameters, resulting a flat void fraction profile in our previous work. Figure 4 presents the vertical component of gas velocity. As one can see, slight over-prediction exists in all three cases. Nevertheless, the profiles are in consistent with those of void fractions, and the general slope in Cases H1-H3 becomes gradually steeper. Note that the deviation in the near-wall region can be ignored, since the void fraction there is nearly zero. The present three sets of simulations using different values of 11 b * , 22 b * , 33 b * conform well with each other. Figure 5 shows the vertical component of liquid velocity. The results of all three cases are in good agreement with the experimental data. The peak of the parabolic profiles increases from Cases H1 to H3, which corresponds to the void fraction distribution presented in Fig. 3. Even though a small discrepancy exists in the pipe center region in Cases H1, H2, and H3, the three simulations with variable distribution of BIT sources are in accordance with each other. It is important to emphasize that the consistency in the prediction of gas volume fraction and velocity profiles shown in Figs. 3-5 is the prerequisite for a fair comparison with the different ij b * approaches in RSM.  We now turn to look at the turbulence parameters. As expected, the TKE provided by the three simulations is identical (see Fig. 6), since we keep the source terms for k in Eq. (4) and ε (involved in Eq. (5)) exactly the same . In all the core-peak cases, the simulation results agree satisfactorily with the experimental data. It shows that the proposed model for bubble-induced production and destruction of TKE in the liquid phase is reliable.
As aforementioned, the present three sets of simulations only differ in partitioning k source into three components of Reynolds normal stresses, i.e., equally, or with fractions of either 1:0.5:0.5 or ij b * suggested by Eq.
(3) in the present work. In the experiment, u u ¢ ¢ was measured approximately twice as large as the other two directions in all these cases over a large range of radial direction (Hosokawa and Tomiyama, 2013), since it is in the slip velocity direction, which contributes to the main BIT source. This component is well reproduced in Fig. 7 by the present work, having the highest value of u u ¢ ¢ . The bubble Reynolds numbers in these three pipe cases are very similar ~900, resulting 11 1 4 b * » . in Eq. (3), which is higher than the values based on the "isotropic" and the "1:0.5:0.5" approaches. As expected, the "isotropic" assumption underpredicts u u ¢ ¢ substantially, and the result obtained by "1:0.5:0.5" is in the middle. Additionally, we show in Appendix A4 a simulation without any BIT sources ( , 0 R ij S = and 0 ε S = ). The results reveal clearly the importance of adding the BIT terms for calculating Reynolds stresses.
To illustrate the turbulence contributed by different sources in Eq. (2) for the streamwise direction, we plot in Fig. 8 the shear-induced production term P 11 and the interfacial term 11 R S , . Although they are terms based on the   Study on bubble-induced turbulence in pipes and containers with Reynolds-stress models 127 RSM modelling equation, it indicates an interesting issue in the present flow configuration. The unladen background flow is indeed laminar in the experiment (Hosokawa and Tomiyama, 2013). By adding bubbles, the flow, however, becomes turbulent and initializes the same time as well the shear-induced contribution if there is a mean flow gradient (in the case for all the considered cases here). The shearinduced contribution compared to the interfacial term is generally small in the present BIT-dominated flows over a wide range across the pipe, but the contribution is not anywhere negligible (e.g., close to the wall where the shear rate is relatively high). Figures 9 and 10, depicting Reynolds normal stresses in the y and z directions, are measured approximately equal in the experiment. This is due to the fact that in the BITdominated cases, the turbulence is generated by the bubble wakes, which individually can be approximated by the axisymmetric turbulence framework to some extent (George and Hussein, 1991;Ma et al., 2020b). Again, the present work gives generally superior results to the other two approaches, which obviously predict significantly higher levels of both components.
We plot in Fig. 11     work  using the Boussinesq hypothesis in the framework of a two-equation eddy-viscosity turbulence model.

A bubble column case
While in the previous subsection we have tested the RSMs for bubbly flow in a pipe with varying gas void fraction but a similar bubble Reynolds number, we consider in the following a bubble column case (Akbar et al., 2012) with different bubble sizes and without bulk velocity.
Three EE-RSM simulations based on different BIT sources were performed for this case. Based on own previous work, the same set of interfacial force models as detailed by Liao et al. (2019) (Table 2) were employed to achieve good agreement of gas void fraction and vertical liquid/gas velocity at all simulations (see Fig. 12). The flatter profiles for the liquid and gas velocities may be due to the slight inaccuracy in the gas void fraction close to the wall. Concerning the liquid Reynolds stresses, as shown in Fig. 13, the present BIT model provides again the most accurate results among these three BIT models, as well as our own previous work using large eddy simulation (Ma et al., 2015b). The ratio of the Reynolds normal stresses 2 5 1 u u v v ¢ ¢ ¢ ¢ : » . : observed in the experiment is very well reproduced for the centre region of the measurement line.

Conclusions
In the present work, a full Reynolds-stress model recently proposed by Ma et al. (2020b) based on the budget analysis of the Reynolds stress and its dispersion from DNS for   bubbly flows is tested rigorously for a large number of bubbly flows. The results of other two standard models (labelled as "isotropic" and "1:0.5:0.5", respectively) are included in the comparison. In all EE-RSM simulations, the lift force coefficient has been carefully fine tuned, so that the simulated gas void fraction fits very well to its experimental distributions. This procedure removes the uncertainty created by the other submodels as much as possible and provides an optimal condition for the comparison of different BIT turbulence models. As expected, the model proposed by Ma et al. (2020b) leads to a substantial improvement for all the considered cases with or without bulk flow involving different void fractions and bubble Reynolds numbers and should be the RSM of choice for dispersed bubbly flow applications. The present BIT model is now tested for steady fully developed bubbly flows with very encouraging results. In the future work, we will implement this BIT model within the scale-resolving simulations framework (Fröhlich and von Terzi, 2008) as suggested by Ullrich et al. (2021), so that the Euler-Euler approach coupled with steady or unsteady RANS modelling can accurately predict the turbulence parameters both for the steady full developed cases, as well as unsteady cases, like bubble plume. = . = . = .

A2 Adjusted lift coefficient
The adjusted lift coefficients for the considered cases are listed in Table A1.

A3 Grid refinement study for Case H3
We present the results for Case H3 using four different meshes in Fig. A1. All results displayed in the main text are from the second mesh with the grid numbers 41 and 801 in radial and axial directions, respectively.

A4 Case H3 without BIT terms
In Fig. A2, we present the results of Case H3 without considering BIT terms ( , 0 R ij S = and 0 ε S = ).

Fig. A1
Comparison of the gas void fraction, mean liquid velocity, and u u ¢ ¢ obtained by four different meshes.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.