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

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.

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 2006) 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 (2015), Michaelides, Crowe & Schwarzkopf (2016) and Sundaresan, Ozel & Kolehmainen (2018). 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 2010; Elghobashi 2019). 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 1994; Morel 1997; Troshko & Hassan 2001; Colombo & Fairweather 2015; Ma et al. 2017). 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 1989; Troshko & Hassan 2001). These studies, as well as own experience (Ma 2017), 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 2011), 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 2010; Akbar et al. 2012; Lai & Socolofsky 2019) and direct numerical simulation (DNS) (Lu & Tryggvason 2008; Bolotnov et al. 2011; Santarelli & Fröhlich 2015), 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. (1990) were the first to address this issue in bubbly flow. They adopted the SMC of Launder, Reece & Rodi (1975) – hereafter referred to as LRR – as the carrier model and combined it with the algebraic BIT expression of Biesheuvel & van Wijngaarden (1984) 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 2002).

Other studies (Masood, Rauh & Delgado 2014; Ullrich et al. 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 (2001) 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. (2001) have also been reproduced very well using scale-resolving simulations without accounting for BIT (Deen et al. 2001; Niceno, Dhotre & Deen 2008; Ma et al. 2015a). 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. 2015a). Further EE second-moment modelling for bubbly flows was performed by Cokljat et al. (2006) and Colombo & Fairweather (2015), and recently by Ullrich (2017). While Colombo & Fairweather (2015) 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 (1992). 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. 2011; Roghair et al. 2011; Dabiri, Lu & Tryggvason 2013; Lu & Tryggvason 2013; Santarelli & Fröhlich 2016; Cifani, Kuerten & Geurts 2018), 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ć (2006) and Erdogan & Wörner (2014) 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 (2016) computed the budget terms in the TKE equation and proposed a BIT closure based on a priori evaluations. Ma et al. (2017) 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 (2015, 2016) 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. 2017). 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 (2008) 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. 2016). 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 (2015, 2016) 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ć 2006; Erdogan & Wörner 2014), 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 (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.)

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.

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

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. (2012), 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 1972)

(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 1989; Kataoka et al. 1992). 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. (2016), 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 2000). 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 (1993) for the channel data of Kim, Moin & Moser (1987) 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 (1993) 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 2006) 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 (1991) is used for the wall force, the model of Burns et al. (2004) for turbulent dispersion, and the model of Auton (1987) 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. (1992) 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 (2011) 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 (1973), 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 2000; Hanjalić & Jakirlić 2002; Hanjalić & Launder 2011). 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 2000).

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. (2017) 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. (2017), 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. 2017).

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 1975), rotation (Launder, Tselepidakis & Younis 1987) and electromagnetic forces (Kenjereš, Hanjalić & Bal 2004) 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 2011).

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 2011). Hence, modelling this term can safely start from the most general formulation in single-phase flow according to the Cayley–Hamilton theorem (Lumley 1978; Fu, Launder & Tselepidakis 1987; Speziale, Sarkar & Gatski 1991)

(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. 1975; Gibson & Launder 1978; Speziale et al. 1991; Ristorcelli, Lumley & Abid 1995). More related models and the corresponding detailed forms of the model coefficients can be found in the summary of Jakirlić & 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 (1951) 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 1978; Fu et al. 1987; Speziale et al. 1991; Hanjalić & Jakirlić 2002). Jakirlić & 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 (1990) and Speziale et al. (1991) – 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 (1970),

(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 4a–c), the nonlinear part of the slow term $\unicode[STIX]{x1D719}_{ij,1,nonlinear}$ is smaller but of the same order as the linear part $\unicode[STIX]{x1D719}_{ij,1,linear}$ and has overall the opposite sign to the latter. The influence of the rapid term $\unicode[STIX]{x1D719}_{ij,2}$ appears to be negligible for this BIT-dominated flow in the channel centre for all components. This is expected, as it is related to the production term $\unicode[STIX]{x1D617}_{ij}$, which is small in the centre due to the vanishing mean flow gradients. All contributions to the pressure–strain term have a much smaller magnitude in the Reynolds shear stress component. Among these, the linear part of the slow term is dominant in the channel centre (see figure 4d). The behaviour of the total term $\unicode[STIX]{x1D719}_{ij}$ will be discussed in a later section.

The analysis of the pressure–strain terms for the other three cases (SmFew, LaMany and BiDisp) reveals the same trends as in figure 4, so that they are not reproduced here. The observations based on the DNS data hence support two guidelines for modelling the pressure–strain term in bubbly flow. First, due to the importance of the nonlinear part $\unicode[STIX]{x1D719}_{ij,1,nonlinear}$, the slow term $\unicode[STIX]{x1D719}_{ij,1}$ requires a nonlinear model. Second, it indicates that the rapid term is very small, so that the specific form of the model is uncritical and the simple model (3.19) sufficient. Furthermore, more elaborate models extending (3.19) are heavily based on the mean flow gradient (Johansson & Hallbäck 1994), 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ć 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. 2012; Santarelli & Fröhlich 2015). 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 2003; Ziegenhein et al. 2017; Liao et al. 2018, 2019). 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 1989) 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. (2016) 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. 2016), 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 (2015). 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. (2006) 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 (2015) 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. (2006) and Colombo & Fairweather (2015) 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 2011). 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 1976; Wilcox 1998) and multiphase flows (Ma 2017), 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. (2012): (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. (2012). 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. 2012). 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. (2019) 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. (2017). 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. 2015b). 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 (1977), 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 1978),

(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 (1977) 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 6, 8b, 8d and 9b) compared to the other two directions over a large range of the channel centre for the bubble-laden cases. Such a major contribution of normal stress in the streamwise direction was also found by a recent DNS study by du Cluzeau, Bois & Toutant (2019), 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 (2016)) and the path of the rising bubbles (Horowitz & Williamson 2010; Ern et al. 2012).

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 2004; Rind & Castro 2012). 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 2012) 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 (2004) and Rind & Castro (2012) 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 1994; Ristorcelli et al. 1995; Jakirlić & 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 (2008)

In this section, the results of the EE SMC using (5.3) are compared to the DNS bubbly channel data of Lu & Tryggvason (2008). 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 2008). 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 (2008).

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

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 (2008)). 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 1998; Hanjalić & Launder 2011). 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 (2008). 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 (2015) 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.