1 Introduction
Disperse turbulent bubbly flows occur in a very large number of situations encountered in process engineering, energy technology, environmental flows, etc. Examples are bubble column reactors, abundantly used in the chemical industry, pipelines, waste-water treatment, gas release at the bottom of the sea, and many others. Such flows can be investigated numerically by various approaches at different levels of detail. For large-scale simulations, the Euler–Euler (EE) approach (Ishii & Hibiki Reference Ishii and Hibiki2006) coupled with steady or unsteady Reynolds-averaged Navier–Stokes (RANS) modelling is the only viable framework. In this case, only continuous statistical quantities are computed, so that beyond the closures for single-phase terms all two-phase phenomena related to the phase boundaries need to be modelled. Extensive reviews of this approach were compiled by Joshi & Nandakumar (Reference Joshi and Nandakumar2015), Michaelides, Crowe & Schwarzkopf (Reference Michaelides, Crowe and Schwarzkopf2016) and Sundaresan, Ozel & Kolehmainen (Reference Sundaresan, Ozel and Kolehmainen2018). Among the various phenomena, the effect that fluid turbulence is generated by bubbles moving relative to the fluid is one of the most important and most delicate to model (Balachandar & Eaton Reference Balachandar and Eaton2010; Elghobashi Reference Elghobashi2019). Providing an improved representation of this mechanism, the so-called bubble-induced turbulence (BIT), is the goal of the present paper.
Over the last two decades, there has been widespread work on supplementing single-phase two-equation eddy-viscosity models with specific source terms representing BIT (Lopez de Bertodano, Lahey & Jones Reference Lopez de Bertodano, Lahey and Jones1994; Morel Reference Morel1997; Troshko & Hassan Reference Troshko and Hassan2001; Colombo & Fairweather Reference Colombo and Fairweather2015; Ma et al. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017). These models take the influence of bubbles into account by including additional source terms in the balance equations for both $k$, the turbulent kinetic energy (TKE), and $\unicode[STIX]{x1D700}$, the turbulent dissipation rate, or another equivalent parameter. This alters the turbulence quantities and, as a result, the effective transport coefficients, such as the eddy viscosity. Albeit widely used, the approach suffers from substantial uncertainties concerning the expression of the eddy viscosity in bubbly flows (Kataoka & Serizawa Reference Kataoka and Serizawa1989; Troshko & Hassan Reference Troshko and Hassan2001). These studies, as well as own experience (Ma Reference Ma2017), show that EE two-equation models should be applied with care for bubbly flows, as discussed in § 2.2 below.
Good alternatives to eddy-viscosity models are differential second-moment closures (SMCs). They employ modelled differential transport equations for all Reynolds-stress components required to close the RANS equations. For single-phase flow, SMCs incorporate more physics than the traditional two-equation models, resulting in a larger range of applicability (Hanjalić & Launder Reference Hanjalić and Launder2011), although the debate on when precisely the additional effort and numerical complexity are justified for single-phase flow, compared to other approaches, is not ultimately settled in the community.
For bubbly flows the situation is somewhat different since the buoyancy of the bubbles generally introduces anisotropy on the small scales. Hence, in many applications, there is considerable interest in accounting for anisotropic velocity fluctuations, which can be found in both experiment (Riboux, Risso & Legendre Reference Riboux, Risso and Legendre2010; Akbar et al. Reference Akbar, Hayashi, Hosokawa and Tomiyama2012; Lai & Socolofsky Reference Lai and Socolofsky2019) and direct numerical simulation (DNS) (Lu & Tryggvason Reference Lu and Tryggvason2008; Bolotnov et al. Reference Bolotnov, Jansen, Drew, Oberai, Lahey and Podowski2011; Santarelli & Fröhlich Reference Santarelli and Fröhlich2015), to name but a few. SMCs employed in the EE approach, however, have been investigated only by comparatively few groups. Lopez de Bertodano et al. (Reference Lopez de Bertodano, Lee, Lahey and Drew1990) were the first to address this issue in bubbly flow. They adopted the SMC of Launder, Reece & Rodi (Reference Launder, Reece and Rodi1975) – hereafter referred to as LRR – as the carrier model and combined it with the algebraic BIT expression of Biesheuvel & van Wijngaarden (Reference Biesheuvel and van Wijngaarden1984) as a source tensor to represent the additional energy produced by bubbles. This BIT expression is based on the assumption of a potential flow in which the influence of the bubbles on the liquid velocity fluctuations primarily results from the displacement of the liquid by the moving bubbles. Such an assumption, however, is not valid for real situations at finite bubble Reynolds number, where the wake contribution is dominant (Risso & Ellingsen Reference Risso and Ellingsen2002).
Other studies (Masood, Rauh & Delgado Reference Masood, Rauh and Delgado2014; Ullrich et al. Reference Ullrich, Krumbein, Maduta and Jakirlić2014) applied different SMC models conceived for single-phase flow in EE simulations of bubbly flow without accounting for any BIT effect. The configuration simulated in these references is the bubble plume experimentally investigated by Deen, Solberg & Hjertager (Reference Deen, Solberg and Hjertager2001) and good agreement was obtained. The good agreement in this particular case, however, does not necessarily mean that single-phase SMC models can be generalized to bubbly flows without any modification. The experimental data of Deen et al. (Reference Deen, Solberg and Hjertager2001) have also been reproduced very well using scale-resolving simulations without accounting for BIT (Deen et al. Reference Deen, Solberg and Hjertager2001; Niceno, Dhotre & Deen Reference Niceno, Dhotre and Deen2008; Ma et al. Reference Ma, Lucas, Ziegenhein, Fröhlich and Deen2015a). The cited papers show that, in this particular flow, the undulatory modulation of the bubble plume is the dominant feature, generating substantially more fluctuations than BIT. Hence, accounting for the latter or not affects the result only to a small extent (Ma et al. Reference Ma, Lucas, Ziegenhein, Fröhlich and Deen2015a). Further EE second-moment modelling for bubbly flows was performed by Cokljat et al. (Reference Cokljat, Slack, Vasquez, Bakker and Montante2006) and Colombo & Fairweather (Reference Colombo and Fairweather2015), and recently by Ullrich (Reference Ullrich2017). While Colombo & Fairweather (Reference Colombo and Fairweather2015) used a larger BIT source term in the streamwise direction, the latter two studies employed an isotropic source tensor to represent BIT. In all these studies, the extra BIT tensor in the SMC model was not derived on the basis of data but resulted from various modelling arguments or ad hoc physical considerations. The resulting expressions were subjected to indirect testing by calculating various bubbly flows with these closures.
The starting point of the present work is constituted by the exact Reynolds-stress equations for two-phase flows rigorously derived by Kataoka, Besnard & Serizawa (Reference Kataoka, Besnard and Serizawa1992). These equations are based on a single-phase representation, including the influence of the bubbles by additional so-called interfacial terms in the balance equations for each Reynolds-stress component and the dissipation rate of TKE. When supplementing the unclosed terms in these equations with adequate models, they constitute an appropriate framework for second-moment modelling of bubbly flows.
Recently, DNS data of sufficiently large-scale turbulent bubbly flows have become available (Bolotnov et al. Reference Bolotnov, Jansen, Drew, Oberai, Lahey and Podowski2011; Roghair et al. Reference Roghair, Mercado, Van Sint Annaland, Kuipers, Sun and Lohse2011; Dabiri, Lu & Tryggvason Reference Dabiri, Lu and Tryggvason2013; Lu & Tryggvason Reference Lu and Tryggvason2013; Santarelli & Fröhlich Reference Santarelli and Fröhlich2016; Cifani, Kuerten & Geurts Reference Cifani, Kuerten and Geurts2018), which can be used as a basis to test the individual model assumptions of EE RANS directly. Furthermore, the data can also be used to develop more elaborate closing approximations for BIT terms than formerly employed. Several works of this type have been accomplished in the framework of EE two-equation RANS modelling. Ilić (Reference Ilić2006) and Erdogan & Wörner (Reference Erdogan and Wörner2014) performed DNS studies with up to eight monodisperse bubbles rising in initially stagnant liquid within a plane channel and evaluated each term in the TKE budget. They demonstrated that the gain of TKE is mainly caused by the interfacial term, while the contribution of the single-phase production term is negligible. Besides evaluating the relative importance of each term in the TKE budget under these conditions, the DNS data were also used for a priori testing of available BIT models to assess the accuracy of them. Employing DNS of a much larger number of bubbles and more realistic physical parameters, such as the density ratio and bubble Reynolds number, Santarelli, Roussel & Fröhlich (Reference Santarelli, Roussel and Fröhlich2016) computed the budget terms in the TKE equation and proposed a BIT closure based on a priori evaluations. Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017) extended this work and for the first time used DNS to develop a complete BIT closure for a two-equation EE RANS approach, employing the full set of equations and performing a posteriori validations of the resulting expression.
To the best of the authors’ knowledge, DNS data so far have not been used to develop a realistic BIT model of bubbly flows in a differential SMC framework. Here, such an approach is used employing the data of Santarelli & Fröhlich (Reference Santarelli and Fröhlich2015, Reference Santarelli and Fröhlich2016) to devise a new SMC for the BIT. This term dominates in many bubbly flows, so that an improvement in this respect has considerable impact, as experienced in earlier studies (Ma et al. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017). In the work reported here, four main issues are investigated. First, the limitation of using the standard eddy-viscosity expression in the EE two-equation models for bubbly flow is discussed in § 2.2. The second issue is the model form of the EE SMC, particularly the pressure–strain term in the application of bubbly flows, which is addressed in § 3. Next, the coefficients in the BIT tensor expressions of the selected SMC are determined via a suitably chosen iterative procedure in § 4. Finally, based on an analysis of anisotropy invariants, the new SMC BIT model is proposed in § 5. The resulting SMC BIT model is then validated by computing the same cases and one case of Lu & Tryggvason (Reference Lu and Tryggvason2008) with the EE approach in § 6. Beyond the new model itself, the paper also furnishes a systematic procedure that is of general use for this type of modelling.
2 Background
2.1 Direct numerical simulation database
To be used for the development of closures, the DNS data employed have to provide resolved information about the respective processes. For the present work, this requires DNS with geometrically resolved bubbles. Extracting statistical information from such simulations about interfacial terms is by no means trivial (Santarelli et al. Reference Santarelli, Roussel and Fröhlich2016). Furthermore, to provide relevant data for modelling, the effect considered should be of sizeable, if not dominating, impact in the reference case selected. This second condition is also met with the flows considered here, which are governed by the balance between production of TKE through BIT and dissipation.
In Santarelli & Fröhlich (Reference Santarelli and Fröhlich2015, Reference Santarelli and Fröhlich2016) bubble-resolving DNS with many thousands of bubbles at low Eötvös number were conducted. Bubbles were considered as spherical objects of fixed shape and a no-slip condition was applied at the phase boundary, matching with the behaviour of air bubbles rising in contaminated water. Compared to other simulations of this type (Ilić Reference Ilić2006; Erdogan & Wörner Reference Erdogan and Wörner2014), these simulations are substantially closer to applications in that they involve turbulent background flow, contaminated fluid, realistic density ratio, higher bubble Reynolds number, a much larger domain and a much larger number of bubbles.
The DNS were conducted for upward vertical flow between two flat walls in a channel, with $x$ the streamwise, $y$ the wall-normal and $z$ the spanwise coordinate. The size of the computational domain is $L_{x}\times L_{y}\times L_{z}=4.41H\times H\times 2.21H$, where $H$ is the distance between the walls. Figure 1 shows the domain and an instantaneous snapshot of the bubbly flow in one of the DNS cases, labelled SmMany. A no-slip condition was applied at the walls and periodic conditions in $x$ and $z$. The gravitational force acts in the negative $x$-direction, and the bulk velocity $U_{b}$ was kept constant by instantaneously adjusting a volume force, equivalent to a pressure gradient, thus imposing a desired bulk Reynolds number $Re_{b}=U_{b}H/\unicode[STIX]{x1D708}$, where $\unicode[STIX]{x1D708}$ is the kinematic viscosity of the liquid. The DNS were all conducted with $Re_{b}=5263$. The data used in this work were obtained for three monodisperse cases (SmMany, SmFew, LaMany) and one bidisperse case labelled BiDisp, of the same void fraction as SmMany and LaMany with half the void fraction consisting of smaller bubbles and the other half of larger bubbles. Additionally, a single-phase simulation labelled Unladen was performed under the same conditions for comparison. Table 1 provides an overview of all cases with the corresponding labels. The data available cover statistical moments of first and second order for liquid and bubbles, as well as two-point correlations. The technically involved numerical procedure to evaluate the TKE budget was presented in Santarelli et al. (Reference Santarelli, Roussel and Fröhlich2016).
While the DNS were performed using a non-dimensional set of parameters, the EE SMC proposed in the present study is related to the size of bubble $d_{p}$, with the corresponding simulations performed in dimensional units. For this reason, the above set-up is converted to a dimensional form using the contaminated air–water system as an example. Based on the discussion in Ern et al. (Reference Ern, Risso, Fabre and Magnaudet2012), the pivoting element is chosen to be the equality of the Archimedes number, $Ar$, defined as
with $g$ the gravity. This leads to an accurate interpretation of the real effect of buoyancy in the considered DNS. Keeping $Ar$ the same as in the DNS and using all the other physical dimensional parameters on the right-hand side of (2.1) yields $d_{p}=1.456~\text{mm}$ for the smaller bubbles and $d_{p}=2.127~\text{mm}$ for the larger bubbles. The ratio $d_{p}/H$ in table 1 for the different cases then results in the extensions $123.6~\text{mm}\times 28.0~\text{mm}\times 61.8~\text{mm}$ of the channel.
2.2 Limitations of eddy-viscosity models for bubbly flow
Figure 2 clarifies a special shortcoming of all EE $k$–$\unicode[STIX]{x1D700}$ type models for bubbly flows constituting one of several motivations to develop the present EE SMC. Figure 2 was constructed from the DNS data (table 1) and shows an evaluation of the constant $C_{\unicode[STIX]{x1D707}}$, which can be expressed using the definition of the turbulent viscosity (Jones & Launder Reference Jones and Launder1972)
in terms of $k=\frac{1}{2}\overline{\overline{u_{i}^{\prime }u_{i}^{\prime }}}$ and $\unicode[STIX]{x1D700}=(1/\unicode[STIX]{x1D70C})\overline{\unicode[STIX]{x1D711}}\,\overline{\overline{\unicode[STIX]{x1D70F}_{ij}^{\prime }(\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j})}}$, with $\unicode[STIX]{x1D70F}_{ij}^{\prime }$ the fluctuating viscous stress tensor (Kataoka & Serizawa Reference Kataoka and Serizawa1989; Kataoka et al. Reference Kataoka, Besnard and Serizawa1992). Here, $\overline{\overline{\cdots \,}}$ denotes the phase-weighted averaging, defined by $\overline{\overline{F}}=\overline{F\unicode[STIX]{x1D711}}/\overline{\unicode[STIX]{x1D711}}$, where $\overline{\cdots \,}$ represents the Reynolds (or statistical) averaging with respect to time, space or ensemble of realizations. In these expressions, $\unicode[STIX]{x1D711}$ is an indicator function for the liquid phase, defined by $\unicode[STIX]{x1D711}(\boldsymbol{x},t)=1$ if $(\boldsymbol{x},t)$ in the liquid phase and equal to zero otherwise. Fluctuations of liquid quantities are defined as $F^{\prime }=F-\overline{\overline{F}}$.
As shown in Santarelli et al. (Reference Santarelli, Roussel and Fröhlich2016), the BIT term dominates for $y^{+}\gtrsim 30$, being equilibrated by the dissipation term only. The mean liquid velocity profile $\overline{\overline{u}}$ becomes flat in the channel centre, so that $\unicode[STIX]{x2202}\overline{\overline{u}}/\unicode[STIX]{x2202}y$ approaches zero and $\unicode[STIX]{x1D708}_{t}$ cannot be determined simply using the Boussinesq hypothesis as $\unicode[STIX]{x1D708}_{t}=-\overline{\overline{u^{\prime }v^{\prime }}}/(\unicode[STIX]{x2202}\overline{\overline{u}}/\unicode[STIX]{x2202}y)$. For this reason, $C_{\unicode[STIX]{x1D707}}$ in figure 2 is only plotted until $y^{+}=150$ for all the considered cases. Away from the wall, $C_{\unicode[STIX]{x1D707}}$ should trend to be a constant in single-phase channel flow (Pope Reference Pope2000). Indeed, in the region $60\leqslant y^{+}\leqslant 150$, the result of $C_{\unicode[STIX]{x1D707}}$ for the Unladen case matches well the result of Rodi & Mansour (Reference Rodi and Mansour1993) for the channel data of Kim, Moin & Moser (Reference Kim, Moin and Moser1987) at $Re_{\unicode[STIX]{x1D70F}}=180$, with $C_{\unicode[STIX]{x1D707}}\approx 0.13$, thus validating the present procedure. At the same time, it is noted that this value is higher than the standard value $C_{\unicode[STIX]{x1D707}}=0.09$ used for high-Reynolds-number single-phase flows.
For the bubble-laden cases, it could be expected that, with small gas void fraction, the flow retains most of the features from single-phase flow. This would imply that $C_{\unicode[STIX]{x1D707}}$ has a similar value as in the unladen flow. Figure 2 shows that this is not the case. For SmFew with the lowest void fraction, the largest deviation from single-phase flow is observed. Its corresponding curve in figure 2 shows $C_{\unicode[STIX]{x1D707}}$ to increase towards the channel centre, with $C_{\unicode[STIX]{x1D707}}\approx 0.3$ at $y^{+}=150$. For the other cases with higher gas void fraction (SmMany, LaMany and BiDisp), the behaviour of $C_{\unicode[STIX]{x1D707}}$ is very different case by case. In particular, for the case SmMany, the value of $C_{\unicode[STIX]{x1D707}}$ observed in the channel centre is approximately $0.13{-}0.2$ over more than four-fifths of the channel width. Surprisingly, $C_{\unicode[STIX]{x1D707}}$ in the most BIT-dominated case LaMany has the lowest value in all the bubble-laden cases and approaches the Unladen case. As expected, the case BiDisp produces the result between SmMany and LaMany.
The behaviour of $C_{\unicode[STIX]{x1D707}}$ can be discussed further via the distributions of $-\overline{\overline{u^{\prime }v^{\prime }}}/k$, the structure parameter, and the ratio of the shear-stress-induced production to the dissipation, $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}$, reflecting the local energy balance in single-phase flow. For channel flow, (2.2) then reads
In single-phase flow, $C_{\unicode[STIX]{x1D707}}=0.09$ corresponds to $-\overline{\overline{u^{\prime }v^{\prime }}}/k=0.3$ and $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}=1$ (local equilibrium). Figures 3(a) and 3(b) show these two parameters. For the Unladen case, the terms are identical to the results reported by Rodi & Mansour (Reference Rodi and Mansour1993) for the low-$Re$ single-phase channel flow, with $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}\approx 0.81$ in the region where $-\overline{\overline{u^{\prime }v^{\prime }}}/k=0.3$. In all the bubble-laden cases, the maxima of $-\overline{\overline{u^{\prime }v^{\prime }}}/k$ are shifted towards the walls and exhibit much lower values compared to the Unladen case. The curve of $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}$ shows that local equilibrium between $\unicode[STIX]{x1D617}_{k}$ and $\unicode[STIX]{x1D700}$ is not achieved for any of the bubble-laden cases. Furthermore, in the bubble-laden cases, $(\overline{\overline{u^{\prime }v^{\prime }}}/k)^{2}$ drops slower than $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}$ towards the channel centre, so that $C_{\unicode[STIX]{x1D707}}=(\overline{\overline{u^{\prime }v^{\prime }}}/k)^{2}/(\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700})$ (equation (2.3)) increases towards the channel centre (figure 2). Attention needs to be directed to the case SmFew, which, compared to the other bubble-laden cases, exhibits a shape similar to the Unladen case in both parameters $-\overline{\overline{u^{\prime }v^{\prime }}}/k$ and $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}$ (figure 3a,b). This is due to its low gas void fraction $(0.29\,\%)$, so that the influence of BIT is smaller. However, it does not necessarily mean that the value of $C_{\unicode[STIX]{x1D707}}$ should be closer to the Unladen case than the other cases, since $C_{\unicode[STIX]{x1D707}}$ is proportional to the ratio of $(\overline{\overline{u^{\prime }v^{\prime }}}/k)^{2}$ to $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}$ according to (2.3). Indeed, for $y^{+}>45$, a value of $C_{\unicode[STIX]{x1D707}}$ higher than in any other case is obtained in the SmFew case.
The above analysis illustrates that an EE $k$–$\unicode[STIX]{x1D700}$ type model for bubbly flows with a constant value of $C_{\unicode[STIX]{x1D707}}$ lacks realism and that, furthermore, employing the single-phase value $C_{\unicode[STIX]{x1D707}}=0.09$ can be off by a considerable factor – up to them in the present cases. As a result, the level of the eddy viscosity can be inappropriate and lead to deficient predictions. At a basic level, the problem is that the physical representation (2.3) for $C_{\unicode[STIX]{x1D707}}$ only includes a purely single-phase parameter, $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}$, to present the state of energy equilibrium, which is not suitable for BIT-dominated flows. As a result, it might be advantageous to optimize the determination of $C_{\unicode[STIX]{x1D707}}$ when employing an eddy-viscosity model.
Apart from the issue of choosing an expression for $\unicode[STIX]{x1D708}_{t}$ and a model constant, eddy-viscosity models cannot account for the anisotropic velocity fluctuations due to the buoyancy-generated rise of bubbles through the liquid. Incorporating this feature in an EE model should improve the results. These observations motivate the development of an EE SMC in the present paper.
3 Form of second-moment closure for bubbly flows
3.1 Basic equations
For incompressible gas–liquid two-phase flow without phase transition, the governing equations of the EE approach (Ishii & Hibiki Reference Ishii and Hibiki2006) are
with all quantities being mean values. Here, the superscript $K$ can be $L$ to designate liquid or $G$ to designate gas. Furthermore, $\overline{\overline{\boldsymbol{u}}}$, $\overline{\overline{\unicode[STIX]{x1D64E}}}$, $\unicode[STIX]{x1D70C}$, $\unicode[STIX]{x1D707}$ and $p$ represent the mean velocity, the mean strain-rate tensor, the density, the molecular dynamic viscosity and the pressure of the respective phase. The term $\overline{\overline{\unicode[STIX]{x1D749}}}_{t}$ results from the unresolved stress tensor. The sum of all interfacial forces acting on phase $K$ is termed $\boldsymbol{M}^{K}$ and needs to be modelled. Established models for the different non-drag interfacial forces are employed here, based on an extensive literature study. The model of Antal, Lahey & Flaherty (Reference Antal, Lahey and Flaherty1991) is used for the wall force, the model of Burns et al. (Reference Burns, Frank, Hamill and Shi2004) for turbulent dispersion, and the model of Auton (Reference Auton1987) for the lift force. Details on the model form and how the various coefficients of the models were selected are given in appendix A and § 4 below, respectively.
Kataoka et al. (Reference Kataoka, Besnard and Serizawa1992) proposed exact balance equations for the Reynolds stresses of the liquid phase, including the terms resulting from the interaction with a transported disperse phase,
where the terms on the right-hand side of (3.3) read
These four terms are analogous to the single-phase terms, just with the average void fraction $\overline{\unicode[STIX]{x1D711}}$ as a prefactor. The additional term
represents the interfacial energy transfer, where the index $L$ indicates that the respective quantity is evaluated at the liquid side of the phase boundary. Finally, $n_{i}$ is the normal vector at the phase boundary directed towards the gas phase and $I$ is the interfacial area concentration, with
3.2 Basic approach for closure
To close the single-phase terms in (3.3), the generalized SMC formulation of Hanjalić & Launder (Reference Hanjalić and Launder2011) is employed for the liquid phase, supplemented with a source tensor $\unicode[STIX]{x1D61A}_{R,ij}^{SMC}$ accounting for production of BIT,
where $\unicode[STIX]{x1D6FC}^{L}=\overline{\unicode[STIX]{x1D711}}$ denotes the mean liquid void fraction. Here and in the following all average quantities refer to the liquid, so that the upper index $L$ is dropped from now on for clarity. The shear-induced production of TKE, $\unicode[STIX]{x1D617}_{ij}^{SMC}$, only comprises Reynolds stresses and mean flow gradients, so that in this framework no closure is required.
The diffusion term $\unicode[STIX]{x1D60B}_{ij}^{SMC}$ is approximated by the gradient diffusion hypothesis of Shir (Reference Shir1973), reading
with $c_{s}=0.22$ and $\unicode[STIX]{x1D708}^{L}$ the liquid molecular kinematic viscosity. For the dissipation term, local isotropy of $\unicode[STIX]{x1D700}_{ij}$ is assumed, which is widely used in SMCs for single-phase flow (Pope Reference Pope2000; Hanjalić & Jakirlić Reference Hanjalić, Jakirlić, Launder and Sandham2002; Hanjalić & Launder Reference Hanjalić and Launder2011). It yields
The modelling strategy, here, is to absorb any departure from isotropy of the dissipation processes in the model for $\unicode[STIX]{x1D719}_{ij}$, as in single-phase flow. To some extent, this is of no consequence, since the anisotropic part $\unicode[STIX]{x1D700}_{ij}-\unicode[STIX]{x1D6FC}^{L}\frac{2}{3}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D700}$ has the same mathematical properties as $\unicode[STIX]{x1D719}_{ij}$ and $\unicode[STIX]{x1D61A}_{R,ij}$, so that it can be absorbed in the models established for them (Pope Reference Pope2000).
The turbulent energy dissipation rate $\unicode[STIX]{x1D700}$ is obtained from its own transport equation,
where all single-phase terms are modelled by the standard form and all constants are given in appendix B. The BIT source term in the $\unicode[STIX]{x1D700}$ equation, $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D700}}^{SMC}$, was modelled in Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017) and is employed here without modification. In this expression, $\unicode[STIX]{x1D61A}_{k}^{SMC}$, the interfacial term for the $k$ equation (i.e. $\unicode[STIX]{x1D61A}_{k}^{SMC}=S_{k}^{k-\unicode[STIX]{x1D700}}$) proposed by Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017), is
with $\boldsymbol{F}_{D}$ the drag and
the time scale characterizing BIT. This time scale is based on the energy spectra analysis and yields more convincing evidence (Ma et al. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017).
The goal of the following is to model the pressure–strain term $\unicode[STIX]{x1D719}_{ij}$ and the interfacial term $\unicode[STIX]{x1D61A}_{R,ij}$ using the DNS data. For the development of the BIT models, it is appropriate to focus on the channel centre, since the interfacial term is balanced by the pressure–strain and dissipation term there, so these three terms can be considered in isolation from other effects, while the shear-production and diffusion terms are negligible.
3.3 Treatment of the pressure–strain correlations
Apart from the interfacial term, the pressure–strain term $\unicode[STIX]{x1D719}_{ij}$ is the only correlation that contains directional information, so that it plays a pivotal role in capturing the Reynolds-stress anisotropy in bubbly flow. In the present work, only turbulence in the liquid phase is considered, with disperse bubbles considered as momentum sources of finite size distributed in the continuous phase. Hence, the closure for the pressure–strain term $\unicode[STIX]{x1D719}_{ij}$ can broadly be taken over from the existing ones for turbulence affected by external force fields, such as buoyancy (Launder Reference Launder1975), rotation (Launder, Tselepidakis & Younis Reference Launder, Tselepidakis and Younis1987) and electromagnetic forces (Kenjereš, Hanjalić & Bal Reference Kenjereš, Hanjalić and Bal2004) in single-phase flow. The term is made up of contributions related to manipulation of the exact Poisson equation for the pressure fluctuations in the short-hand form
with each term on the right-hand side related to a different physical process of isotropization: turbulence self-interactions $\unicode[STIX]{x1D719}_{ij,1}$ (the slow term), strain production $\unicode[STIX]{x1D719}_{ij,2}$ (the rapid term), bubble-induced force production $\unicode[STIX]{x1D719}_{ij,3}$ and wall blocking $\unicode[STIX]{x1D719}_{ij,w}$ (Hanjalić & Launder Reference Hanjalić and Launder2011).
Further analysis focuses on the two dominant contributors $\unicode[STIX]{x1D719}_{ij,1}$ and $\unicode[STIX]{x1D719}_{ij,3}$ remote from the wall in the present BIT-dominated cases. Evidently, the related physical character of the slow term $\unicode[STIX]{x1D719}_{ij,1}$ in bubbly flow should not differ from that in single-phase flow. It redistributes energy among the components and diminishes any shear stress, thus causing turbulence ‘slowly’ to approach its isotropic state (Hanjalić & Launder Reference Hanjalić and Launder2011). Hence, modelling this term can safely start from the most general formulation in single-phase flow according to the Cayley–Hamilton theorem (Lumley Reference Lumley1978; Fu, Launder & Tselepidakis Reference Fu, Launder and Tselepidakis1987; Speziale, Sarkar & Gatski Reference Speziale, Sarkar and Gatski1991)
where $\unicode[STIX]{x1D622}_{ij}$ is the Reynolds-stress anisotropy tensor
and $A_{2}=\unicode[STIX]{x1D622}_{ji}\unicode[STIX]{x1D622}_{ij}$ is the second invariant of $\unicode[STIX]{x1D622}_{ij}$. The factors $c_{1}$ and $c_{1}^{\prime }$ are two model coefficients, being either constants or functions of the turbulence Reynolds number and stress anisotropy invariants. Most model proposals for the slow term in the literature have the structure of (
3.17), differing in the model coefficients $c_{1}$ and $c_{1}^{\prime }$ (see e.g. Launder et al. Reference Launder, Reece and Rodi1975; Gibson & Launder Reference Gibson and Launder1978; Speziale et al. Reference Speziale, Sarkar and Gatski1991; Ristorcelli, Lumley & Abid Reference Ristorcelli, Lumley and Abid1995). More related models and the corresponding detailed forms of the model coefficients can be found in the summary of Jakirlić & Hanjalić (Reference Jakirlić and Hanjalić2013). Among these, most models for the slow term consider only the linear part, i.e. $c_{1}^{\prime }=0$ in (
3.17), which was first proposed by Rotta (Reference Rotta1951) and later adopted in the popular LRR model. However, the importance of the nonlinear part of the slow term has been mentioned by many studies (e.g. by Lumley Reference Lumley1978; Fu et al. Reference Fu, Launder and Tselepidakis1987; Speziale et al. Reference Speziale, Sarkar and Gatski1991; Hanjalić & Jakirlić Reference Hanjalić, Jakirlić, Launder and Sandham2002). Jakirlić & Hanjalić (Reference Jakirlić and Hanjalić2013) scrutinized the model coefficients comprehensively using single-phase DNS channel data and found that the standard coefficients $c_{1}=1.7$ and $c_{1}^{\prime }=1.05$ proposed by Sarkar & Speziale (Reference Sarkar and Speziale1990) and Speziale et al. (Reference Speziale, Sarkar and Gatski1991) – hereafter referred to as SSG – are in very good agreement with the DNS, so that these are used here.
To assess the relative importance of the nonlinear part in the slow term for modelling bubbly flow, both the nonlinear term $\unicode[STIX]{x1D719}_{ij,1,nonlinear}$ and the linear term $\unicode[STIX]{x1D719}_{ij,1,linear}$ were evaluated from (3.17) in an a priori manner using the DNS data of the case SmMany to compute $\unicode[STIX]{x1D6FC}^{L}$, $\unicode[STIX]{x1D700}$, $\unicode[STIX]{x1D622}_{ij}$ and $A_{2}$. The results are shown in figure 4.
The rapid term $\unicode[STIX]{x1D719}_{ij,2}$ is modelled here with the general isotropization-of-production (IP) model of Naot, Shavit & Wolfshtein (Reference Naot, Shavit and Wolfshtein1970),
with $c_{2}=0.6$. This expression can be evaluated as well using $\unicode[STIX]{x1D617}_{ij}$ from the DNS data, and the result of this expression for SmMany is reported in figure 4 to give an impression of the corresponding magnitude in the present bubbly flow. A perfect reference to assess the relative contribution of the particular terms $\unicode[STIX]{x1D719}_{ij,1,nonlinear}$, the linear term $\unicode[STIX]{x1D719}_{ij,1,linear}$ and $\unicode[STIX]{x1D719}_{ij,2}$ requires DNS data for $\unicode[STIX]{x1D719}_{ij}$. However, as mentioned in § 2, no Reynolds-stress budget is available from the DNS, but a TKE budget is available. For this reason, the interfacial term of the TKE equation, $S_{k}\equiv \frac{1}{2}\unicode[STIX]{x1D61A}_{R,ii}$, is used here as a reference to decide about the relevance of these terms in the Reynolds-stress budget (3.3). The term $S_{k}$ is the main energy input in the present flow, so that its value is of the same order as $\unicode[STIX]{x1D719}_{ij}$ and $\unicode[STIX]{x1D700}_{ij}$.
For the channel centre in all normal components (figure 4a–c), the nonlinear part of the slow term $\unicode[STIX]{x1D719}_{ij,1,nonlinear}$ is smaller but of the same order as the linear part $\unicode[STIX]{x1D719}_{ij,1,linear}$ and has overall the opposite sign to the latter. The influence of the rapid term $\unicode[STIX]{x1D719}_{ij,2}$ appears to be negligible for this BIT-dominated flow in the channel centre for all components. This is expected, as it is related to the production term $\unicode[STIX]{x1D617}_{ij}$, which is small in the centre due to the vanishing mean flow gradients. All contributions to the pressure–strain term have a much smaller magnitude in the Reynolds shear stress component. Among these, the linear part of the slow term is dominant in the channel centre (see figure 4d). The behaviour of the total term $\unicode[STIX]{x1D719}_{ij}$ will be discussed in a later section.
The analysis of the pressure–strain terms for the other three cases (SmFew, LaMany and BiDisp) reveals the same trends as in figure 4, so that they are not reproduced here. The observations based on the DNS data hence support two guidelines for modelling the pressure–strain term in bubbly flow. First, due to the importance of the nonlinear part $\unicode[STIX]{x1D719}_{ij,1,nonlinear}$, the slow term $\unicode[STIX]{x1D719}_{ij,1}$ requires a nonlinear model. Second, it indicates that the rapid term is very small, so that the specific form of the model is uncritical and the simple model (3.19) sufficient. Furthermore, more elaborate models extending (3.19) are heavily based on the mean flow gradient (Johansson & Hallbäck Reference Johansson and Hallbäck1994), which is of less interest here.
Finally, the idea underlying the IP model can also be applied to $\unicode[STIX]{x1D719}_{ij,3}$ to consider the influence of the bubble-induced force field on the pressure–strain correlation,
with $c_{3}$ a model coefficient and $\unicode[STIX]{x1D61A}_{R,ij}^{SMC}$ the BIT source term in the Reynolds-stress equation (3.10). This term tends to redistribute the action of the interfacial term, reducing the net BIT production in the rich component.
3.4 Closure of the interfacial term
Closure of (3.10) as well as (3.20) requires a model for $\unicode[STIX]{x1D61A}_{R,ij}$. The first step is to express the model of the interfacial term, $\unicode[STIX]{x1D61A}_{R,ij}^{SMC}$, in terms of the algebraic expression $b_{ij}$ and half the trace $\unicode[STIX]{x1D61A}_{k}^{SMC}\equiv \frac{1}{2}\unicode[STIX]{x1D61A}_{R,ii}^{SMC}$ given by (3.14), i.e.
with $b_{ij}$ to be determined to consider the Reynolds-stress anisotropy in bubbly flow. The second modelling step is to assume that the diagonal BIT terms $\unicode[STIX]{x1D61A}_{R,ij}$, $i=j$, dominate over the off-diagonal terms, $\unicode[STIX]{x1D61A}_{R,ij}$, $i\neq j$, so that the latter can be set to zero. This is based on the fact that accurate reproduction of Reynolds normal stress $\overline{\overline{u_{i}^{\prime }u_{j}^{\prime }}}$, $i=j$, is the main prerequisite to capture the Reynolds-stress anisotropy. Reynolds shear stresses, on the other hand, while themselves being a measure of turbulence anisotropy, only reflect the effect of shear straining, which is of lower importance (Jakirlić & Hanjalić Reference Jakirlić and Hanjalić2013). The third step is to observe that the source term $\unicode[STIX]{x1D61A}_{R,ij}$ differs substantially between the vertical direction, here $i=1$, i.e. the direction of bubble rise, and the other two directions. This is backed by experiments and simulations (Akbar et al. Reference Akbar, Hayashi, Hosokawa and Tomiyama2012; Santarelli & Fröhlich Reference Santarelli and Fröhlich2015). Furthermore, in the channel far from the walls, it is reasonable to assume isotropy in both cross-stream directions for geometrical reasons. This justifies setting $b_{22}=b_{33}$, which also is in line with the cited reference data.
Finally, having the same attribute as the slow term and the rapid term, $\unicode[STIX]{x1D719}_{ij,3}^{SMC}$ employed with IP model (3.20) is also traceless, i.e. $\unicode[STIX]{x1D719}_{ii,3}^{SMC}=0$. When considering $\unicode[STIX]{x1D61A}_{R,ij}^{SMC}$ based on the second modelling step being only with the diagonal terms, it is possible to assign $\unicode[STIX]{x1D719}_{ij,3}^{SMC}$ to $\unicode[STIX]{x1D61A}_{R,ij}^{SMC}$ and, considered as a whole, defined as ‘effective BIT source’, i.e.
The models for $\unicode[STIX]{x1D719}_{ij,1}$ and $\unicode[STIX]{x1D719}_{ij,2}$ remain untouched. Again, $b_{22}^{\ast }=b_{33}^{\ast }$. Introducing the information about the trace, here the facts that $b_{ii}\equiv 2$ and $\unicode[STIX]{x1D719}_{ii,3}^{SMC}=0$ result in $b_{ii}^{\ast }=2$. This yields
All the assumptions made above and the concept to import the ‘effective BIT source’ $\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff}$ succeeds in substantially reducing the complexity of the resulting formulations. In fact, providing a BIT closure now amounts to providing an expression for $b_{11}^{\ast }$. Determining $c_{3}$ in (3.20) is no longer needed.
3.5 Proposal for multi-disperse bubble swarms
One of the applications below features a bidisperse bubble swarm, so that the extension of the proposed approach to this situation and to multi-disperse bubbles in general is briefly discussed.
For BIT modelling of bidisperse bubble swarms in the framework of the two-equation RANS model, linear superposition of the source term for each gas group is used in the traditional way, in both the TKE equation and the dissipation equation (e.g. by Politano, Carrica & Converti Reference Politano, Carrica and Converti2003; Ziegenhein et al. Reference Ziegenhein, Rzehak, Ma and Lucas2017; Liao et al. Reference Liao, Ma, Liu, Ziegenhein, Krepper and Lucas2018, Reference Liao, Ma, Krepper, Lucas and Fröhlich2019). The current SMC for the multi-disperse case extends this concept with the assumption that the contributions from the two bubble sizes sum up to the total contribution. In the present BiDisp case (table 1), this reads
for the effective BIT source terms in the Reynolds-stress equations and
for the source term in the dissipation equation, respectively.
4 Determination of model parameters for the different cases
4.1 Procedure for determination of model coefficients
Based on the analysis of figure 4, EE SMC simulations with $\unicode[STIX]{x1D719}_{ij}^{SMC}$ from the SSG model (appendix B) were run in the same domain as the DNS and discretized with $56\times 60\times 51$ grid points in the $x$-, $y$- and $z$-directions, respectively, according to the mesh studies to supply the best prediction. A uniform grid is selected here, considering the different gas distributions in the present DNS cases. At the walls, a no-slip boundary condition was applied for the continuous phase and a free-slip condition for the disperse phase. The other boundary conditions were identical to those of the DNS. The liquid and gas are defined as homogeneously distributed at the beginning in the entire computed domain, representing the DNS set-up, e.g. with $\unicode[STIX]{x1D6FC}^{L}=0.9786$ and $\unicode[STIX]{x1D6FC}^{G}=0.0214$ for the case SmMany. Then a constant mass flow rate for SmMany, $0.283~\text{kg}~\text{s}^{-1}$, is calculated from ${\dot{m}}=U_{b}(\unicode[STIX]{x1D70C}^{L}\unicode[STIX]{x1D6FC}^{L}+\unicode[STIX]{x1D70C}^{G}\unicode[STIX]{x1D6FC}^{G})A_{yz}$, with $A_{yz}=L_{y}\times L_{z}$ being the cross-section of the channel. The same procedure is performed for the other cases as well. To perform the simulations, the commercial code ANSYS CFX v17.2 was used. For the spatial discretization, a high-resolution scheme (Barth & Jesperson Reference Barth and Jesperson1989) is employed, and a second-order backward Euler scheme is used in time. The new model was implemented by means of a set of user-defined functions (UDFs). (These UDFs are available from the authors upon request.)
The task now is to determine the coefficient $b_{11}^{\ast }$ in (3.22) for each DNS dataset, which later on should reveal the general trend. To accomplish this, one prerequisite is to determine the effective BIT source $\unicode[STIX]{x1D61A}_{R,ij}^{eff}$ from the DNS data constituting the target for the corresponding modelling term $\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff}$. The data of Santarelli et al. (Reference Santarelli, Roussel and Fröhlich2016) comprise mean flow, Reynolds stresses and budget terms of $k$, but no budget of $\overline{\overline{u_{i}^{\prime }u_{j}^{\prime }}}$. Hence, $\unicode[STIX]{x1D61A}_{R,ij}^{eff}$ is evaluated indirectly from these data focusing in the core region of the channel where BIT effects dominate. In this region, (3.3) reduces to
The modified pressure–strain term $\unicode[STIX]{x1D719}_{ij}^{mod}=\unicode[STIX]{x1D719}_{ij}-\unicode[STIX]{x1D719}_{ij,3}=\unicode[STIX]{x1D719}_{ij,1}+\unicode[STIX]{x1D719}_{ij,2}$ is represented here by the sum of $\unicode[STIX]{x1D719}_{ij,1}^{SMC}$ (3.17) and $\unicode[STIX]{x1D719}_{ij,2}^{SMC}$ (3.19), where $\unicode[STIX]{x1D6FC}^{L}$, $\unicode[STIX]{x1D700}$, $\unicode[STIX]{x1D622}_{ij}$, $A_{2}$ and $\unicode[STIX]{x1D617}_{ij}$ are determined from the DNS data,
using the standard coefficients $c_{1}=1.7$, $c_{1}^{\prime }=1.05$ and $c_{2}=0.6$. In an analogous manner, the dissipation term is evaluated as
Equation (
4.1) then provides the effective BIT source $\unicode[STIX]{x1D61A}_{R,ij}^{eff}=-(\unicode[STIX]{x1D700}_{ij}+\unicode[STIX]{x1D719}_{ij}^{mod})$.
The present approach strives to develop the BIT model not as a stand-alone procedure by means of a priori tests, as done by (Santarelli et al. Reference Santarelli, Roussel and Fröhlich2016), but in the framework of the complete EE SMC model. This procedure is now described. The entire EE model is complex, with different submodels influencing the result. When working on one of these submodels, the uncertainty generated by the other terms has to be minimized. For this reason, it is mandatory that the bubble distribution over the domain corresponds to the reference, for example, since otherwise the BIT model cannot be assessed. This is achieved by optimizing $C_{L}$ and $C_{D}$ in the EE SMC model. For this purpose, the drag coefficient $C_{D}$ (table 1) is determined for each simulation from the balance between drag force and buoyancy force, according to
The relative velocity $u_{r}$ is obtained from the difference of the mean gas velocity and the mean liquid velocity as justified at length by Santarelli & Fröhlich (Reference Santarelli and Fröhlich2015). The same is not possible for the lift coefficient $C_{L}$. Instead, this coefficient is optimized such that the resulting void fraction distribution matches the DNS target, while the other interfacial force models are employed as they are (appendix A). Since these computations are fast, this is done manually using a bisection procedure. Once the correct void distribution is achieved, the model coefficient of the SMC, $b_{11}^{\ast }$, is determined. Changing the BIT model, however, generally results in a modified void distribution, so that the previous steps are repeated until convergence, as summarized in table 2. The originality here is that the targets of the optimization procedure are not statistics of the Reynolds stresses but rather the DNS values of the particular terms to be closed, like $\unicode[STIX]{x1D61A}_{R,11}^{eff}$ and $S_{\unicode[STIX]{x1D700}}$.
4.2 Closure for the case SmMany
For the case SmMany, the procedure of table 2 yields $C_{L}=0.06$ and $b_{11}^{\ast }=1.8$. The results obtained with these values match the DNS data very well. The void fraction in figure 5(a) is identical in the centre and only the near-wall peaks are somewhat more rounded. The mean velocities of gas and liquid also agree very well. Very close to the wall the DNS data exhibit a discretization phenomenon as the gas velocity was determined for the centre of the bubbles, which is at least the bubble radius off the wall. Figure 6 shows that also the Reynolds stresses from the EE SMC agree very well with the reference data in the channel centre. It is important to note that there is a factor of approximately eight between the magnitude of the streamwise and the wall-normal fluctuations and another factor of approximately eight to the turbulent shear stress. This very satisfactory agreement validates not only the choice of the iteratively determined model constants but also the functional form (3.22) and all the other ingredients in the model equations (3.10) and (3.13). Indeed, the achievement can be highlighted by comparison with the results obtained using other standard models from the literature under exactly the same conditions. Cokljat et al. (Reference Cokljat, Slack, Vasquez, Bakker and Montante2006) proposed an isotropic BIT source term, i.e.
Colombo & Fairweather (Reference Colombo and Fairweather2015) suggested
However, these relations do not account for the influence of the bubble-induced force field on the pressure–strain correlation (3.20). Their expressions are used here employing the same overall procedure of table 2 for the optimal comparison. It is obvious that, although the shapes of the profiles are similar, the level of the normal stresses is substantially different and bubble-induced anisotropy is not correctly captured.
The reasons for the different behaviour of the BIT models can be seen in the following figures. Using the DNS data, figure 7 shows the essential contributions to the budgets of the Reynolds normal stress $\overline{\overline{u^{\prime }u^{\prime }}}$ with its estimated dissipation term (4.3), modified pressure–strain term (4.2) and effective interfacial term (4.1). Here, the modified pressure–strain term and the dissipation rate constitute the two main sinks for the $\overline{\overline{u^{\prime }u^{\prime }}}$ component. In the channel centre, the budget of $\overline{\overline{u^{\prime }u^{\prime }}}$ reduces to $0\approx \unicode[STIX]{x1D700}_{11}+\unicode[STIX]{x1D719}_{11}^{mod}+\unicode[STIX]{x1D61A}_{R,11}^{eff}$, equation (4.1). This is backed up by figure 7(a), with $\unicode[STIX]{x1D61A}_{R,11}^{eff}$ the dominant source to compensate nearly the sum of the dissipation and the modified pressure–strain. For modelling $\unicode[STIX]{x1D61A}_{R,11}^{eff}$, the target is defined in the core region of the channel as remarked before, since capturing this correctly is a prerequisite for all models. With the present value of $b_{11}^{\ast }$, $\unicode[STIX]{x1D61A}_{R,11}^{SMC\text{-}eff}$ matches the value of the effective interfacial term in the DNS extremely well. The BIT models of Cokljat et al. (Reference Cokljat, Slack, Vasquez, Bakker and Montante2006) and Colombo & Fairweather (Reference Colombo and Fairweather2015) obviously predict significantly lower levels of $\unicode[STIX]{x1D61A}_{R,11}^{eff}$, which in turn leads to the underprediction of $\overline{\overline{u^{\prime }u^{\prime }}}$ and overprediction of $\overline{\overline{v^{\prime }v^{\prime }}}$ observed in figure 6.
Validating individual terms of the $\unicode[STIX]{x1D700}$ equation (3.13) is not possible since these are not available from the DNS data. Still, some validation can be performed by observing that
in the channel centre for the same reasons as with the budget of $\overline{\overline{u^{\prime }u^{\prime }}}$. Hence, $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D700}}$ can be determined indirectly from the DNS data. For this purpose, it is represented here by the expression
using the standard coefficient $C_{\unicode[STIX]{x1D700}2}=1.83$ (Hanjalić & Launder Reference Hanjalić and Launder2011). Equation (4.7) then results in an expression for $S_{\unicode[STIX]{x1D700}}$, with $S_{\unicode[STIX]{x1D700}}=-\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D700}}$. Figure 7(c) shows that the modelled BIT source of the $\unicode[STIX]{x1D700}$ equation, $S_{\unicode[STIX]{x1D700}}^{SMC}$, indeed, provides a very accurate result.
Note, as an aside, that figure 7(b), depicting the dominant components of the budget of the Reynolds normal stress $\overline{\overline{v^{\prime }v^{\prime }}}$, confirms the discussion on the determination of $b_{ij}^{\ast }$, $i=j$, above. Here, the value $b_{22}^{\ast }=0.1$ is obtained not from matching $\unicode[STIX]{x1D61A}_{R,22}^{eff}$, but rather by imposing $b_{ii}^{\ast }=2$ and $b_{22}^{\ast }=b_{33}^{\ast }$, resulting in $b_{22}^{\ast }=\frac{1}{2}(2-b_{11}^{\ast })$. It is worth underlining that the increase in $\overline{\overline{v^{\prime }v^{\prime }}}$ is due to energy being extracted from the $\overline{\overline{u^{\prime }u^{\prime }}}$ component and fed to the $\overline{\overline{v^{\prime }v^{\prime }}}$ component by pressure–strain interaction. Here, in contrast to the $\overline{\overline{u^{\prime }u^{\prime }}}$ component, the modified pressure–strain term, $\unicode[STIX]{x1D719}_{22}^{mod}$, appears as a source in the budget for $\overline{\overline{v^{\prime }v^{\prime }}}$, and the corresponding effective interfacial term $\unicode[STIX]{x1D61A}_{R,22}^{eff}$ exhibits very small influence in the centre.
4.3 Closure for the other monodisperse channel flow cases
The procedure of table 2 is now employed for the other two monodisperse bubble-laden channel flows of table 1, namely the cases SmFew and LaMany. Overall, this procedure achieves good agreement for the gas void fraction (figure 8a,c), the Reynolds stresses away from the wall (figure 8b,d) and the other parameters (not shown here to avoid redundancy), yielding the corresponding values of $b_{11}^{\ast }$ and $C_{L}$ reported in table 3. Close to the walls, the Reynolds normal stress profiles are underpredicted by the present SMC computations. This deficiency remains for description of channel flows that are intrinsic to high-$Re$ RANS models for both single-phase (Hanjalić & Launder Reference Hanjalić and Launder1976; Wilcox Reference Wilcox1998) and multiphase flows (Ma Reference Ma2017), and, hence, cannot be overcome by the mere introduction of a more complex BIT model.
4.4 Closure for bidisperse case
In particular for the case BiDisp, it is important to reproduce the profiles of the gas void fraction for smaller bubbles and larger bubbles, individually, not only the profile of the total gas void fraction. This requires that the simulation should be run with different momentum equations describing the predefined gas groups, one group for the smaller bubbles and the other for the larger bubbles. The EE BIT modelling for BiDisp is performed as described in § 3.5 with linear superposition of the BIT source term for each gas group, in both $\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff}$ (3.24) and $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D700}}^{SMC}$ (3.25).
The interfacial force setting follows § 4.1 and the optimization procedure (table 2) is carried out separately for each phase. In particular, the drag coefficients for small and large bubbles are calculated separately using (4.4) based on DNS data for each bubble size and the lift coefficients $C_{L}$ are adjusted to fit the void fraction profiles for smaller and larger bubbles, individually (figure 9a). The agreement achieved in the two void fractions is very good (figure 9a), so that also the total void fraction is well matched (not shown here). The good results for the Reynolds stresses (figure 9b) confirm that the BIT modelling according to § 3.5 is valid also for this bidisperse case. The corresponding values of $b_{11}^{\ast }$ and $C_{L}$ are listed in table 3 as well, together with those of the other three monodisperse cases.
4.5 A case without bulk velocity
To enhance the dataset over a wider range of bubble Reynolds numbers, further experimental data are used for a case that is characterized by the dominance of BIT as well, here the case reported by Akbar et al. (Reference Akbar, Hayashi, Hosokawa and Tomiyama2012). The Akbar case features a rectangular water–air bubble column, with a gas superficial velocity of $3~\text{mm}~\text{s}^{-1}$. Its height $(x)$, width $(y)$ and depth $(z)$ are 800 mm, 240 mm 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 %, which is around half of the value in the case SmMany. However, the bubble diameter is much larger. Furthermore, this case is very close to monodisperse with $d_{p}=4.37~\text{mm}$ (Akbar et al. Reference Akbar, Hayashi, Hosokawa and Tomiyama2012). More details are provided in the cited reference.
The present EE SMC simulation was performed on a $175\times 240\times 36$ grid in the $x$-, $y$- and $z$-directions, respectively. Based on our own previous work, the same set of interfacial force models as detailed by Liao et al. (Reference Liao, Ma, Krepper, Lucas and Fröhlich2019) were employed to achieve good agreement of gas void fraction and relative velocity. However, the experimental database does not provide any information about budgets of turbulence parameter, so that the procedure of table 2 cannot be applied directly here to yield $b_{11}^{\ast }$. Instead, this is achieved by adjustments of $b_{11}^{\ast }$ in the source term to match the Reynolds stresses measured in the experiment. Such a procedure is allowed for the present BIT closure at second-moment level for the following reason. If the complete set of the SMCs (six Reynolds-stress equations and one dissipation equation) for a BIT-dominated flow are fixed, with $b_{11}^{\ast }$ the only unknown, since $b_{22}^{\ast }=b_{33}^{\ast }=\frac{1}{2}(1-b_{11}^{\ast })$ is imposed, there appears to be only one value of $b_{11}^{\ast }$. For the Akbar case this procedure yields the value $b_{11}^{\ast }=1.4$ suitable to fit the resulting Reynolds stress $\overline{\overline{u^{\prime }u^{\prime }}}$. Since the form of the BIT source in the $\unicode[STIX]{x1D700}$ equation is fixed, this results in corresponding fixed $\unicode[STIX]{x1D700}_{ij}$ in the Reynolds-stress equations.
Note that the present approach is exempt from the pitfalls identified by the two-equation RANS model using ad hoc models targeting TKE only. The key difference is that there the BIT source of the $\unicode[STIX]{x1D700}$ equation is not fixed as detailed in Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017). Figure 10(a,b) shows a comparison of the simulated and experimental gas distribution, vertical liquid velocity and vertical gas velocity. All profiles were 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~\text{mm}$). The agreement of the mean flow profiles obtained with the present SMC closure is very satisfactory and improves substantially upon other models, such as our own previous work (Ma et al. Reference Ma, Ziegenhein, Lucas, Krepper and Fröhlich2015b). The flatter profile for the liquid velocity in figure 10(b) may be due to the slight inaccuracy in the gas void fraction close to the wall and using the high-$Re$ SMC version. Concerning the liquid Reynolds stresses, the present BIT model provides accurate results, reproducing the ratio of the Reynolds normal stresses $\overline{\overline{u^{\prime }u^{\prime }}}:\overline{\overline{v^{\prime }v^{\prime }}}\sim 2.5:1$ observed in the experiment for the centre region of the measurement line, as shown in figure 10(c,d).
5 Proposed second-moment closure for bubble-induced turbulence
5.1 Anisotropy invariants
It is common practice when developing SMC models to analyse the anisotropy of the flows. As argued by Lumley & Newman (Reference Lumley and Newman1977), the local state of the Reynolds-stress anisotropy tensor $\unicode[STIX]{x1D622}_{ij}$ can be usefully characterized by its second ($A_{2}=\unicode[STIX]{x1D622}_{ji}\unicode[STIX]{x1D622}_{ij}$) and third ($A_{3}=\unicode[STIX]{x1D622}_{ij}\unicode[STIX]{x1D622}_{jk}\unicode[STIX]{x1D622}_{ki}$) invariants. These two anisotropy invariants can certainly be combined in various ways to generate further invariants, for example, a useful third norm, known as Lumley’s flatness parameter (Lumley Reference Lumley1978),
In isotropic turbulence, both $A_{2}$ and $A_{3}$ vanish, and $A=1$. The other extreme condition corresponds to the two-component limit, where $A_{2}=A_{3}+8/9$, so that $A=0$.
To give some impression of how the above invariants may be distributed in a bubbly flow, figure 11 shows a comparison between the case SmMany and the case Unladen, both obtained from the DNS data. The flatness parameter $A$ vanishes at the wall in both cases as turbulence is in the two-component limit for geometrical reasons. In contrast to the single-phase situation, when $A$ increases towards the centre and approaches unity, the value of $A$ is around $0.25$ almost over the entire channel in the case SmMany, indicating a very anisotropic turbulence state that is generated by the bubbles. This state can also be identified by a much larger second invariant $A_{2}$ observed in this case, which gives a direct measure of the magnitude of the stress anisotropy. The third invariant $A_{3}$ has a similar behaviour as the corresponding value of $A_{2}$ in both cases and broader flatter profiles are observed in the SmMany case. Similar trends as for SmMany are identified in the other bubble-laden cases as well (not shown here), implying a very different nature of stress anisotropy compared to the corresponding single-phase flow.
For a more detailed analysis, it was shown by Lumley & Newman (Reference Lumley and Newman1977) that the realizable states for the stress anisotropy $\unicode[STIX]{x1D622}_{ij}$ should satisfy the following inequalities:
These inequalities give rise to the anisotropy-invariant map shown in figure 12. The upper line corresponds to the two-component limit (the last equality in (5.2)), the left-hand line to the axisymmetric contraction limit (the first equality in (5.2), with $A_{3}<0$), and the right-hand line to the axisymmetric expansion limit (the first equality in (5.2), with $A_{3}>0$).
Figure 12 shows the points for $\unicode[STIX]{x1D622}_{ij}$ derived from DNS and experiment on the anisotropy-invariant map for different $y^{+}$ positions over all the cases considered (the centre point of the measurement line is plotted for the Akbar case). Profiles for the bubble-laden cases and the single-phase case are plotted separately for better comparison. Additionally, the points obtained from the case SmFew would overlap the points from SmMany, so that it is plotted in figure 12(b) separately, as well. A first observation is that all states from both single-phase and bubbly flows, indeed, stay within the Lumley triangle, as is required by realizability constraints. The main difference between the two is that the path crosses the channel in the wall-normal direction. Among them, the Unladen case (figure 12c) shows a traverse from the wall towards the channel centre – starting in the two-component limit, progressing along the axisymmetric expansion line, and ending in the vicinity of the point $A_{2}=A_{3}=0$, which corresponds to isotropic turbulence. In contrast, the present bubbly flows achieve a general trend that they begin from the two-component limit at the wall, traverse with a relatively short path, and then directly end very close to the axisymmetric expansion limit. Moreover, most points aggregate close to this right-hand line to axisymmetric limit for the bubble-laden DNSs, quite the opposite of what is found in the single-phase channel. The origin of this difference is due to the dominant streamwise velocity fluctuations (see figures 6, 8b, 8d and 9b) compared to the other two directions over a large range of the channel centre for the bubble-laden cases. Such a major contribution of normal stress in the streamwise direction was also found by a recent DNS study by du Cluzeau, Bois & Toutant (Reference du Cluzeau, Bois and Toutant2019), who considered this feature later in an algebraic expression for BIT.
Another important observation can be made about the locations of the points aggregated close to the line corresponding to the axisymmetric expansion limit for each particular bubble-laden case. These points are from the centre region of the channel as remarked before; the points distribution on this right-hand limit line towards the isotropic point – direction decreasing with $A_{2}$ – is in the sequence SmMany (SmFew), BiDisp, LaMany to Akbar, which corresponds to increasing bubble Reynolds number, $Re_{p}$, as illustrated in figure 12(a,b). This interesting finding can also be identified by figures 6, 8(b), 8(d), 9(b) and 10, with a gradually less dominant streamwise Reynolds normal stress for theses cases in the order mentioned above. The trend reflects that $Re_{p}$ appears to have a decisive effect on the anisotropy state of the BIT-dominated flows. Such an effect may originate from two parts: the $Re_{p}$-dependent wake structure (see the supplementary material to Santarelli & Fröhlich (Reference Santarelli and Fröhlich2016)) and the path of the rising bubbles (Horowitz & Williamson Reference Horowitz and Williamson2010; Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012).
5.2 The model and its realizability
The final step is to propose a general model for the BIT terms. Based on previous experience, this is done by providing a functional relation for the coefficient $b_{11}^{\ast }$. The important observation in § 5.1 that the magnitude of anisotropy decreases with the bubble Reynolds number suggests modelling $b_{11}^{\ast }$ as a function of $Re_{p}$. Pairs of discrete values of $b_{11}^{\ast }$ and $Re_{p}$ are displayed in figure 13. For the case BiDisp, two data points are plotted, one for the smaller bubbles and one for the larger bubbles, with the assumption that the contributions from the two bubble sizes sum up to the total contribution as discussed in § 4.4. Obviously, $b_{11}^{\ast }$ decreases with $Re_{p}$ – implying a trend towards more isotropy, which supports the observations made before. Hence, an expression for $b_{11}^{\ast }$ as a function of $Re_{p}$ seems suitable to generate a model for $\unicode[STIX]{x1D61A}_{R,ij}^{eff}$. By curve fitting (figure 13), the following effective BIT source was determined:
Here $b_{11}^{\ast }\leqslant 2$ is imposed to fulfill $b_{ii}^{\ast }=2$ and $b_{22}^{\ast }=b_{33}^{\ast }\geqslant 0$. The above expression for $b_{11}^{\ast }$ included in figure 13 provides an excellent approximation for the results over the range considered in the available data with a maximum deviation of 1.6 % at the data points. In the limit of high particle Reynolds number, the expression for $b_{11}^{\ast }$ in (5.3) approaches the value $1.34$, which implies that the Reynolds-stress anisotropy does not change much when reaching a certain $Re_{p}$ (recall that for $Re_{p}=1080$ in the case Akbar $b_{11}^{\ast }=1.4$). This is reasonable on the background of an axisymmetric wake embedded in turbulence at high enough particle Reynolds number (Bagchi & Balachandar Reference Bagchi and Balachandar2004; Rind & Castro Reference Rind and Castro2012). The former observed a ratio of the streamwise velocity fluctuation to the cross-streamwise velocity fluctuation of approximately $1.4$ for $Re_{p}=610$. That is similar to the ratio of approximately $1.3$ obtained by the latter (Rind & Castro Reference Rind and Castro2012) for an axisymmetric wake with an equivalent $Re_{p}$ (based on the wake half-width and the centreline deficit velocity) in excess of $10\,000$. Here, the relatively small change of this ratio over a large range of $Re_{p}$ supports an almost constant value for $b_{11}^{\ast }$ appearing in the source for the Reynolds-stress anisotropy at relatively large $Re_{p}$. However, the absolute values obtained by Bagchi & Balachandar (Reference Bagchi and Balachandar2004) and Rind & Castro (Reference Rind and Castro2012) cannot be directly used for correlation, since in the present study bubble swarms are considered. It might appear more direct to construct an expression for $b_{11}^{\ast }$ as a function of $A_{2}$ and $A_{3}$. Such stress invariants find many uses in SMC for single-phase flows (e.g. by Launder & Li Reference Launder and Li1994; Ristorcelli et al. Reference Ristorcelli, Lumley and Abid1995; Jakirlić & Hanjalić Reference Jakirlić and Hanjalić2002). However, this provides pitfalls for modelling $\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff}$ in bubbly flows, e.g. in cases when bubbles rise in a highly turbulent background flow. In this situation, $A_{2}$ represents both the property of the background flow as well as the property of the BIT. In contrast, the present $Re_{p}$-dependent expression for $b_{11}^{\ast }$ avoids such a pitfall.
Figure 14 shows the anisotropy-invariant map with the points derived from the EE SMC simulations. All points lie within the triangular domain, indicating that the present EE SMC employed with (5.3) provides realizable results for all the cases. Moreover, it is interesting to compare the results of the SMC with the DNS in figure 12. The assumption $b_{22}=b_{33}$, which results in $b_{22}^{\ast }=b_{33}^{\ast }$ as described in § 3.4, imposes that the results of the EE SMC yield axisymmetric turbulence in BIT-dominated cases. The DNS show that in all cases there are several points very close to the two-component limit. These points correspond to the viscous sublayer that is closest to the wall. As mentioned before, the present EE SMC in a high-$Re$ version is not intended to capture this near-wall region, hence, the two-component limit is not approached. Going further away from the wall, the model results follow the same $Re_{p}$-dependent trend as observed by DNS along the axisymmetric expansion line very well, within reasonable tolerance in the absolute values of the stress invariants.
6 Performance in the case in Lu & Tryggvason (Reference Lu and Tryggvason2008)
In this section, the results of the EE SMC using (5.3) are compared to the DNS bubbly channel data of Lu & Tryggvason (Reference Lu and Tryggvason2008). The configuration is an upward-directed channel similar to SmMany, with a smaller domain $\unicode[STIX]{x03C0}h\times 2h\times (\unicode[STIX]{x03C0}/2)h$ in streamwise $(x)$, wall-normal $(y)$ and spanwise $(z)$ directions, where $h$ is the channel half-width. Of interest in examining the proposed BIT expression (5.3), only the deformable bubble case in their study is considered, since there the Reynolds stresses are considerably larger than that in the single-phase simulation with a similar flow rate (see Lu & Tryggvason Reference Lu and Tryggvason2008). Table 4 summarizes the non-dimensional parameters used for the considered DNS case.
The number of grid points used was $56\times 60\times 51$ in the $x$-, $y$- and $z$-directions, respectively. The drag coefficient $C_{D}\approx 1.4$ was determined from the relative velocity in DNS and an optimized lift coefficient was chosen to reproduce the void fraction in DNS as well as possible, while keeping the other interfacial force modelling the same as for SmMany. This allows one to investigate how the present SMC employed with the BIT model performs in the considered case, without including compounding errors that originate from other submodels.
The results of the EE SMC are presented in figure 15 and compared with the DNS data. Indeed, it shows that the resulting $C_{L}=-0.005$ and $C_{D}$ based on the DNS data yield the correct void fraction profile (figure 15a) and both liquid and gas velocities (figure 15b). Note that complete vanishing of the gas phase close to the walls as in the DNS (see figure 15a) cannot be achieved generally using the EE approach. In figure 15(c), the streamwise and wall-normal velocity fluctuations are shown. The quantitative accuracy of the model in the channel centre is encouraging; however, it is not extremely good as in the other cases shown in § 4. To analyse the origin of this difference, attention needs to be directed to the performance of the original SSG model in the single-phase flow with a similar flow rate in figure 16 ($Re_{b}\approx 4000$, DNS from Lu & Tryggvason (Reference Lu and Tryggvason2008)). The corresponding $Re_{\unicode[STIX]{x1D70F}}=127$ is the same as in the deformable bubble case, indicating a very low turbulence level for the single-phase flow. It is clear from figure 16(a) that the high-$Re$ version of SSG fails to reproduce the sharp peak of $u^{\prime +}$ near the walls and the level of $u^{\prime +}$ is overpredicted in the channel centre. As mentioned before, this is a feature of high-$Re$ SMCs (Wilcox Reference Wilcox1998; Hanjalić & Launder Reference Hanjalić and Launder2011). This point is amplified by calculation of the present flow, with such an extremely small $Re_{\unicode[STIX]{x1D70F}}$. Checking now more closely, the slight overprediction of $u^{\prime +}$ in the channel centre for the bubble-laden case in figure 15(c) results rather from the inaccuracy of the high-$Re$ SMC performance in the corresponding single-phase flow but is not due to the BIT model. Indeed, the BIT model works very well, considering the relative attenuation for $u^{\prime +}$ in the bubble-laden case compared to that in the single-phase flow by comparing the results from the DNS and SMC, respectively.
7 Conclusions
In the present paper, a full second-moment closure (SMC) for bubble-laden flows in the framework of the Euler–Euler (EE) approach was proposed for the channel centre. Several important conclusions based on DNS data for bubbly flow were obtained that led to the development of a new bubble-induced turbulence (BIT) model at the second-moment level (5.3). The important findings are as follows.
(i) From the evaluation of the DNS data for bubbly channel flow, it was found that there are departures from the general formulation for the slow part of the pressure–strain correlation (3.17): its nonlinear term plays an important role, which requires the SSG model rather than Rotta’s linear model. Further, it confirms that the rapid part caused by mean strain in the pressure–strain correlation is negligible in the BIT-dominated region, so it does not matter what form is taken for this part.
(ii) The concept to import the definition of the effective BIT source (3.22) greatly simplifies the modelling work, so that the term $\unicode[STIX]{x1D719}_{ij,3}$ does not need to be modelled explicitly. Moreover, in the framework of this concept, the Reynolds-stress budgets have been estimated from the one-point statistics and TKE budgets using the DNS bubbly channel data.
(iii) A suitably chosen iterative procedure employing the full EE SMC provides suitable model coefficients for the closure of the terms resulting from BIT while largely removing the influence of other submodels. At the same time these results validate the closure, exhibiting very good agreement with the DNS data and better performance than the standard closures. The obtained model coefficients indicate that the isotropic assumption for the interfacial term turned out to be too simplistic and a value for $b_{11}^{\ast }$ independent of $Re_{p}$ is not general.
(iv) For the first time, an anisotropy-invariant map of bubbly flows was shown based on the available DNS data. The map identifies that $Re_{p}$ emerges as the key parameter for the anisotropy of this type flow and this is confirmed by the resulting $b_{11}^{\ast }$. A new BIT formulation for EE SMC has been developed based on this anisotropy-invariant analysis. It is worth noting that the proposal is derived from a small number of cases, but the information in each case is very sound and covers a range of $Re_{p}$ from $233$ to $1080$, which is highly relevant for practical applications.
(v) Besides the BIT-dominated flows, the present EE SMC employed with the proposed BIT model has been tested against the DNS bubbly channel data of Lu & Tryggvason (Reference Lu and Tryggvason2008). Admittedly, though the overall agreement with DNS data is encouraging, the case may not represent an optimal validation case for BIT model assessment. It has been shown that $u^{\prime +}$ in the channel centre for the bubble-laden case is only approximately twice that for the single-phase flow with a similar mass flow rate. Hence, by calculating such bubbly flows, the performance of the BIT model is contaminated by the intrinsic problem of the high-$Re$ RANS version for simulating single-phase flow. Caution should be taken when interpreting the results achieved by testing the BIT models in this type of flow. For example, in the considered bubble-laden case, using $b_{11}^{\ast }=1$ as proposed by Colombo & Fairweather (Reference Colombo and Fairweather2015) instead of the present model (5.3), a better fit of $u^{\prime +}$ profile in the channel centre can be achieved than figure 15(c). Checking closely, however, this results from underpredicting the interfacial term in the streamwise component, along with the inherent overprediction problem of the high-$Re$ RANS closure in the channel centre.
The BIT model can now be applied in EE simulations. This closure employs $Re_{p}$, which is available in any EE SMC simulation. The proposed expression can as well be combined with any $\unicode[STIX]{x1D61A}_{k}^{SMC}$ and used for any similar Reynolds-stress model. Further tests should focus on the performance of the model over a wider range of bubble Reynolds numbers and on more complex flow fields, such as flows in impeller rotation. For the latter purpose, in appendix C we have extended the BIT term for the use in multi-dimensional flows so that it is independent of the choice of the coordinate direction.
Acknowledgements
The authors would like to acknowledge Y. Liao for an inspiring discussion on (3.22) and P. Shi for enlightening suggestions for appendix C. We also gratefully acknowledge Professor A. Bragg for his helpful comments at the stage of revising the paper. The DNS data used were made available by C. Santarelli.
Appendix A. Modelling interfacial forces $\boldsymbol{M}^{K}$
In the EE model (3.2), different interfacial forces $\boldsymbol{M}^{K}$ are considered, namely drag force $\boldsymbol{F}_{D}$, lift force $\boldsymbol{F}_{L}$, wall force $\boldsymbol{F}_{W}$ and turbulent dispersion $\boldsymbol{F}_{TD}$, and given by
A.1 Drag force
The drag force results from the viscous force acting on the bubble surface and the pressure differences caused by the bubble shape. It acts as the resistance between the gas phase and the liquid phase and generates a momentum exchange due to the slip velocity between phases. The corresponding gas-phase momentum is defined as
The validity of most of the correlations for the drag coefficient $C_{D}$ is limited to restricted bubble shape regimes, e.g. Schiller & Naumann (Reference Schiller and Naumann1933) for the present bubble Reynolds number regime ($Re_{p}<800$) in contaminated air–water systems:
This model would work well in the SmFew case, where the swarm effect is negligible due to the low void fraction. However, for the denser bubble-laden cases, SmMany, LaMany and BiDisp, having more intense swarm effect, (A 3) fails to reproduce the relative velocity from DNS. To minimize the influence of the drag modelling, we use $C_{D}$ defined by (4.4) in the present study based on the relative velocity as justified by Santarelli & Fröhlich (Reference Santarelli and Fröhlich2015).
A.2 Lift force
In shear flows, bubbles experience a force perpendicular to the direction of the relative velocity. This effect generally is referred to as lift force and described by the expression (Auton Reference Auton1987)
In EE simulations, the empirical correlation of Tomiyama et al. (Reference Tomiyama, Tamai, Zun and Hosokawa2002) is frequently used. In the present study, this coefficient is optimized such that the resulting void fraction distribution matches the DNS target, while the other non-drag force models are employed as they are (see the procedure in table 2).
A.3 Wall force
The wall force behaves analogously to a lubrication force and acts on a bubble near walls to prevent bubbles from touching the wall (Antal et al. Reference Antal, Lahey and Flaherty1991). The wall force has the general form
where $C_{W}=\max (0,C_{W1}/d_{p}+C_{W2}/y)$, with $C_{W1}=-0.01$ and $C_{W2}=0.03$ in the present study; and $\hat{\boldsymbol{y}}$ is the unit normal vector perpendicular to the wall directed towards the fluid. The wall force coefficient $C_{W}$ depends on the distance $y$ to the wall and is expected to be positive, hence, the bubble is driven away from the wall.
A.4 Turbulent dispersion force
The turbulent dispersion force is the result of the turbulent fluctuations of the liquid phase. In dispersed bubbly flows, turbulence in the continuous liquid phase causes bubbles to be transferred from regions of high concentration to regions of low concentration. In the present study, the model of Burns et al. (Reference Burns, Frank, Hamill and Shi2004) derived by Favre averaging of the drag force is used,
with $\unicode[STIX]{x1D70E}_{TD}$ the turbulent Schmidt number, typically taken to be $0.9$.
Appendix B. The present Reynolds-stress model
The transport equations for the Reynolds stresses read
with $c_{s}=1.63$. The modelled modified pressure–strain term $\unicode[STIX]{x1D719}_{ij}^{SMC\text{-}mod}$ is identical to the SSG model (two-phase version) and can be expressed in terms of the stress production tensor $\unicode[STIX]{x1D617}_{ij}$, its complement $\unicode[STIX]{x1D60C}_{ij}$ and the mean rate of strain $\unicode[STIX]{x1D60D}_{ij}$:
where
The constants in (B 2) are (Hanjalić & Launder Reference Hanjalić and Launder2011):
The effective BIT source term $\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff}$ is
while $\unicode[STIX]{x1D61A}_{k}^{SMC}$, the interfacial term for the $k$ equation, is adopted from Ma et al. (Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017) as
The turbulent dissipation rate, $\unicode[STIX]{x1D700}$, in (
B 1) is obtained from its own transport equation:
where $C_{\unicode[STIX]{x1D700}1}=1.45$, $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D700}}=1.36$ and $C_{\unicode[STIX]{x1D700}2}=1.83$. Here $\unicode[STIX]{x1D70F}$ is a time scale, reading $\unicode[STIX]{x1D70F}=d_{p}/u_{r}$.
Appendix C. Extending the BIT source for multi-dimensional flows
The present effective BIT source term (B 5) is applicable for the case where the relative velocity is aligned with the vertical direction, such as the present DNS channel, vertical bubble columns as well simple pipe flows. In multi-dimensional flows, e.g. flows in stirred tanks (Shi & Rzehak Reference Shi and Rzehak2018), however, no such simple correspondence is valid. Here, we use the idea based on proper Euler angles (figure 17) to cast the BIT term so that it is independent of the choice of coordinate directions. We define a fixed coordinate system $xyz$ and a rotating (body-fixed) coordinate system $XYZ$, in which $X$ is aligned with the local relative velocity between two phases. Any body-fixed coordinate system $XYZ$ can be achieved by composing three elemental rotations, starting from an original fixed coordinate system $xyz$. Here, we use a common sequence of rotations – the so-called $x$–$y^{\prime }$–$x^{\prime \prime }$ rotation, based on the orientations denoted as follows:
(i) $x$–$y$–$z$ (initial);
(ii) $x^{\prime }$–$y^{\prime }$–$z^{\prime }$ (after first rotation around the $x$-axis with the angle $\unicode[STIX]{x1D6FC}$);
(iii) $x^{\prime \prime }$–$y^{\prime \prime }$–$z^{\prime \prime }$ (after second rotation around the $y^{\prime }$-axis with the angle $\unicode[STIX]{x1D6FD}$);
(iv) $X$–$Y$–$Z$ (after third rotation around the $x^{\prime \prime }$-axis with the angle $\unicode[STIX]{x1D6FE}$).
Since the third rotation occurs about $X$, it does not change the orientation of $X$. Hence $X$ coincides with $x^{\prime \prime }$.
The present BIT term takes into account that the source in parallel with and perpendicular to the direction of the relative velocity is different. We write it in the currently defined $XYZ$-frame as a tensor:
The term is aligned with the values in its parallel direction,
and in its perpendicular direction,
where $b_{11}^{\ast }$ is from (B 5).
The task now is to transform $(\unicode[STIX]{x1D64E}_{R}^{SMC\text{-}eff})_{XYZ}$ in the fixed $xyz$-frame using the $x^{\prime \prime }$–$y^{\prime }$–$x$ rotation (the reverse one compared to the $x$–$y^{\prime }$–$x^{\prime \prime }$ rotation mentioned before). This can be achieved by
where $\unicode[STIX]{x1D64C}$ is a rotation matrix and can be decomposed as the product of the latter two elemental rotation matrices ($y^{\prime }$–$x$ rotations):
The first rotation $x^{\prime \prime }$ with the related Euler angle $\unicode[STIX]{x1D6FE}$ indeed is not needed in the present study, since we are only interested in the orientation of $X$ (the direction of the local relative velocity).
The final step is to find the Euler angles $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$. This can be achieved based on the known unit relative velocity vectors $(1,0,0)$ in the $XYZ$-frame and $(\hat{u} _{1},\hat{u} _{2},\hat{u} _{3})$ in the $xyz$-frame ($\hat{\boldsymbol{u}}_{r}=\boldsymbol{u}_{r}/|\boldsymbol{u}_{r}|$). After some projections (shown in figure 17) and linear algebra, the results are