Hostname: page-component-78c5997874-94fs2 Total loading time: 0 Render date: 2024-11-08T05:34:59.448Z Has data issue: false hasContentIssue false

Progress in the second-moment closure for bubbly flow based on direct numerical simulation data

Published online by Cambridge University Press:  25 November 2019

Tian Ma*
Affiliation:
Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden–Rossendorf, 01328 Dresden, Germany
Dirk Lucas
Affiliation:
Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden–Rossendorf, 01328 Dresden, Germany
Suad Jakirlić
Affiliation:
Institute for Fluid Mechanics and Aerodynamics, Technische Universität Darmstadt, 64287 Darmstadt, Germany
Jochen Fröhlich
Affiliation:
Institute of Fluid Mechanics, Technische Universität Dresden, 01062 Dresden, Germany
*
Email address for correspondence: [email protected]

Abstract

Data from direct numerical simulations (DNS) of disperse bubbly flow in an upward vertical channel are used to develop a new second-moment closure for bubble-induced turbulence (BIT) in the Euler–Euler framework. The closure is an extension of a BIT model originally proposed by Ma et al. (Phys. Rev. Fluids, vol. 2, 2017, 034301) for two-equation eddy-viscosity models and focuses on the core region of the channel, where the interfacial term and dissipation term are in balance. Particular attention in this study is given to the treatment of the pressure–strain term for bubbly flows and the form of the interfacial term to account for BIT. For the latter, the concept of an effective BIT source is proposed, which leads to a significant simplification of the modelling work for both the pressure–strain correlation and the interfacial term itself. The anisotropy of bubbly flow is analysed with the aid of the anisotropy-invariant map obtained from the DNS data, and a parameter governing this issue is established. The complete second-moment closure is tested against reference data for different bubbly channel flows and a case of a bubble column. The agreement achieved with the DNS data is very good and the performance of the new model is better than obtained with the standard procedure. Furthermore, the model is shown to be robust and to fulfil the requirements of realizability.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2019. Published by Cambridge University Press

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.

Figure 1. Instantaneous DNS data for the case SmMany of Santarelli & Fröhlich (Reference Santarelli and Fröhlich2015). The bubble size is to scale and the vertical plane with the contour plot shows the instantaneous streamwise liquid velocity. (Main part of the picture reprinted from Santarelli (Reference Santarelli2015) with permission from TUDpress.)

Table 1. Parameters of the cases used for the present study according to Santarelli & Fröhlich (Reference Santarelli and Fröhlich2015, Reference Santarelli and Fröhlich2016). The labels Sm and La used for the case with the bidisperse bubble swarm, BiDisp, indicate averaging over the small and large bubbles of the case BiDisp, respectively. Here, $N_{p}$ is the number of bubbles, $\unicode[STIX]{x1D6FC}$ the void fraction, $d_{p}$ the particle diameter and $Ar$ the Archimedes number. The values of $Re_{\unicode[STIX]{x1D70F}}$, the friction Reynolds number, $Re_{p}$, the particle Reynolds number based on $d_{p}$ and the relative velocity, and $C_{D}$, the drag coefficient obtained according to (4.4) below, are results of the simulations.

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

(2.1)$$\begin{eqnarray}Ar=\frac{|\unicode[STIX]{x1D70C}^{G}-\unicode[STIX]{x1D70C}^{L}|gd_{p}^{3}}{\unicode[STIX]{x1D70C}^{L}\unicode[STIX]{x1D708}^{2}},\end{eqnarray}$$

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. Distribution of $C_{\unicode[STIX]{x1D707}}$ in a $k$$\unicode[STIX]{x1D700}$ type model computed using DNS data.

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)

(2.2)$$\begin{eqnarray}C_{\unicode[STIX]{x1D707}}=\frac{\unicode[STIX]{x1D708}_{t}\unicode[STIX]{x1D700}}{k^{2}}\end{eqnarray}$$

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.

Figure 3. Shear stress and shear-stress-induced production of TKE for the DNS. (a) Distribution of the structure parameter, $-\overline{\overline{u^{\prime }v^{\prime }}}/k$. (b) Distribution of the ratio of production by mean flow deformation to dissipation of TKE, $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}$.

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

(2.3)$$\begin{eqnarray}C_{\unicode[STIX]{x1D707}}=\frac{(\overline{\overline{u^{\prime }v^{\prime }}}/k)^{2}}{\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}}.\end{eqnarray}$$

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

(3.1)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D6FC}^{K}\unicode[STIX]{x1D70C}^{K})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}^{K}\unicode[STIX]{x1D70C}^{K}\overline{\overline{\boldsymbol{u}}}^{K})=0, & \displaystyle\end{eqnarray}$$
(3.2)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}(\unicode[STIX]{x1D6FC}^{K}\unicode[STIX]{x1D70C}^{K}\overline{\overline{\boldsymbol{u}}}^{K})}{\text{D}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(2\unicode[STIX]{x1D6FC}^{K}\unicode[STIX]{x1D707}^{K}\overline{\overline{\unicode[STIX]{x1D64E}}}^{K})-\unicode[STIX]{x1D6FC}^{K}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FC}^{K}\unicode[STIX]{x1D70C}^{K}\boldsymbol{g}+\boldsymbol{M}^{K}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}^{K}\overline{\overline{\unicode[STIX]{x1D749}}}_{t}^{K}), & \displaystyle\end{eqnarray}$$

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,

(3.3)$$\begin{eqnarray}\frac{\text{D}}{\text{D}t}(\overline{\unicode[STIX]{x1D711}}\overline{\overline{u_{i}^{\prime }u_{j}^{\prime }}})=\unicode[STIX]{x1D617}_{ij}+\unicode[STIX]{x1D60B}_{ij}+\unicode[STIX]{x1D719}_{ij}+\unicode[STIX]{x1D700}_{ij}+\unicode[STIX]{x1D61A}_{R,ij},\end{eqnarray}$$

where the terms on the right-hand side of (3.3) read

(3.4)$$\begin{eqnarray}\displaystyle & \displaystyle \text{production}\quad \unicode[STIX]{x1D617}_{ij}=-\overline{\unicode[STIX]{x1D711}}\overline{\overline{u_{i}^{\prime }u_{k}^{\prime }}}\frac{\unicode[STIX]{x2202}\overline{\overline{u_{j}}}}{\unicode[STIX]{x2202}x_{k}}-\overline{\unicode[STIX]{x1D711}}\overline{\overline{u_{j}^{\prime }u_{k}^{\prime }}}\frac{\unicode[STIX]{x2202}\overline{\overline{u_{i}}}}{\unicode[STIX]{x2202}x_{k}}, & \displaystyle\end{eqnarray}$$
(3.5)$$\begin{eqnarray}\displaystyle & \displaystyle \text{diffusion}\quad \unicode[STIX]{x1D60B}_{ij}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{k}}\left(\frac{1}{\unicode[STIX]{x1D70C}}\overline{\unicode[STIX]{x1D711}}\overline{\overline{p^{\prime }(\unicode[STIX]{x1D6FF}_{jk}u_{i}^{\prime }+\unicode[STIX]{x1D6FF}_{ik}u_{j}^{\prime })}}+\overline{\unicode[STIX]{x1D711}}\overline{\overline{u_{i}^{\prime }u_{j}^{\prime }u_{k}^{\prime }}}-\frac{1}{\unicode[STIX]{x1D70C}}\overline{\unicode[STIX]{x1D711}}(\overline{\overline{u_{j}^{\prime }\unicode[STIX]{x1D70F}_{ik}^{\prime }}}+\overline{\overline{u_{i}^{\prime }\unicode[STIX]{x1D70F}_{jk}^{\prime }}})\right), & \displaystyle\end{eqnarray}$$
(3.6)$$\begin{eqnarray}\displaystyle & \displaystyle \text{pressure{-}strain}\quad \unicode[STIX]{x1D719}_{ij}=\frac{1}{\unicode[STIX]{x1D70C}}\overline{\unicode[STIX]{x1D711}}\overline{\overline{p^{\prime }\left(\frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}^{\prime }}{\unicode[STIX]{x2202}x_{i}}\right)}}, & \displaystyle\end{eqnarray}$$
(3.7)$$\begin{eqnarray}\displaystyle & \displaystyle \text{dissipation}\quad \unicode[STIX]{x1D700}_{ij}=-\frac{1}{\unicode[STIX]{x1D70C}}\overline{\unicode[STIX]{x1D711}}\overline{\overline{\unicode[STIX]{x1D70F}_{ik}^{\prime }\left(\frac{\unicode[STIX]{x2202}u_{j}^{\prime }}{\unicode[STIX]{x2202}x_{k}}\right)}}-\frac{1}{\unicode[STIX]{x1D70C}}\overline{\unicode[STIX]{x1D711}}\overline{\overline{\unicode[STIX]{x1D70F}_{jk}^{\prime }\left(\frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{k}}\right)}}. & \displaystyle\end{eqnarray}$$

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

(3.8)$$\begin{eqnarray}\unicode[STIX]{x1D61A}_{R,ij}=-\frac{1}{\unicode[STIX]{x1D70C}}(\overline{p_{L}^{\prime }u_{L,j}^{\prime }n_{i}I}+\overline{p_{L}^{\prime }u_{L,i}^{\prime }n_{j}I})+\frac{1}{\unicode[STIX]{x1D70C}}(\overline{\unicode[STIX]{x1D70F}_{L,ik}^{\prime }u_{L,j}^{\prime }n_{k}I}+\overline{\unicode[STIX]{x1D70F}_{L,jk}^{\prime }u_{L,i}^{\prime }n_{k}I})\end{eqnarray}$$

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.9)$$\begin{eqnarray}\unicode[STIX]{x2202}\unicode[STIX]{x1D711}/\unicode[STIX]{x2202}x_{i}=-In_{i}.\end{eqnarray}$$

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,

(3.10)$$\begin{eqnarray}\frac{\text{D}(\unicode[STIX]{x1D6FC}^{L}\overline{\overline{u_{i}^{\prime }u_{j}^{\prime }}})}{\text{D}t}=\unicode[STIX]{x1D617}_{ij}^{SMC}+\unicode[STIX]{x1D60B}_{ij}^{SMC}+\unicode[STIX]{x1D719}_{ij}^{SMC}+\unicode[STIX]{x1D700}_{ij}^{SMC}+\unicode[STIX]{x1D61A}_{R,ij}^{SMC},\end{eqnarray}$$

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

(3.11)$$\begin{eqnarray}\unicode[STIX]{x1D60B}_{ij}^{SMC}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{k}}\left(\unicode[STIX]{x1D6FC}^{L}(\unicode[STIX]{x1D708}^{L}+c_{s}\unicode[STIX]{x1D708}_{t})\frac{\unicode[STIX]{x2202}\overline{\overline{u_{i}^{\prime }u_{j}^{\prime }}}}{\unicode[STIX]{x2202}x_{k}}\right),\end{eqnarray}$$

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

(3.12)$$\begin{eqnarray}\unicode[STIX]{x1D700}_{ij}^{SMC}=-\unicode[STIX]{x1D6FC}^{L}{\textstyle \frac{2}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D700}.\end{eqnarray}$$

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,

(3.13)

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

(3.14)$$\begin{eqnarray}\unicode[STIX]{x1D61A}_{k}^{SMC}=\min (0.18Re_{p}^{0.23},1)\boldsymbol{F}_{D}(\boldsymbol{u}^{G}-\boldsymbol{u}^{L}),\end{eqnarray}$$

with $\boldsymbol{F}_{D}$ the drag and

(3.15)$$\begin{eqnarray}\unicode[STIX]{x1D70F}=\frac{d_{p}}{u_{r}}\end{eqnarray}$$

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

(3.16)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{ij}=\unicode[STIX]{x1D719}_{ij,1}+\unicode[STIX]{x1D719}_{ij,2}+\unicode[STIX]{x1D719}_{ij,3}+\unicode[STIX]{x1D719}_{ij,w},\end{eqnarray}$$

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)

(3.17)

where $\unicode[STIX]{x1D622}_{ij}$ is the Reynolds-stress anisotropy tensor

(3.18)$$\begin{eqnarray}\unicode[STIX]{x1D622}_{ij}=\frac{\overline{\overline{u_{i}^{\prime }u_{j}^{\prime }}}}{k}-\frac{2}{3}\unicode[STIX]{x1D6FF}_{ij}\end{eqnarray}$$

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.

Figure 4. Comparison between the interfacial term $S_{k}=\frac{1}{2}\unicode[STIX]{x1D61A}_{R,ii}$ according to (3.8) and the a priori evaluation of the pressure–strain terms split into a linear slow term $\unicode[STIX]{x1D719}_{ij,1,linear}$, a nonlinear slow term $\unicode[STIX]{x1D719}_{ij,1,nonlinear}$ and a rapid term $\unicode[STIX]{x1D719}_{ij,2}$ for SmMany, all normalized with $U_{b}^{3}/H$: (a$\unicode[STIX]{x1D719}_{11}$ component, (b$\unicode[STIX]{x1D719}_{22}$ component, (c$\unicode[STIX]{x1D719}_{33}$ component and (d$\unicode[STIX]{x1D719}_{12}$ component.

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),

(3.19)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{ij,2}^{SMC}=-c_{2}(\unicode[STIX]{x1D617}_{ij}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D617}_{kk}),\end{eqnarray}$$

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 4ac), 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,

(3.20)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{ij,3}^{SMC}=-c_{3}(\unicode[STIX]{x1D61A}_{R,ij}^{SMC}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D61A}_{R,kk}^{SMC}),\end{eqnarray}$$

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.

(3.21)

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.

(3.22)

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

(3.23)$$\begin{eqnarray}b_{22}^{\ast }=b_{33}^{\ast }={\textstyle \frac{1}{2}}(2-b_{11}^{\ast }).\end{eqnarray}$$

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

(3.24)$$\begin{eqnarray}\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff}=\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff(Sm)}+\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff(La)}\end{eqnarray}$$

for the effective BIT source terms in the Reynolds-stress equations and

(3.25)$$\begin{eqnarray}\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D700}}^{SMC}=\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D700}}^{SMC(Sm)}+\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D700}}^{SMC(La)}\end{eqnarray}$$

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

(4.1)

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,

(4.2)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{ij}^{mod} & = & \displaystyle -c_{1}\unicode[STIX]{x1D6FC}^{L,DNS}\unicode[STIX]{x1D700}^{DNS}\unicode[STIX]{x1D622}_{ij}^{DNS}+c_{1}^{\prime }\unicode[STIX]{x1D6FC}^{L,DNS}\unicode[STIX]{x1D700}^{DNS}(\unicode[STIX]{x1D622}_{ik}^{DNS}\unicode[STIX]{x1D622}_{kj}^{DNS}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}A_{2}^{DNS})\nonumber\\ \displaystyle & & \displaystyle -\,c_{2}(\unicode[STIX]{x1D617}_{ij}^{DNS}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D617}_{kk}^{DNS}),\end{eqnarray}$$

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

(4.3)$$\begin{eqnarray}\unicode[STIX]{x1D700}_{ij}=-\unicode[STIX]{x1D6FC}^{L,DNS}{\textstyle \frac{2}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D700}^{DNS}.\end{eqnarray}$$

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

(4.4)$$\begin{eqnarray}C_{D}^{DNS}=\left(\frac{4}{3}\unicode[STIX]{x1D6FC}^{L}\frac{d_{p}(\unicode[STIX]{x1D70C}^{L}-\unicode[STIX]{x1D70C}^{G})g}{\unicode[STIX]{x1D70C}^{L}u_{r}^{2}}\right)^{DNS}.\end{eqnarray}$$

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}}$.

Table 2. Optimization procedure for the evaluation of $b_{11}^{\ast }$.

4.2 Closure for the case SmMany

Figure 5. Mean flow statistics from the present EE SMC and DNS data for the case SmMany: (a) gas void fraction; (b) liquid streamwise velocity and gas streamwise velocity.

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.

(4.5)$$\begin{eqnarray}b_{11}^{C}=b_{22}^{C}=b_{33}^{C}=2/3.\end{eqnarray}$$

Colombo & Fairweather (Reference Colombo and Fairweather2015) suggested

(4.6a,b)$$\begin{eqnarray}b_{11}^{CF}=1,\quad b_{22}^{CF}=b_{33}^{CF}=0.5.\end{eqnarray}$$

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.

Figure 6. Reynolds-stress components from the present EE SMC and DNS data for the case SmMany: (a) streamwise Reynolds normal stress $\overline{\overline{u^{\prime }u^{\prime }}}$; (b) wall-normal Reynolds normal stress $\overline{\overline{v^{\prime }v^{\prime }}}$; (c) Reynolds shear stress $\overline{\overline{u^{\prime }v^{\prime }}}$. All components are normalized by $U_{b}^{2}$. The label Cokljat refers to EE SMC with $b_{ij}^{C}$ from (4.5), while the label Colombo refers to EE SMC with $b_{ij}^{CF}$ from (4.6) as detailed in the text.

Figure 7. Comparison of key contributions in the present EE SMC with the DNS data for the case SmMany as described in the text. (a) Selected terms in the budget of $\overline{\overline{u^{\prime }u^{\prime }}}$ with all terms normalized by $U_{b}^{3}/H$. (b) Selected terms in the budget of $\overline{\overline{v^{\prime }v^{\prime }}}$ with all terms normalized by $U_{b}^{3}/H$. (c) Selected terms in the budget of $\unicode[STIX]{x1D700}$ with all terms normalized by $U_{b}^{4}/H^{2}$.

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

(4.7)$$\begin{eqnarray}0\approx \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D700}}+S_{\unicode[STIX]{x1D700}}\end{eqnarray}$$

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

(4.8)$$\begin{eqnarray}\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D700}}=-\unicode[STIX]{x1D6FC}^{L,DNS}C_{\unicode[STIX]{x1D700}2}\frac{(\unicode[STIX]{x1D700}^{DNS})^{2}}{k^{DNS}},\end{eqnarray}$$

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

Figure 8. Gas void fraction and Reynolds-stress components from the present EE SMC and DNS data for the other two monodisperse cases: (a,b) SmFew; and (c,d) LaMany.

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).

Figure 9. (a) Gas void fraction and (b) Reynolds-stress components from the present EE SMC and DNS data for the case BiDisp.

Table 3. Values of $b_{11}^{\ast }$ and $C_{L}$ obtained by the iterative process for all cases.

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

Figure 10. Comparison of the present SMC with experimental data of Akbar et al. (Reference Akbar, Hayashi, Hosokawa and Tomiyama2012): (a) gas void fraction, (b) liquid streamwise velocity and gas streamwise velocity, (c) liquid Reynolds normal stresses and (d) liquid Reynolds shear stress.

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),

(5.1)$$\begin{eqnarray}A=1-9(A_{2}-A_{3})/8.\end{eqnarray}$$

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$.

Figure 11. Anisotropy invariants: (a) the single-phase case Unladen; and (b) the case SmMany.

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:

(5.2)$$\begin{eqnarray}6^{1/3}{|A_{3}|}^{2/3}\leqslant A_{2}\leqslant {\textstyle \frac{8}{9}}+A_{3}.\end{eqnarray}$$

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. Anisotropy-invariant map of the DNS data for all cases considered. (a) The cases SmMany, LaMany, BiDisp and Akbar; $Re_{p}$ of the case BiDisp is calculated from the averaged $u_{r}$ and the averaged $d_{p}$ from the DNS of BiDisp. (b) The case SmFew and (c) the single-phase case Unladen.

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 68b8d 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

Figure 13. Values of $b_{11}^{\ast }$ as a function of $Re_{p}$ for all cases considered. The dashed line shows the fit.

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:

(5.3)

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. Anisotropy-invariant map of the present SMC for all bubble-laden cases considered.

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.

Table 4. DNS parameters of the deformable bubble case in Lu & Tryggvason (Reference Lu and Tryggvason2008).

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.

Figure 15. Comparison of the present SMC with DNS data for the deformable bubble case in Lu & Tryggvason (Reference Lu and Tryggvason2008): (a) gas void fraction; (b) liquid streamwise velocity and gas streamwise velocity; (c) liquid streamwise and wall-normal velocity fluctuations normalized by $u_{\unicode[STIX]{x1D70F}}$; and (d) liquid Reynolds shear stress normalized by $u_{\unicode[STIX]{x1D70F}}^{2}$.

Figure 16. Comparison of the SSG with DNS data for the single-phase case in Lu & Tryggvason (Reference Lu and Tryggvason2008): (a) streamwise and wall-normal velocity fluctuations normalized by $u_{\unicode[STIX]{x1D70F}}$; and (b) Reynolds shear stress normalized by $u_{\unicode[STIX]{x1D70F}}^{2}$.

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.

  1. (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.

  2. (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.

  3. (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.

  4. (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.

  5. (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)$$\begin{eqnarray}\boldsymbol{M}^{G}=-\boldsymbol{M}^{L}=\boldsymbol{F}_{D}+\boldsymbol{F}_{L}+\boldsymbol{F}_{W}+\boldsymbol{F}_{TD}.\end{eqnarray}$$

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

(A 2)$$\begin{eqnarray}\boldsymbol{F}_{D}=\frac{3}{4d_{p}}C_{D}\unicode[STIX]{x1D70C}^{L}\unicode[STIX]{x1D6FC}^{G}|\boldsymbol{u}^{G}-\boldsymbol{u}^{L}|(\boldsymbol{u}^{G}-\boldsymbol{u}^{L}).\end{eqnarray}$$

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:

(A 3)$$\begin{eqnarray}C_{D}=\frac{24}{Re_{p}}(1+0.15Re_{p}^{0.678}).\end{eqnarray}$$

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)

(A 4)$$\begin{eqnarray}\boldsymbol{F}_{L}=-C_{L}\unicode[STIX]{x1D70C}^{L}\unicode[STIX]{x1D6FC}^{G}(\boldsymbol{u}^{G}-\boldsymbol{u}^{L})\times \text{rot}(\boldsymbol{u}^{L}).\end{eqnarray}$$

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

(A 5)$$\begin{eqnarray}\boldsymbol{F}_{W}=\frac{2}{d_{p}}C_{W}\unicode[STIX]{x1D70C}^{L}\unicode[STIX]{x1D6FC}^{G}|\boldsymbol{u}^{G}-\boldsymbol{u}^{L}|^{2}\hat{\boldsymbol{y}},\end{eqnarray}$$

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,

(A 6)$$\begin{eqnarray}\boldsymbol{F}_{D}=-\frac{3}{4d_{p}}C_{D}\unicode[STIX]{x1D70C}^{L}\unicode[STIX]{x1D6FC}^{G}|\boldsymbol{u}^{G}-\boldsymbol{u}^{L}|\frac{\unicode[STIX]{x1D708}_{t}}{\unicode[STIX]{x1D70E}_{TD}}\left(\frac{1}{\unicode[STIX]{x1D6FC}^{L}}+\frac{1}{\unicode[STIX]{x1D6FC}^{G}}\right)\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}^{G},\end{eqnarray}$$

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

(B 1)

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}$:

(B 2)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{ij}^{SMC\text{-}mod} & = & \displaystyle -c_{1}\unicode[STIX]{x1D6FC}^{L}\unicode[STIX]{x1D700}\unicode[STIX]{x1D622}_{ij}+c_{1}^{\prime }\unicode[STIX]{x1D6FC}^{L}\unicode[STIX]{x1D700}(\unicode[STIX]{x1D622}_{ik}\unicode[STIX]{x1D622}_{kj}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}A_{2})\nonumber\\ \displaystyle & & \displaystyle -\,c_{2}^{\ast }(\unicode[STIX]{x1D617}_{ij}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D617}_{kk})-c_{3}^{\ast }(\unicode[STIX]{x1D60C}_{ij}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D60C}_{kk})-c_{4}^{\ast }\unicode[STIX]{x1D6FC}^{L}k\unicode[STIX]{x1D60D}_{ij}-c_{5}^{\ast }\unicode[STIX]{x1D622}_{ij}\unicode[STIX]{x1D617}_{kk},\end{eqnarray}$$

where

(B 3a,b)$$\begin{eqnarray}\unicode[STIX]{x1D60C}_{ij}=-\unicode[STIX]{x1D6FC}^{L}\left(\overline{\overline{u_{i}^{\prime }u_{k}^{\prime }}}\frac{\unicode[STIX]{x2202}\overline{\overline{u_{k}}}}{\unicode[STIX]{x2202}x_{j}}+\overline{\overline{u_{j}^{\prime }u_{k}^{\prime }}}\frac{\unicode[STIX]{x2202}\overline{\overline{u_{k}}}}{\unicode[STIX]{x2202}x_{i}}\right),\quad \unicode[STIX]{x1D60D}_{ij}=\frac{1}{2}\left(\frac{\unicode[STIX]{x2202}\overline{\overline{u_{i}}}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}\overline{\overline{u_{j}}}}{\unicode[STIX]{x2202}x_{i}}\right).\end{eqnarray}$$

The constants in (B 2) are (Hanjalić & Launder Reference Hanjalić and Launder2011):

(B 4)$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}c_{1}=1.7,\quad c_{1}^{\prime }=1.05,\quad c_{2}^{\ast }=0.4125,\quad c_{3}^{\ast }=0.2125,\\ c_{4}^{\ast }=0.033+0.65A_{2}^{1/2},\quad c_{5}^{\ast }=0.45.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The effective BIT source term $\unicode[STIX]{x1D61A}_{R,ij}^{SMC\text{-}eff}$ is

(B 5)

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

(B 6)$$\begin{eqnarray}\unicode[STIX]{x1D61A}_{k}^{SMC}=\min (0.18Re_{p}^{0.23},1)\boldsymbol{F}_{D}(\boldsymbol{u}^{G}-\boldsymbol{u}^{L}).\end{eqnarray}$$

The turbulent dissipation rate, $\unicode[STIX]{x1D700}$, in (

B 1

) is obtained from its own transport equation:

(B 7)

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

Figure 17. Geometrical definition of proper Euler angles (standard $y$-convention). The $xyz$ (fixed) system is shown in blue; the $XYZ$ (rotated) system aligned with the local relative velocity in $X$-axis is shown in red. The line of nodes ($N$) as the intersection of the planes $yz$ and $YZ$ is shown in green.

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:

  1. (i) $x$$y$$z$ (initial);

  2. (ii) $x^{\prime }$$y^{\prime }$$z^{\prime }$ (after first rotation around the $x$-axis with the angle $\unicode[STIX]{x1D6FC}$);

  3. (iii) $x^{\prime \prime }$$y^{\prime \prime }$$z^{\prime \prime }$ (after second rotation around the $y^{\prime }$-axis with the angle $\unicode[STIX]{x1D6FD}$);

  4. (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:

(C 1)$$\begin{eqnarray}(\unicode[STIX]{x1D64E}_{R}^{SMC\text{-}eff})_{XYZ}=\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D61A}_{\Vert }^{SMC\text{-}eff} & 0 & 0\\ 0 & \unicode[STIX]{x1D61A}_{\bot }^{SMC\text{-}eff} & 0\\ 0 & 0 & \unicode[STIX]{x1D61A}_{\bot }^{SMC\text{-}eff}\end{array}\right)_{XYZ}.\end{eqnarray}$$

The term is aligned with the values in its parallel direction,

(C 2)$$\begin{eqnarray}\unicode[STIX]{x1D61A}_{\Vert }^{SMC\text{-}eff}=b_{11}^{\ast }\unicode[STIX]{x1D61A}_{k}^{SMC},\end{eqnarray}$$

and in its perpendicular direction,

(C 3)$$\begin{eqnarray}\unicode[STIX]{x1D61A}_{\bot }^{SMC\text{-}eff}={\textstyle \frac{1}{2}}(2-b_{11}^{\ast })\unicode[STIX]{x1D61A}_{k}^{SMC},\end{eqnarray}$$

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

(C 4)$$\begin{eqnarray}(\unicode[STIX]{x1D64E}_{R}^{SMC\text{-}eff})_{xyz}=\unicode[STIX]{x1D64C}\boldsymbol{\cdot }(\unicode[STIX]{x1D64E}_{R}^{SMC\text{-}eff})_{XYZ}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}^{\text{T}},\end{eqnarray}$$

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):

(C 5)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64C} & = & \displaystyle \left(\begin{array}{@{}ccc@{}}1 & 0 & 0\\ 0 & \cos \unicode[STIX]{x1D6FC} & -\!\sin \unicode[STIX]{x1D6FC}\\ 0 & \sin \unicode[STIX]{x1D6FC} & \cos \unicode[STIX]{x1D6FC}\end{array}\right)\left(\begin{array}{@{}ccc@{}}\cos \unicode[STIX]{x1D6FD} & 0 & \sin \unicode[STIX]{x1D6FD}\\ 0 & 1 & 0\\ -\!\sin \unicode[STIX]{x1D6FD} & 0 & \cos \unicode[STIX]{x1D6FD}\end{array}\right)\nonumber\\ \displaystyle & = & \displaystyle \left(\begin{array}{@{}ccc@{}}\cos \unicode[STIX]{x1D6FD} & 0 & \sin \unicode[STIX]{x1D6FD}\\ \sin \unicode[STIX]{x1D6FC}\sin \unicode[STIX]{x1D6FD} & \cos \unicode[STIX]{x1D6FC} & -\!\sin \unicode[STIX]{x1D6FC}\cos \unicode[STIX]{x1D6FD}\\ -\!\cos \unicode[STIX]{x1D6FC}\sin \unicode[STIX]{x1D6FD} & \sin \unicode[STIX]{x1D6FC} & \cos \unicode[STIX]{x1D6FC}\cos \unicode[STIX]{x1D6FD}\end{array}\right).\end{eqnarray}$$

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

(C 6a,b)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\arccos \left(-\hat{u} _{3}/\sqrt{1-{\hat{u} _{1}}^{2}}\right),\quad \unicode[STIX]{x1D6FD}=\arccos (\hat{u} _{1}).\end{eqnarray}$$

Footnotes

Present address: Department of Civil and Environmental Engineering, Duke University, Durham, NC 27708, USA.

References

Akbar, M. H. M., Hayashi, K., Hosokawa, S. & Tomiyama, A. 2012 Bubble tracking simulation of bubble-induced pseudo turbulence. Multiphase Sci. Technol. 24, 197222.CrossRefGoogle Scholar
Antal, S. P., Lahey, R. T. & Flaherty, J. E. 1991 Analysis of phase distribution in fully developed laminar bubbly two-phase flow. Intl J. Multiphase Flow 17, 635.CrossRefGoogle Scholar
Auton, T. R. 1987 The lift force on a spherical body in a rotational flow. J. Fluid Mech. 183, 199218.CrossRefGoogle Scholar
Bagchi, P. & Balachandar, S. 2004 Response of the wake of an isolated particle to an isotropic turbulent flow. J. Fluid Mech. 518, 95123.CrossRefGoogle Scholar
Balachandar, S. & Eaton, J. K. 2010 Turbulent dispersed multiphase flow. Annu. Rev. Fluid Mech. 42, 111133.CrossRefGoogle Scholar
Barth, T. & Jesperson, D. 1989 The design and application of upwind schemes on unstructured meshes. AIAA Paper 89, 0366.Google Scholar
Biesheuvel, A. & van Wijngaarden, L. 1984 Two-phase flow equations for a dilute dispersion of gas bubbles in liquid. J. Fluid Mech. 148, 301318.CrossRefGoogle Scholar
Bolotnov, I., Jansen, K., Drew, D., Oberai, A., Lahey, R. & Podowski, M. 2011 Detached direct numerical simulations of turbulent two-phase bubbly channel flow. Intl J. Multiphase Flow 37 (6), 647659.CrossRefGoogle Scholar
Burns, A., Frank, T., Hamill, I. & Shi, J.-M. 2004 The Favre averaged drag model for turbulent dispersion in Eulerian multi-phase flows. In 5th International Conference on Multiphase Flow, Yokohama, Japan, Paper No. 392. ICMF.Google Scholar
Cifani, P., Kuerten, J. G. M. & Geurts, B. J. 2018 Highly scalable DNS solver for turbulent bubble-laden channel flow. Comput. Fluids 172, 6783.CrossRefGoogle Scholar
du Cluzeau, A., Bois, G. & Toutant, A. 2019 Analysis and modelling of Reynolds stresses in turbulent bubbly up-flows from direct numerical simulations. J. Fluid Mech. 866, 132168.CrossRefGoogle Scholar
Cokljat, D., Slack, M., Vasquez, S. A., Bakker, A. & Montante, G. 2006 Reynolds-stress model for Eulerian multiphase. Prog. Comput. Fluid Dyn. Intl J. 6, 168178.CrossRefGoogle Scholar
Colombo, M. & Fairweather, M. 2015 Multiphase turbulence in bubbly flows: RANS simulations. Intl J. Multiphase Flow 77, 222243.CrossRefGoogle Scholar
Dabiri, S., Lu, J. & Tryggvason, G. 2013 Transition between regimes of a vertical channel bubbly upflow due to bubble deformability. Phys. Fluids 25, 102110.CrossRefGoogle Scholar
Deen, N. G., Solberg, T. & Hjertager, B. H. 2001 Large eddy simulation of the gas–liquid flow in a square cross-sectioned bubble column. Chem. Engng Sci. 56, 63416349.CrossRefGoogle Scholar
Elghobashi, S. 2019 Direct numerical simulation of turbulent flows laden with droplets or bubbles. Annu. Rev. Fluid Mech. 51 (1), 217244.CrossRefGoogle Scholar
Erdogan, S. & Wörner, M. 2014 Analysis and modelling of liquid phase turbulence kinetic energy equation in bubbly flows via direct numerical simulation. In 2nd International Symposium on Multiscale Multiphase Process Engineering, Hamburg, September 24–27. MMPE.Google Scholar
Ern, P., Risso, F., Fabre, D. & Magnaudet, J. 2012 Wake-induced oscillatory paths of bodies freely rising or falling in fluids. Annu. Rev. Fluid Mech. 44, 97121.CrossRefGoogle Scholar
Fu, S., Launder, B. E. & Tselepidakis, D. P. 1987 Accommodating the effects of high strain rates in modelling the pressurestrain correlation. In Thermofluids Report TFD/87/5, UMIST, Manchester.Google Scholar
Gibson, M. M. & Launder, B. E. 1978 Ground effects on pressure fluctuations in the atmospheric boundary layer. J. Fluid Mech. 86 (3), 491511.CrossRefGoogle Scholar
Hanjalić, K. & Jakirlić, S. 2002 Second-moment turbulence closure modelling. In Closure Strategies for Turbulent and Transitional Flows (ed. Launder, B. E. & Sandham, N. H.), pp. 47101. Cambridge University Press.Google Scholar
Hanjalić, K. & Launder, B. 2011 Modelling Turbulence in Engineering and the Environment. Cambridge University Press.CrossRefGoogle Scholar
Hanjalić, K. & Launder, B. E. 1976 Contribution towards a Reynolds-stress closure for low-Reynolds-number turbulence. J. Fluid Mech. 74 (4), 593610.CrossRefGoogle Scholar
Horowitz, M. & Williamson, C. H. K. 2010 The effect of Reynolds number on the dynamics and wakes of freely rising and falling spheres. J. Fluid Mech. 651, 251294.CrossRefGoogle Scholar
Ilić, M.2006 Statistical analysis of liquid phase turbulence based on direct numerical simulations of bubbly flows. PhD thesis, Forschungszentrum Karlsruhe, Germany.Google Scholar
Ishii, M. & Hibiki, T. 2006 Thermo-fluid Dynamics of Two-phase Flow. Springer.CrossRefGoogle Scholar
Jakirlić, S. & Hanjalić, K. 2002 A new approach to modelling near-wall turbulence energy and stress dissipation. J. Fluid Mech. 459, 139166.CrossRefGoogle Scholar
Jakirlić, S. & Hanjalić, K. 2013 A direct numerical simulation-based re-examination of coefficients in the pressure strain models in second-moment closures. Fluid Dyn. Res. 45, 055509.Google Scholar
Johansson, A. V. & Hallbäck, M. 1994 Modelling of rapid pressure strain in Reynolds-stress closures. J. Fluid Mech. 269, 143168.CrossRefGoogle Scholar
Jones, W. P. & Launder, B. E. 1972 The prediction of laminarization with a two-equation model of turbulence. Intl J. Heat Mass Transfer 15, 301314.CrossRefGoogle Scholar
Joshi, J. B. & Nandakumar, K. 2015 Computational modeling of multiphase reactors. Annu. Rev. Chem. Biomol. 6 (1), 347378.CrossRefGoogle ScholarPubMed
Kataoka, I., Besnard, D. C. & Serizawa, A. 1992 Basic equation of turbulence and modeling of interfacial transfer terms in gas–liquid two-phase flow. Chem. Engng Commun. 118, 221236.CrossRefGoogle Scholar
Kataoka, I. & Serizawa, A. 1989 Basic equations of turbulence in gas–liquid two-phase flow. Intl J. Multiphase Flow 15, 843855.CrossRefGoogle Scholar
Kenjereš, S., Hanjalić, K. & Bal, D. 2004 A direct-numerical-simulation-based second-moment closure for turbulent magnetohydrodynamic flows. Phys. Fluids 16, 12291241.CrossRefGoogle Scholar
Kim, J., Moin, P. & Moser, R. 1987 Turbulence statistics in fully developed channel flow at low Reynolds number. J. Fluid Mech. 177, 133166.CrossRefGoogle Scholar
Lai, C. C. K. & Socolofsky, S. A. 2019 The turbulent kinetic energy budget in a bubble plume. J. Fluid Mech. 865, 9931041.CrossRefGoogle Scholar
Launder, B. E. 1975 On the effects of a gravitational field on the turbulent transport of heat and momentum. J. Fluid Mech. 67 (3), 569581.CrossRefGoogle Scholar
Launder, B. E. & Li, S. P. 1994 On the elimination of wall-topography parameters from second-moment closure. Phys. Fluids 6 (2), 9991006.CrossRefGoogle Scholar
Launder, B. E., Reece, G. J. & Rodi, W. 1975 Progress in the development of a Reynolds-stress turbulence closure. J. Fluid Mech. 68, 537566.CrossRefGoogle Scholar
Launder, B. E., Tselepidakis, D. P. & Younis, B. A. 1987 A second-moment closure study of rotating channel flow. J. Fluid Mech. 183, 6375.CrossRefGoogle Scholar
Liao, Y., Ma, T., Krepper, E., Lucas, D. & Fröhlich, J. 2019 Application of a novel model for bubble-induced turbulence to bubbly flows in containers and vertical pipes. Chem. Engng Sci. 202, 5569.CrossRefGoogle Scholar
Liao, Y., Ma, T., Liu, L., Ziegenhein, T., Krepper, E. & Lucas, D. 2018 Eulerian modelling of turbulent bubbly flow based on a baseline closure concept. Nucl. Engng Des. 337, 450459.CrossRefGoogle Scholar
Lopez de Bertodano, M., Lahey, R. T. & Jones, O. C. 1994 Development of a k-epsilon model for bubbly two-phase flow. Trans ASME J. Fluids Engng 116, 128.CrossRefGoogle Scholar
Lopez de Bertodano, M., Lee, S.-J., Lahey, R. T. & Drew, D. A. 1990 The prediction of two-phase turbulence and phase distribution phenomena using a Reynolds stress model. Trans ASME J. Fluids Engng 112, 107.CrossRefGoogle Scholar
Lu, J. & Tryggvason, G. 2008 Effect of bubble deformability in turbulent bubbly upflow in a vertical channel. Phys. Fluids 20, 040701.CrossRefGoogle Scholar
Lu, J. & Tryggvason, G. 2013 Dynamics of nearly spherical bubbles in a turbulent channel upflow. J. Fluid Mech. 732, 166189.CrossRefGoogle Scholar
Lumley, J. L. 1978 Computational modeling of turbulent flows. Adv. Appl. Mech. 18, 123176.CrossRefGoogle Scholar
Lumley, J. L. & Newman, G. R. 1977 The return to isotropy of homogeneous turbulence. J. Fluid Mech. 82 (1), 161178.CrossRefGoogle Scholar
Ma, T.2017 A contribution to turbulence modelling in bubbly flows. PhD thesis, Institute of Fluid Mechanics, Technische Universität Dresden, TUDpress, Dresden.Google Scholar
Ma, T., Lucas, D., Ziegenhein, T., Fröhlich, J. & Deen, N. G. 2015a Scale-adaptive simulation of a square cross-sectional bubble column. Chem. Engng Sci. 131, 101108.CrossRefGoogle Scholar
Ma, T., Santarelli, C., Ziegenhein, T., Lucas, D. & Fröhlich, J. 2017 Direct numerical simulation-based Reynolds-averaged closure for bubble-induced turbulence. Phys. Rev. Fluids 2, 034301.CrossRefGoogle Scholar
Ma, T., Ziegenhein, T., Lucas, D., Krepper, E. & Fröhlich, J. 2015b Euler–Euler large eddy simulations for dispersed turbulent bubbly flows. Intl J. Heat Fluid Flow 56, 5159.CrossRefGoogle Scholar
Masood, R. M. A., Rauh, C. & Delgado, A. 2014 CFD simulation of bubble column flows: An explicit algebraic Reynolds stress model approach. Intl J. Multiphase Flow 66, 1125.CrossRefGoogle Scholar
Michaelides, E., Crowe, C. & Schwarzkopf, J. 2016 Multiphase Flow Handbook, 2nd edn. CRC Press.CrossRefGoogle Scholar
Morel, C.1997 Turbulence modeling and first numerical simulations in turbulent two-phase flows. Tech. Rep. CEA/Grenoble, France.Google Scholar
Naot, D., Shavit, A. & Wolfshtein, M. 1970 Interaction between components of the turbulent velocity correlation tensor due to pressure fluctuations. Israel J. Technol. 8, 259269.Google Scholar
Niceno, B., Dhotre, M. T. & Deen, N. G. 2008 One-equation sub-grid scale (SGS) modelling for Euler–Euler large eddy simulation (EELES) of dispersed bubbly flow. Chem. Engng Sci. 63, 39233931.CrossRefGoogle Scholar
Politano, M. S., Carrica, P. M. & Converti, J. 2003 A model for turbulent polydisperse two-phase flow in vertical channels. Intl J. Multiphase Flow 29, 11531182.CrossRefGoogle Scholar
Pope, S. B. 2000 Turbulent Flows, 1st edn. Cambridge University Press.CrossRefGoogle Scholar
Riboux, G., Risso, F. & Legendre, D. 2010 Experimental characterization of the agitation generated by bubbles rising at high Reynolds number. J. Fluid Mech. 643, 509539.CrossRefGoogle Scholar
Rind, E. & Castro, I. P. 2012 Direct numerical simulation of axisymmetric wakes embedded in turbulence. J. Fluid Mech. 710, 482504.CrossRefGoogle Scholar
Risso, F. & Ellingsen, K. 2002 Velocity fluctuations in a homogeneous dilute dispersion of high-Reynolds-number rising bubbles. J. Fluid Mech. 453, 395410.CrossRefGoogle Scholar
Ristorcelli, J. R., Lumley, J. L. & Abid, R. 1995 A rapid-pressure covariance representation consistent with the Taylor–Proudman theorem materially frame indifferent in the two-dimensional limit. J. Fluid Mech. 292, 111152.CrossRefGoogle Scholar
Rodi, W. & Mansour, N. N. 1993 Low Reynolds number k-epsilon modelling with the aid of direct simulation data. J. Fluid Mech. 250, 509529.CrossRefGoogle Scholar
Roghair, I., Mercado, J. M., Van Sint Annaland, M., Kuipers, H., Sun, C. & Lohse, D. 2011 Energy spectra and bubble velocity distributions in pseudo-turbulence: numerical simulations versus experiments. Intl J. Multiphase Flow 37, 10931098.CrossRefGoogle Scholar
Rotta, J. C. 1951 Statistische Theorie nichthomogener Turbulenz. Z. Phys. 129, 547572 (in German).Google Scholar
Santarelli, C.2015 Direct numerical simulations of bubbles in turbulent flows with heat transfer. PhD thesis, Institute of Fluid Mechanics, Technische Universität Dresden, TUDpress, Dresden.Google Scholar
Santarelli, C. & Fröhlich, J. 2015 Direct numerical simulations of spherical bubbles in vertical turbulent channel flow. Intl J. Multiphase Flow 75, 174193.CrossRefGoogle Scholar
Santarelli, C. & Fröhlich, J. 2016 Direct numerical simulations of spherical bubbles in vertical turbulent channel flow. Influence of bubble size and bidispersity. Intl J. Multiphase Flow 81, 2745.CrossRefGoogle Scholar
Santarelli, C., Roussel, J. & Fröhlich, J. 2016 Budget analysis of the turbulent kinetic energy for bubbly flow in a vertical channel. Chem. Engng Sci. 141, 4662.CrossRefGoogle Scholar
Sarkar, S. & Speziale, C. G. 1990 A simple nonlinear model for the return to isotropy in turbulence. Phys. Fluids A 2, 8493.CrossRefGoogle Scholar
Schiller, L. & Naumann, A. 1933 Über die grundlegenden berechungen bei der schwerkraftaufbereitung. Verein. Deutsch. Ing. 77, 318320.Google Scholar
Shi, P. & Rzehak, R. 2018 Bubbly flow in stirred tanks: Euler–Euler/RANS modeling. Chem. Engng Sci. 190, 419435.CrossRefGoogle Scholar
Shir, C. C. 1973 A preliminary numerical study of atmospheric turbulent flows in the idealized turbulent boundary layer. J. Atmos. Sci. 30, 13271339.2.0.CO;2>CrossRefGoogle Scholar
Speziale, C., Sarkar, S. & Gatski, T. 1991 Modelling the pressure–strain correlation of turbulence: an invariant dynamical systems approach. J. Fluid Mech. 227, 245272.CrossRefGoogle Scholar
Sundaresan, S., Ozel, A. & Kolehmainen, J. 2018 Toward constitutive models for momentum, species, and energy transport in gas–particle flows. Annu. Rev. Chem. Biomol. 9 (1), 6181.CrossRefGoogle ScholarPubMed
Tomiyama, A., Tamai, H., Zun, I. & Hosokawa, S. 2002 Transverse migration of single bubbles in simple shear flows. Chem. Engng Sci. 57, 18491858.CrossRefGoogle Scholar
Troshko, A. A. & Hassan, Y. A. 2001 A two-equation turbulence model of turbulent bubbly flows. Intl J. Multiphase Flow 27, 19652000.CrossRefGoogle Scholar
Ullrich, M.2017 Second-moment closure modeling of turbulent bubbly flows within the two-fluid model framework. PhD thesis, Technische Universität Darmstadt.Google Scholar
Ullrich, M., Krumbein, B., Maduta, R. & Jakirlić, S. 2014 An eddy-resolving Reynolds stress model for the turbulent bubbly flow in a square cross-sectioned bubble column. In ASME 2014 International Mechanical Engineering Congress and Exposition, Montreal, Canada, Paper No. 38054. ASME.Google Scholar
Wilcox, D. C. 1998 Turbulence Modeling for CFD, 2nd edn. DCW Industries, Inc.Google Scholar
Ziegenhein, T., Rzehak, R., Ma, T. & Lucas, D. 2017 Towards a unified approach for modelling uniform and non-uniform bubbly flows. Can. J. Chem. Engng 95 (1), 170179.CrossRefGoogle Scholar
Figure 0

Figure 1. Instantaneous DNS data for the case SmMany of Santarelli & Fröhlich (2015). The bubble size is to scale and the vertical plane with the contour plot shows the instantaneous streamwise liquid velocity. (Main part of the picture reprinted from Santarelli (2015) with permission from TUDpress.)

Figure 1

Table 1. Parameters of the cases used for the present study according to Santarelli & Fröhlich (2015, 2016). The labels Sm and La used for the case with the bidisperse bubble swarm, BiDisp, indicate averaging over the small and large bubbles of the case BiDisp, respectively. Here, $N_{p}$ is the number of bubbles, $\unicode[STIX]{x1D6FC}$ the void fraction, $d_{p}$ the particle diameter and $Ar$ the Archimedes number. The values of $Re_{\unicode[STIX]{x1D70F}}$, the friction Reynolds number, $Re_{p}$, the particle Reynolds number based on $d_{p}$ and the relative velocity, and $C_{D}$, the drag coefficient obtained according to (4.4) below, are results of the simulations.

Figure 2

Figure 2. Distribution of $C_{\unicode[STIX]{x1D707}}$ in a $k$$\unicode[STIX]{x1D700}$ type model computed using DNS data.

Figure 3

Figure 3. Shear stress and shear-stress-induced production of TKE for the DNS. (a) Distribution of the structure parameter, $-\overline{\overline{u^{\prime }v^{\prime }}}/k$. (b) Distribution of the ratio of production by mean flow deformation to dissipation of TKE, $\unicode[STIX]{x1D617}_{k}/\unicode[STIX]{x1D700}$.

Figure 4

Figure 4. Comparison between the interfacial term $S_{k}=\frac{1}{2}\unicode[STIX]{x1D61A}_{R,ii}$ according to (3.8) and the a priori evaluation of the pressure–strain terms split into a linear slow term $\unicode[STIX]{x1D719}_{ij,1,linear}$, a nonlinear slow term $\unicode[STIX]{x1D719}_{ij,1,nonlinear}$ and a rapid term $\unicode[STIX]{x1D719}_{ij,2}$ for SmMany, all normalized with $U_{b}^{3}/H$: (a$\unicode[STIX]{x1D719}_{11}$ component, (b$\unicode[STIX]{x1D719}_{22}$ component, (c$\unicode[STIX]{x1D719}_{33}$ component and (d$\unicode[STIX]{x1D719}_{12}$ component.

Figure 5

Table 2. Optimization procedure for the evaluation of $b_{11}^{\ast }$.

Figure 6

Figure 5. Mean flow statistics from the present EE SMC and DNS data for the case SmMany: (a) gas void fraction; (b) liquid streamwise velocity and gas streamwise velocity.

Figure 7

Figure 6. Reynolds-stress components from the present EE SMC and DNS data for the case SmMany: (a) streamwise Reynolds normal stress $\overline{\overline{u^{\prime }u^{\prime }}}$; (b) wall-normal Reynolds normal stress $\overline{\overline{v^{\prime }v^{\prime }}}$; (c) Reynolds shear stress $\overline{\overline{u^{\prime }v^{\prime }}}$. All components are normalized by $U_{b}^{2}$. The label Cokljat refers to EE SMC with $b_{ij}^{C}$ from (4.5), while the label Colombo refers to EE SMC with $b_{ij}^{CF}$ from (4.6) as detailed in the text.

Figure 8

Figure 7. Comparison of key contributions in the present EE SMC with the DNS data for the case SmMany as described in the text. (a) Selected terms in the budget of $\overline{\overline{u^{\prime }u^{\prime }}}$ with all terms normalized by $U_{b}^{3}/H$. (b) Selected terms in the budget of $\overline{\overline{v^{\prime }v^{\prime }}}$ with all terms normalized by $U_{b}^{3}/H$. (c) Selected terms in the budget of $\unicode[STIX]{x1D700}$ with all terms normalized by $U_{b}^{4}/H^{2}$.

Figure 9

Figure 8. Gas void fraction and Reynolds-stress components from the present EE SMC and DNS data for the other two monodisperse cases: (a,b) SmFew; and (c,d) LaMany.

Figure 10

Figure 9. (a) Gas void fraction and (b) Reynolds-stress components from the present EE SMC and DNS data for the case BiDisp.

Figure 11

Table 3. Values of $b_{11}^{\ast }$ and $C_{L}$ obtained by the iterative process for all cases.

Figure 12

Figure 10. Comparison of the present SMC with experimental data of Akbar et al. (2012): (a) gas void fraction, (b) liquid streamwise velocity and gas streamwise velocity, (c) liquid Reynolds normal stresses and (d) liquid Reynolds shear stress.

Figure 13

Figure 11. Anisotropy invariants: (a) the single-phase case Unladen; and (b) the case SmMany.

Figure 14

Figure 12. Anisotropy-invariant map of the DNS data for all cases considered. (a) The cases SmMany, LaMany, BiDisp and Akbar; $Re_{p}$ of the case BiDisp is calculated from the averaged $u_{r}$ and the averaged $d_{p}$ from the DNS of BiDisp. (b) The case SmFew and (c) the single-phase case Unladen.

Figure 15

Figure 13. Values of $b_{11}^{\ast }$ as a function of $Re_{p}$ for all cases considered. The dashed line shows the fit.

Figure 16

Figure 14. Anisotropy-invariant map of the present SMC for all bubble-laden cases considered.

Figure 17

Table 4. DNS parameters of the deformable bubble case in Lu & Tryggvason (2008).

Figure 18

Figure 15. Comparison of the present SMC with DNS data for the deformable bubble case in Lu & Tryggvason (2008): (a) gas void fraction; (b) liquid streamwise velocity and gas streamwise velocity; (c) liquid streamwise and wall-normal velocity fluctuations normalized by $u_{\unicode[STIX]{x1D70F}}$; and (d) liquid Reynolds shear stress normalized by $u_{\unicode[STIX]{x1D70F}}^{2}$.

Figure 19

Figure 16. Comparison of the SSG with DNS data for the single-phase case in Lu & Tryggvason (2008): (a) streamwise and wall-normal velocity fluctuations normalized by $u_{\unicode[STIX]{x1D70F}}$; and (b) Reynolds shear stress normalized by $u_{\unicode[STIX]{x1D70F}}^{2}$.

Figure 20

Figure 17. Geometrical definition of proper Euler angles (standard $y$-convention). The $xyz$ (fixed) system is shown in blue; the $XYZ$ (rotated) system aligned with the local relative velocity in $X$-axis is shown in red. The line of nodes ($N$) as the intersection of the planes $yz$ and $YZ$ is shown in green.