Assessment of turbulence models for bubbly flows: toward a Five-Equation turbulence model

Journal Pre-proofs Assessment of turbulence models for bubbly flows: toward a Five-Equation turbulence model Ghazi Bellakhal Ghazi Bellakhal, Fathia Chaibina, Jamel Chahed PII: DOI: Reference:

S0009-2509(19)30915-7 https://doi.org/10.1016/j.ces.2019.115425 CES 115425

To appear in:

Chemical Engineering Science

Received Date: Revised Date: Accepted Date:

12 July 2019 7 December 2019 9 December 2019

Please cite this article as: G. Bellakhal Ghazi Bellakhal, F. Chaibina, J. Chahed, Assessment of turbulence models for bubbly flows: toward a Five-Equation turbulence model, Chemical Engineering Science (2019), doi: https:// doi.org/10.1016/j.ces.2019.115425

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Ltd.

Assessment of turbulence models for bubbly flows: toward a FiveEquation turbulence model Ghazi Bellakhal*, Fathia Chaibina and Jamel Chahed University of Tunis El Manar, National Engineering School of Tunis, BP N37, Le Belvedere 1002, Tunis, Tunisia Corresponding

author Email addresses:[email protected] (Ghazi Bellakhal), [email protected] (Fathia Chaibina), [email protected] (Jamel Chahed)

Abstract The paper presents an assessment of turbulence models in dispersed bubbly flows based on an analysis of their application to uniform and uniformly sheared homogeneous turbulence. This analysis leads to the development of a five-equation model where the total turbulent fluctuations in bubbly flow is submitted to a double decomposition: First, the turbulent fluctuations are split into Shear-Induced-Turbulence (SIT) and BubbleInduced-Turbulence (BIT) components. Secondly, The BIT contribution is split into two components corresponding to non-dissipative pseudo-turbulent part and dissipative turbulent part produced by local shear in bubble’s wake. Consequently, the total turbulent energy is split into three contributions : two dissipative quantities k t and k s produced respectively by gradient of mean velocity in the liquid phase and by local shear in bubble’s wake both are described by their modeled transport equations as well as their dissipation rates  t and s ; and a non-dissipative component k p attributed to fluctuations of potential flow around bubbles described by its modeled transport equation. This model was implemented in CFD code and validated in homogeneous bubbly turbulence situations. The numerical results show a satisfactory agreement with experimental data of turbulence.

-1-

Keywords Bubble-Induced-Turbulence ; Shear-Induced-Turbulence ; Bubbly flows ; Homogeneous Turbulence ; Five-Equation turbulence model ; Multiphase CFD

1. Introduction Multiphase flows are encountered in a wide variety of industrial processes as well as natural phenomena. Among them, bubbly flows are largely involved in industrial applications because they combine the advantages of providing large interfacial area by unit of volume and of inducing important liquid agitation. These factors enhance the interfacial mass transfer as well as the chemical reagents mixing and accelerate the chemical reactions. Direct Numerical Simulation (DNS) is yet difficult to apply to industrial scale twophase systems and Computational Fluid Dynamics (CFD) has been and continues to be an important tool to predetermine complex turbulent flows. When applied to two-phase flows, CFD is based on an Eulerian-Eulerian approach, that uses Ishii formulation (Ishii, 1975; Drew, 1983; Ishii and Hibiki, 2011). According to this approach, each phase is considered as a continuous phase through the whole flow domain for which the dynamic is calculated using Reynolds-Averaged Navier–Stokes equations (RANS) of both twophases. The resulting “two-fluid models” require closure of the interfacial terms in the momentum equations traducing the interfacial momentum transfer between the two phases as well as the turbulence in both the two phases (Lopez De Bertodano et al., 1990 - 1994 ; Chahed et al., 2002 - 2003 ; Rzehak and Kriebitzsch, 2015 ; Rsehak et al., 2017 ; Guan et al., 2015 ; ... ). Many experimental works (Serizawa et al., 1975a b c ; Wang et al., 1987 ; Liu and Bankoff, 1992a) showed the important effects of interfacial transfer due to the bubble’s presence on the turbulence of the liquid phase. On the other hand, many authors evoked the role of the turbulence on the bubbles size, shape, distribution and dynamics (Zerizawa et al., 1975b ; Liu and Bankoff, 1992b ; Kawamura and Kodama, 2002 ; Kato et al., 1999; Colin et al., 2012 ; Chahed et al., 2016 ; Rezig et al., 2017 ; Aouadi et al., 2019). On the

-2-

whole, these results indicate that interfacial transfer and turbulence are fundamentally coupled in bubbly flows. Turbulence modeling is thus an important issue in developing accurate Eulerian-Eulerian two-fluid models for bubbly flows. Dispersed bubbly flows could roughly be classified under two categories: bubbly flows where the liquid phase flow is driven mechanically and bubbly flows where the fluid motion is induced by the bubbles (bubble-driven flows). The first category is often characterized by relatively high average liquid velocity, generated by external mechanical support (e.g. pipe bubbly flows, bubbly flows in agitated reactors). Similarly to the singlephase turbulent flow, turbulence production occurs in the sheared zones of the flow (Shear-Induced-Turbulence, SIT). On the other hand, the bubbles, in their relative motion, generate, at their own scale, specific turbulent fluctuations (Bubble-Induced-Turbulence, BIT), which may alter the mean velocity and the turbulence structure of the continuous phase (Balachandar and Eaton 2010). In the case where Shear-Induced-Turbulence (SIT) (the one that could exist without bubbles) is dominant (this occurs in highly sheared bubbly flows with low void fraction and small bubbles), the turbulence is mainly controlled by mechanisms similar to those that occur in single-phase flows. On the other hand, the main source of momentum in bubble-driven flows is the buoyancy effect due to bubble’s presence (e.g. bubble-column, air-lift). Without mechanical support, liquid velocities in these flows are often relatively low and fluctuations are mainly induced by bubbles relative motion (BIT). Nonetheless, nonuniform void fraction distribution induces a contrast of mixture density which may generate recirculating flow with mean velocity gradient in the continuous phase (Hassan et al., 1992 ; Mudde, 2005 ; Simiano et al., 2009 ; Simiano and Lakehal, 2012) leading to shear production of turbulence (SIT). Here again the two kind of the turbulent fluctuations (SIT and BIT) coexist in turbulent bubbly flows and may interact (Colombet et al., 2015). We deduce therefore that turbulence in dispersed bubbly flows is governed by different but so strongly coupled dynamics and involves different scales. Indeed, many experimental results studies showed that the presence of the bubbles alters considerably the liquid turbulence structure in homogeneous uniform and uniformly sheared bubbly flows (Lance and Bataille (1991) ; Lance et al. (1991)). Bubbles agitation enhance considerably the turbulence of the liquid in homogeneous uniform bubbly flows and

-3-

provoke an attenuation of the shear stress in homogeneous uniformly sheared turbulence. Similar results are reported from wall-bounded bubbly flows experiments where the turbulence of the liquid is amplified in the low sheared regions (away from the wall), when, in the highly sheared zones, a reduction of the shear stress is observed with, in some cases, an attenuation of the turbulent intensity in comparison with the equivalent single-phase flows (Moursali et al., 1995 ; Wang et al., 1987). 1.1. Case studies of uniform and uniformly sheared homogeneous bubbly flows Homogeneous turbulent bubbly flows (uniform or uniformly sheared) make it possible to analyse the effect of bubbles on turbulence while avoiding the difficulties related to the non-uniform distribution of void fraction as well as to the coupling between mean velocity and turbulent fluctuations in liquid phase. The experiments carried out by Lance and Bataille (1991) and Lance et al. (1991) are, in this respect, of central importance. Concerning the homogeneous turbulence in an upward uniform bubbly flow behind a grid, the authors have investigated, single-phase and three bubbly cases (with characteristic void fractions : 0.5% ; 1% and 2%) where the mean liquid velocity is maintained constant equal to 0.6 m/s. In these experiments the Shear-Induced-Turbulence (SIT) is generated by means of a grid placed at the bottom of a vertical square test section. Bubbles are injected uniformly through the grid by means of equally spaced injectors (Lance and Bataille, 1991). The injection system generates a mono-disperse millimetric population of bubbles characterised by mean diameter equal to 5 mm. As regards homogeneous turbulence with uniform shear, the experiments produce a constant uniform shear, equal to 2.9 s-1, both in single-phase and in three bubbly cases (with void fractions of 1% ; 1.4% and 2%) (Lance et al., 1991). These experiments allowed the authors to analyse the effect of the bubble’s presence on the turbulent shear production as well as on the anisotropy of the liquid turbulence. The results obtained from the uniform homogeneous bubbly turbulence experiments indicate a similar decay of the turbulence intensity behind the grid in bubbly cases as that observed in the single-phase case but starting from a more important level as it is shown in Figure-1. These experimental results suggest that the decays of turbulence behind the grid in bubbly flow cases occur with the same characteristic time scale as in single-phase case.

-4-

Figure-1 : experimental data of turbulence intensity behind a grid in uniform homogeneous turbulent flow in single-phase case and three bubbly cases (Lance and Bataille, 1991) Experimental data obtained in uniformly sheared homogeneous bubbly turbulence show two main results:

a significant decrease of turbulent shear in bubbly cases when

compared to the corresponding single-phase case and a more pronounced trend toward isotropy in two-phase cases as shown in Figure-2 and Figure-3.

-5-

l Figure-2: Experimental data of Reynolds stress tensor longitudinal evolution in uniform sheared homogeneous turbulent flow in singlephase case (Lance et al., 1991)

Figure-3: Experimental data of Reynolds stress tensor longitudinal evolution in uniform sheared homogeneous turbulent flow in twophase case (α = 1.4%) (Lance et al., 1991)

-6-

Before to be applied in simulating complex two-phase dispersed turbulent bubbly flows, turbulence modeling in bubbly flows should be able to reproduce adequately the behaviors observed in the basic situation of homogeneous bubbly turbulence. In particular, turbulence closures should be able to reproduce a decay of turbulence behind grid in homogeneous uniform bubbly flows as well as the decrease of the turbulent shear observed in homogeneous uniformly sheared bubbly in comparison to the corresponding single-phase case. 1.2. Turbulence modeling approaches in bubbly flows Modeling of turbulence in bubbly flows appears to be so far from an extent of the single-phase turbulence analysis and modeling. From literature, we may distinguish three main orientations in developing turbulence models for bubbly flows giving rise to three classes of turbulence models. The first and second classes of turbulence models do not distinguish between SIT and BIT contributions and gives rise to turbulence models of the same type as those developed for single-phase flows. In the resulting two-equation (k, ε) models, the transport equations of the kinetic turbulent energy (k) and of its dissipation rate (ε) are modelled as in single phase flow. To take account for the effect of the bubbles, additional source terms are introduced in transport equations of k and ε. However, whereas the turbulent viscosity in the first class of turbulence models is calculated as in single-phase flow (Kataoka et al., 1989 ; Lee et al., 1989 ; Simonin et Viollet, 1990a-b ; Simonin, 1990 ; Rzehak and Krepper, 2013 ; Colombo and Fairweather, 2015; …), a supplementary turbulent viscosity associated to bubbles agitation is introduced in the second class of turbulence models introducing, along with, a time scale related to BIT (Sato et al, 1981). The third class of turbulence models is based on the decomposition of the turbulent fluctuations into two contributions : the turbulent components of the Reynolds stress tensor are calculated as the sum of a Shear-Induced-Turbulence (SIT) and Bubble-Induced-Turbulence (BIT) (Lance et al., 1991 ; Lopez de Bertodano et al., 1994 ; Chahed et al., 2003 ; Guan et al., 2015 ; Bellakhal et al., 2004b ; Risso, 2018). In a recent review, Risso (2018) highlighted that turbulence closure separating the Shear-Induced-Turbulence (SIT) and Bubble-Induced-Turbulence (BIT) are more likely in a position to develop models able to describe the physical phenomena and the specific scales relating to each of the contributions as well to account for their possible couplings.

-7-

However, the development and implementation of such turbulence models remains limited and the main difficulty lies in properly establishing criteria of separation. Chahed et al. (2003) proposed to split the Reynolds stress tensor in turbulent dissipative part and pseudo-turbulent non-dissipative part. The turbulent dissipative part, which has the fluctuating character of ordinary turbulence subject to spectral cascade transfer leading to a dissipation process, includes the SIT and part of the BIT (the part produced by the local shear and dissipated in the bubbles wakes). The pseudo-turbulent non dissipative part of the BIT refers to the velocity fluctuations in the liquid resulting from the displacement of liquid by the moving bubbles which does not participate to the turbulent cascade and do not contribute to dissipation. Similar decompositions of BIT are proposed to distinguish between turbulent and non-turbulent contributions associated to the fluctuations induced by the moving bubbles (Risso, 2018 ; Du Cluzeau et al., 2019). There are broadly two methods used in modeling for splitting the Reynolds stress tensor into two components for which specific closures are developed. There is on the one hand the decomposition of the Reynolds stress tensor into pseudo-turbulent nondissipative and turbulent dissipative parts. According to this separating method adopted by Lance et al. (1991), Chahed et al. (2003), Bellakhal et al. (2004a,b), Chaibina et al. (2019), the turbulent part contains two components related to two different time-scales: the Shear-Induced-Turbulence (SIT) (the one that could exist without bubbles) and the turbulence produced by the local shear in the bubble’s wake. The two components correspond to two different scales. On the other hand, there is the decomposition method based on the distinction between SIT and BIT contributions adopted by Lopez de Bertodano et al. (1994) and Guan et al. (2015). BIT includes, in this case turbulent dissipative and pseudo-turbulent non-dissipative contributions. Although related to the same timescale (that of the bubble relative motion), turbulent and pseudo-turbulent fluctuations are controlled by two fundamentally different phenomena. Within these hypotheses, several researches lead to the development of three transport equations turbulence models (Lopez de Bertodano et al., 1994; Bellakhal et al, 2004b ; Guan et al., 2015). Moreover, it is so hard to produce experimental data allowing to separate the two kinds of turbulent fluctuations in a turbulent bubbly flow. However, recent and ongoing experimental and numerical researches on Bubble-Induced-

-8-

Turbulence (Amoura et al., 2010 ; Mrabtini et al., 2017 ; Risso, 2018 ; Du Cluzeau et al., 2019) provide new insights not only for turbulence models assessments but also for supporting new ideas in developing more accurate turbulence models for bubbly flows. 1.3. Objectives of this work In line with this perspective, the present paper proposes a new turbulence model based on a double decomposition of the turbulence in bubbly flows. Turbulence in bubbly flows is first split into SIT and BIT contributions. The latter is then decomposed into two parts: a turbulent dissipative part and pseudo-turbulent non-dissipative part. The turbulent kinetic energy is thus divided into three separated contributions: two turbulent-dissipative components ( k T and k S ) produced respectively by the mean velocity gradient of the liquid flow and by the local shear in bubbles’ wakes ; and a pseudo-turbulent part k P corresponding to the non-dissipative contribution of BIT. Following a similar approach to that adopted in developing three-equation turbulence model (Chahed et al., 2003 ; Bellakhal et al., 2004 a-b ; Chaibina et al., 2019), we propose to develop separate closures for the transport equations of each component of the turbulent kinetic energy. This gives rise to a five-equation turbulence model: two k-epsilon (two-equation) turbulence models for the two turbulent dissipative parts and a one-equation turbulence model for the pseudo-turbulent non-dissipative part. As mentioned above, turbulence and phase distribution are coupled. To avoid difficulties induced by non-uniform bubbles distribution, we consider here homogeneous turbulent bubbly flows with specific reference to the results obtained in uniform bubbly flow (Lance and Bataille, 1991) and in uniformly sheared bubbly blows (Lance et al., 1991). Recall that an increase or a reduction of turbulent statistical qualities were observed in uniform and uniformly sheared homogeneous bubbly flows turbulence in comparison with the equivalent single-phase liquid flows. In the absence of heterogeneity related to the void fraction the experiments of homogeneous turbulence constitute an important step in the analysis and modeling of bubbly turbulence. The diffusion may be neglected in this case, and the analysis can thus be focused on the mechanisms that control turbulence budget of each components of the turbulent fluctuations and their interactions.

-9-

2. Assessment of turbulence models for bubbly flows: case study of homogeneous turbulence 2.1. Turbulence modeling in bubbly turbulent flows: a classification proposal As mentioned above, three main orientations may be distinguished in developing turbulence models for bubbly flows giving rise to three classes of models. For the two first classes, turbulence kinetic energy and its dissipation rate are maintained as a unique entities (Kataoka et al., 1989 ; Lee et al., 1989 ; Simonin et Viollet, 1990a-b ; Simonin, 1990 ; Morel, 1995 ; Ben Magolan et al, 2019). The common forms of their transport equations are identical to the single-phase case where the coefficients keep the same values. They write as follows:

 1    k    . 1    v k     t      . 1      L  t k  

    k   1     P     Sk  

(1)

 1      . 1    v     t       . 1     L  t     1       C1 P  C2   S k   k  

(2)

Where the turbulence production rate term by gradient of mean liquid velocity P is expressed using turbulent viscosity concept as:

P 

 

 

t   t  v .   v   v   

(3)

The BIT is taken into account implicitly by means of two additive source terms Sk and S . The Sk term describes the turbulence production by interfacial interactions. It is interpreted as the transfer rate of bubble energy to the turbulent kinetic energy in its wake by means of the drag force and is modeled as:

Sk







vG  v  FDrag



3 3 CD  vG  v 4 dB

-10-

(4)

Where FDrag is the drag force and C D is its coefficient (Tomiyama et al., 1998). It should be noted that Simonin (1990) has suggested a distinguished model for the production source term Sk involving the drift velocity. It may be stated as follows :

Sk 

3 CD  vR 4 dB

v

d

vR  vR  v



(5)

Where vd is the drift velocity interpreted as a correlation between liquid velocity fluctuations and the instantaneous gas phase presence G . It is given by the following relation :

 vd 

 G v 

with < . > an appropriate averaging operator

(6)

Hence, the relative velocity for this modeling approach is expressed as follows :

vR



v

G



 vL  vd

(7)

as regards the dissipative source term S , it may be interpreted as a supplementary dissipation induced by bubbles agitations. It is modeled for the whole of modeling approaches following a similar approach as in single-phase case by relating it to Sk via the characteristic time scale  

S



C3

k : 

Sk 

(8)

The coefficient C 3  allows to calibrate interfacial production and dissipation of turbulence for this class of models. It was proposed by Lee et al. (1989) that this coefficient would be adjusted to the coefficient C 2  as :

C3



C2



1.92

(9)

While for the first class of models the turbulent viscosity is still modeled similarly to single-phase case as follows (Lee et al., 1989 ; Morel, 1995 - 1997 ; Rzehak and Krepper,

-11-

2013 ; Politano et al, 2003 ; Colombo and Fairweather, 2015 ; Ma et al., 2017 ; Ziegenhein et al., 2017) : t



c

k2 

(10)

It will be modeled for the second class following the formulation suggested by Sato et al. (1981) (even though this author kept for the kinetic turbulent energy and its dissipation rate modeling transport equations identical to those used in single-phase case). The BIT is accounted for at turbulent viscosity modeling by adding to the effective turbulent viscosity (modelled as in single-phase case) a supplementary term attributed to the dispersed phase characteristics. The adopted formulation of the viscosity turbulent is expressed as follows (Sato et al., 1981 ; Krepper et al., 2005 ; Lahey, 2005 ; Tabib et al, 2008 ; Rzehak and Krepper, 2013 ; Vaidheeswaran and Hibiki, 2017 ; Magolan et al., 2019) :

t 

SIt  b

(11)

Where : SIt



c

k2 

and

b 

Cb dB vG  v

(12)

The coefficient Cb was adjusted by Sato et al. (1981) at the value 0.6. Regarding of several reviewing works (Vaidheeswaran and Hibiki, 2017 ; Ziegenhein et al. 2017 ; Magolan et al., 2019) , the first class of models is widely used in two-phase CFD practical simulations. Nonetheless, the second class is also adopted in many CFD simulation and investigation developments of bubbly flows in both opposite situations of highly sheared bubbly flows and bubble-driven flows. Although it gives good agreement with experimental data in many cases as pipe flows (Krepper et al., 2005) or column bubbles (Tabib et al., 2008), it underestimates the kinetic turbulent energy in similar cases where experimental conditions were modified (Krepper et al., 2011 ; Ziegenhein et al., 2015). The third class of turbulence models is founded upon the decomposition of total turbulent fluctuations into two parts: SIT and BIT. Nonetheless, regarding the two kinds of fluctuations contributing in BIT as evoked at first by Lance and Bataille (1991) and

-12-

then followed by many authors (Chahed et al., 2003 ; Risso, 2018 ; du Cluzeau at al., 2019). The latter entity may be seen as the sum of pseudo-turbulent non-dissipative fluctuations and local shear produced dissipative fluctuations in bubble’s wake. Hence, the Reynolds stress tensor characteristic of BIT may be divided into two statistically independent parts corresponding to the two contributions as follows:

TRe_BIT 

P S TRe_BIT + TRe_BIT

(13)

Bisheuvel and Van Vijingaarden (1984) have developed an analytical expression for the pseudo-turbulent contribution in homogeneous situation by averaging a potential flow around a group of spheres. This anisotropic expression is written as:

T



P Re_ BIT H



3 1   uR  uR  uR 20  20

2

 I 

(14)

According Chahed et al. (2003) modeling approach, a phenomenological separating mode is adopted for which all turbulent fluctuations produced by shear are included in a unique entity even they may be described by different characteristic time scales. Therefore, the two dissipative quantities of fluctuations produced by local shear in bubble’s wake and SIT contribution produced in continuous phase by gradient of mean liquid velocity are included in the same component as explained schematically in Figure-4.

-13-

Figure 4 : Schematic decomposition of total turbulent fluctuations proposed by Chahed et al. (2003) Lance and Bataille (1991) have analyzed the alteration of the spectral transfer of total turbulent energy homogeneous turbulent bubbly flow. They have deduced that there is an equilibrium between local turbulence production and dissipation in the bubbles’ wakes. Based on this experimental result and according to the turbulence model proposed by Chahed et al. (2003), the dissipative part of BIT does not appear in the transport equation of SIT component, which contains the turbulence, produced in the bubbles’ wake. Consequently, at first order, the total turbulent energy is split into two contributions : dissipative turbulent energy k T and non-dissipative pseudo-turbulent energy k P as follows:

k 

k P  kT

(15)

Such modeling approach had given rise to a three-equation turbulence model (Bellakhal et al., 2004b ; Chaibina et al., 2019) : two transport equations of the dissipative turbulent part of kinetic energy k T (equation (16)) and its dissipation rate  T (equation (17)) modeled similarly as in single phase case (without any additional interfacial term source):

-14-

 1    k T    . 1    v k T     t   . 1    t k T  b k P   k T   1     PT  T   1    T    . 1    v T     t   . 1    t k T  b k P   T  

1    t

C1 PT

 C2 T 

(16)

(17)

Where :

PT 

 

 

t   t  v .   v   v   

(18)

The transport equation of non-dissipative pseudo-turbulent part k P (equation (20)) is modeled as a generalized form of the potential solution in homogeneous case proposed by Bisheuvel and Van Vijingaarden (1984). In fact, according to the algebraic expression of the pseudo-turbulent Reynolds stress tensor established by those authors in homogeneous situation given by equation (13), the pseudo-turbulent energy kPH in homogeneous turbulent bubbly flow situation may be deduced therefore analytically as follows:

k PH





1 P trace TRe_ BIF 2

 H



1  uR 4

2

(19)

In inhomogeneous situation, the pseudo-turbulent energy kP transport equation modeling was achieved by adding a supplementary diffusive term attributed to the nonhomogeneous character of flow. Such equation is thus written as follows:

 1    k P    . 1    v k P     t 





 k   . 1    t k T  b k P   k P    PH   . v k PH   t 

-15-

(20)

The turbulent viscosity formulation is deduced by a reduction of second order turbulence closure (Bellakhal et al., 2004a). It involves, two characteristic time scales of SIT component t and bubbles agitation b :

t 

kT T

and

b  CR

dB

(21)

uR

The turbulent viscosity is expressed as follows:

t



c

k T2 T

cb k P c 0 k T  1  t b

1

(22)

The diffusion terms given by first ones in right hand side of equations (16), (17) and (20) are modeled using a gradient law (following Launder et al. (1975) approach in singlephase case) with diffusion coefficient that includes two effects : turbulent diffusion expressed by t kT and a supplementary diffusion induced by bubbles agitation expressed by bkP . The Reynolds stress tensor is calculated using the following Boussinesq closure (Chaibina et al., 2019) :

 

t 2   TRe  TRe_ BIT + TRe_SIT  t   v   v    kP  kT  I 3  

(23)

Another modeling approach founded upon a different separating mode was proposed by Lopez de Bertodano et al. (1994) and Guan et al. (2015). According to this approach, the local wake shear produced dissipative fluctuations are included in BIT contribution as explained schematically in Figure-5:

-16-

Figure 5 : Schematic decomposition of total turbulent fluctuations proposed by Lopez de Bertodano et al. (1994) Such separating mode includes, in the same variable, two kinds of turbulent fluctuations: non-dissipative pseudo-turbulent fluctuations and dissipative fluctuations produced in the bubbles’ wake. Even though these two quantities may be associated to the same characteristic time scale, they are characterized by two different production mechanisms as well as different anisotropic structures. According to this modeling approach, the total turbulent energy is split into two dissipative contributions as follows:

k 

kBI  kSI

(24)

The kSI part and its dissipation rate  SI are modeled using standard single-phase transport equations:

 1    kSI   . 1    v kSI     t       . 1     L  t   kSI   1     PSI  SI  k    

-17-

(25)

 1    SI    . 1    v SI     t          . 1      L  t   SI   1     SI   C1 PSI  C2 SI  k    k SI   

(26)

Where the production by gradient of mean velocity term PSI is calculated using turbulent viscosity as follows:

PSI 

 

 

t   t  v .   v   v   

(27)

Concerning the k BI part, its transport equation is modeled by Lopez de Bertodano et al. I

(1994) by means of an additive source term Sk BI traducing turbulence interfacial production and given by equation (29). The Bubble Induced Turbulent energy is supposed a

to reach the asymptotic value in homogeneous situation k BI using a relaxation term associated to the bubble characteristic time scale B . The asymptotic value (noticed k BI ) a

corresponds to the pseudo-turbulent energy in homogeneous situation kPH proposed by Bisheuvel and Van Vijingaarden (1984) as described above, (equation 19).     1    k BI    . 1    v k BI    . 1      L  t   t k  

  I   k BI   Sk BI  

(28)

Where :

SkI BI

1 a  kBI  k BI B





;

B

 2  CVM   dB    3  CD   v  v  G

    

(29)

The coefficient CVM appearing in B corresponds to the virtual added mass force coefficient which takes the standard value equal to 1 2 . According to Guan et al. (2015), the transport equation of k BI is modeled by means of two additive source terms traducing interfacial production and dissipation of turbulence. This transport equation is written as follows:

-18-

 1    k BI   . 1    v k BI     t       . 1     L  t   k BI   PBI  1    BI k    

(30)

Where the interfacial production term PBI is interpreted as the transfer rate of bubble energy to the turbulent kinetic energy in its wake by means of the drag force and is modeled similarly as mentioned for the first class of models (equation (4)) as follows:

PBI



3 CD   vG  v 4 dB

3

(31)

The interfacial dissipation rate  BI , is modeled using an algebraic formulation given by :

BI 

CD

 kBI 

3 2

(32)

dB

In both models of Lopez de Bertodano et al. (1994) and Guan et al. (2015), the turbulent viscosity is modeled as the sum of two contributions :

t 

SIt  BI t

(33)

Where t is attributed to SIT contribution and expressed as: SI

SIt



c

kSI2 SI

(34)

The Bubble Induced Turbulent viscosity t is modeled by Lopez de Bertodano et al., BI

(1994) using Sato et al. (1981) model (equation (12)), while Guan et al. (2015) has proposed an expression involving a velocity scale deduced from Shear-InducedTurbulent energy :

 BI t



cW

k SI d B

(35)

The Reynolds stress tensor is calculated by Guan et al. (2015) using an isotropic Boussinesq closure similar to the precedent relation (23) written as:

-19-

 

t 2   TRe  TRe_ BIT + TRe_SIT  t   v   v    kBI  kSI  I 3  

(36)

However, Lopez de Bertodano et al. (1994) has taken account for the anisotropic character of the pseudo-turbulent fluctuations described in homogeneous potential situation by solution (14) of Bisheuvel and Van Vijingaarden (1984). To describe this anisotropy, we introduce a spherical tensor A. Considering the z-axis as the upward vertical direction, the spherical tensor may be written as:

0  3 5 0   A   0 35 0   0 0 4 5  

(37)

Lopez de Bertodano et al (1994) express the Reynolds stress tensor by:

 

t 2   TRe  TRe_ BIF + TRe_SIT   t   v   v   kSI I  k BI A 3  

(38)

Such formulation extends the anisotropic structure of pseudo-turbulent fluctuations to the whole of Bubble-Induced-Fluctuations which would be in incoherence with the anisotropic structure of Bubble-Induced-Turbulence in bubbles rising in quiescent liquid (du Cluseau et al., 2019). Table-1 summarizes the characteristic time scales related to each class of turbulence model .

-20-

Table 1: Classification of turbulence model for bubbly flow SIT/BIT decomposition criteria

Time scales criteria

No

One-time scale

No

Two-time scales

Yes

Two-time scales

characteristic time scales



t 

t 

k 

k d ;  b  B  vr

SIT d ;  b  B  vr

2.2. Reduced forms of turbulence models in homogeneous turbulent bubbly flows We analyze in this section the reduced forms of each class of turbulence models described above in homogeneous bubbly turbulence with constant shear rate S (including the uniform situation where S  0 s 1 ). The reduced form of the two first classes of models may be written as follows:

d 1  k   tS2  1    Sk  dt

d  1    dt

1 

(39)

 C1tS2  C2  S k

(40)

This reduced formulation is unable to produce a turbulence decay behind the grid in uniform homogeneous bubbly flow with the same characteristic time scale as in singlephase case. whatever the model adopted for the interfacial source term of Sk . In fact, in single-phase uniform case ( S  0 s 1 and   0 ), equations (39) and (40) are reduced as follows:

-21-

dk   dt d dt

(41)

  C2

2 k

(42)

Which allow us to deduce the following equation linking k and  :

d  dk   C2 dt k dt

(43)

This equation may be integrated after separating variables to establish that:  0

 k      k0 

C2 

(44)

Where k0 and 0 are respectively the turbulent energy generated by the grid at the flow inlet and its dissipation rate. Equation (41) may therefore be written as: dk dt

 k    0    k0 

C2 

(45)

Which is integrated to obtain:

k k0

  1    1   C  1 0 2  k0 

   t  

1 C2  1

(46)

Knowing that C2  1.92  2 , we deduce therefore from (46) that the turbulence decay behind the grid for the homogeneous uniform turbulent single-phase flow may be described by the following approximate expression :

k k0

  1    1  0  k 0 

   t    exp     0  t  

;

0 

-22-

k0 0

(47)

This relation states that turbulence in this case decays exponentially with characteristic time scale 0 . On the other hand, in the corresponding bubbly case, equations (39) and (40) may be generalized as:

Sk dk    dt 1  

(48)

S d 2   C2  dt k 1   

(49)

Taking into account the relation (8) between Sk and S , a similar developpment as in single-phase case allow us to establish the following equation : d k     dt  k 0 

   k  Sk   0   1    k 0  k 0   k 0 

C 2

(50)

If we introduce the function  Sk ,  such as:

Sk 1    k 0

2      Sk ,    0   k0 

(51)

Using the approximate value of C 2  2 , equation (49) may be written as: 2   dk 2  k  0         dt  k0    k0   k0

(52)

Which may be integrated to deduce:

k k0

  1    t  1    exp  2   0    1        1    t  1    exp  2   0     1    

(53)

Such relation leads us to describe the turbulence decay in uniform homogeneous turbulence in bubbly case by the following approximate expression:

-23-

k k0

  t   1     (Sk , ) 1  2   exp  2 (Sk , )   0    1    

(54)

This relation implies a turbulence exponentially decay for the bubbly case with the following characteristic time scale:

1 

0 2  (Sk , )

(55)

According to Lance et al (1991) experimental results, this turbulence decay should be similar to that observed in single-phase. This occurs only if  (Sk , ) 

1 which is 2

obviously not the case since it depends on  . We deduce then that the two first classes of turbulence models conceived without decomposition of the total turbulent fluctuations are structurally incapable to reproduce the turbulence decay behind the grid in uniform homogeneous bubbly turbulence with the same characteristic time scale as in single-phase case. Indeed, Figure-6 illustrates the obtained analytical solution of turbulent intensity behind the grid compared to experimental data of Lance and Bataille (1991). It shows a so important disagreement between this solution and experimental results in bubbly flows even it is in satisfactory agreement with experimental data in single-phase case.

-24-

Figure 6 : Analytical solution of turbulent intensities behind the grid in homogeneous uniform turbulent bubbly flow. Comparison with experimental data of Lance and Bataille (1991) A similar result was pointed out by Chaibina et al. (2019) when analyzing for the same turbulence case (uniform homogeneous bubbly flow) the numerical results obtained by CFD approach from several proposed turbulence models. Such limitation will be addressed with the third-class of models. In fact, the reduced form of turbulence model equations conceived according to Chahed et al. (2003) approach in homogeneous uniform bubbly turbulence ( shear rate S  0 s 1 ) may be deduced from equations (16), (17) and (20). The transport equations of the dissipative part of the turbulent energy kT and its dissipation rate T are written as:

d kT dt

  T

d T dt

  C2

(56)

T 2 kT

(57)

While the transport equation (20) of k P is reduced to the following algebraic expression:

-25-

kP



k PH



1  uR 4

2

(58)

Similarly, the reduced forms of transport equations of kSI and SI for Lopez de Bertodano et al. (1994) and Guan et al. (2015) may be deduced from (25) and (26) as:

d kSI dt

  SI

(59)

d SI dt

SI2   C2 k SI

(60)

The transport equation of kBI is reduced respectively for Lopez de Bertodano et al. (1994) and Guan et al. (2015) as:

k BI



k aBI



d  1    k BI   dt

1  uR 4

2

3 CD  vG  v 4 dB

(61) 3

 1     BI

(62)

The set of equations (56)-(57) and (59)-(60) are equivalent to reduced forms in singlephase case (equations (41) and (42)). Therefore, when affecting the single-phase inlet boundary conditions of turbulence ( k 0 and  0 ) to the corresponding ones of SIT contributions, the quantities kT and kSI issued from the two separating approaches will be identified to single-phase turbulent energy in homogeneous uniform bubbly flow situation. On the other hand, the equations (58) and (61) show that when bubbles reach the limit slip velocity in their rising movement (controlled by balance between gravity and drag force), the quantities k P and k BI will be maintained constant. We then deduce that total turbulent energy produced by both models of Thee-Equation (Bellakhal et al, 2004b) and of Lopez de Bertodano et al. (1994) will decay correctly (as observed experimentally) in homogeneous uniform bubbly turbulence in a similar way as in singlephase case. We may also note that the reduced form of transport equation (62) of k BI for Guan et al. (2015) model is unable to produce a still constant value of k BI and consequently an adequate decay of turbulence behind the grid in bubbly cases.

-26-

However, when considering the uniformly sheared bubbly flows, the models proposed by Lopez de Bertodano et al. (1994) and Guan et al. (2015) are uncapable to produce the turbulent shear decrease observed in homogeneous uniformly sheared bubbly flows experiments in comparison to the corresponding single-phase flow. These models are conceived using the additive concept of turbulent viscosity modeling of Sato et al. (1981) as explained by equation (30), turbulent shear may thus be calculated as follows:

T 

Re yz



T



 SIt S   BI t S

Re_ SIT yz

(63)

As the shear rate S was maintained constant, equal to 2.9 s-1, both in single-phase and several bubbly cases in experiments carried out by Lance et al. (1991), the first term in rhs in this relation may be identified to the turbulent shear in single-phase flow, while the second one is a supplementary turbulent shear generated by bubbles agitation. Consequently, such relation is unable to produce the decrease of turbulent shear induced by the bubbles’ presence as observed in the experiments of Lance et al. (1991). Such limitation has been overcome with the first order turbulence model Three-Equation (Bellakhal et al, 2004b) by means of appropriate formulation for turbulent viscosity (equation 22). In fact, similarly to the precedent analysis, the pseudo-turbulent shear is also zero and the turbulent shear is written according to (22) as follows:

T 

Re yz



T 

Re_ T yz

cb k P   1  c 0 k T  k2    c T    T   1   t  b 

  S   

(64)

This relation may be rewritten as:

T 

Re yz



one  phase   TRe



(65)

yz

Where  is a dimensionless coefficient given by:

cb k P c 0 k T  1  t b

1  

(66)

-27-

one  phase and  TRe  the single-phase turbulent shear given by: yz

T



one  phase Re yz

 k2    c T  S  T 

(67)

where k T and  T may be identified to the turbulent energy and its dissipation rate in single-phase case of Lance et al. (1991) experiments. Recall that the coefficient  that appears in turbulent viscosity model (22) as well as in (65) results from the reduction of second order turbulence closures proposed by Chahed et al. (2003) as detailed in Bellakhal et al. (2004a). The coefficient  involves two ratios k P to k T (the pseudoturbulent energy to the turbulent one) and  t to  b (respectively the characteristic timescales related to the turbulence and pseudo-turbulenc). This represents competitive effects related to two antagonist aspects induced by the presence of bubbles: enhancement of turbulent viscosity by bubbles’ agitation and random eddies stretching in liquid phase which reduces anisotropic structure of turbulence and as consequence the turbulent shear as well the turbulent viscosity Bellakhal et al. (2004b) and Chaibina et al. (2019). According to this separation approach, the turbulence model is structurally able to produce the characteristic behaviors of bubbly homogeneous turbulence uniform and uniformly sheared as described above. However, when applied to homogeneous bubbly flows and confronted to experimental data obtained by Lance and Bataille (1991) and Lance et al. (1991), the authors were constrained to overestimate the relative velocity (by adjusting the drag force coefficient) in order to get a satisfactory concordance with experiments. Indeed, Numerical results produce a relative velocity around 0.4 m/s while corresponding experimental data is around 0.27 m/s. This means that it is necessary to adjust the pseudo-turbulence at the inlet in order to obtain numerical results that are in concordance with experimental data. In our opinion, the overestimating of relative velocity in homogeneous bubbly turbulence induces an increasing of pseudo-turbulent energy k P to compensate the amount of turbulent energy produced by local shear in bubble’s wake. We remind that this turbulent contribution does not appear in three-equation turbulence model because turbulence in the bubble’s wakes is supposed to be in equilibrium. The turbulence

-28-

produced in the bubbles wakes is thus not explicitly calculated. The turbulence and pseudo-turbulence are thus adjusted at the inlet so that the total turbulence level is well reproduced. Another important issue is that by construction of the tree-equation turbulence model, the turbulent dissipative part contains the turbulence produced in the bubbles’ wakes. The latter is associated with the specific bubbles’ scales while transport equation of the turbulent dissipative part involves a single scale as a whole for the all turbulent dissipative contributions. Further progress in bubbly flows turbulence modeling could be achieved by developing specific closures for the turbulence produced in bubbles’ wakes especially for high void fraction bubbly flows. To this respect, the construction of five-equation turbulence model is appropriate to remedy this incoherence by modeling specific transport equations of the dissipative turbulent energy at the bubble scale as this turbulence is produced and dissipated in the bubbles wakes. 3. Five-equation turbulence model for bubbly flow proposal In five-equation turbulence model, the total turbulent fluctuations in bubbly flow is submitted to a double decomposition. It is first split into SIT and BIT components produced respectively by gradient of mean velocity in continuous phase and bubble’s presence. The BIT is then split into two components corresponding to non-dissipative pseudo-turbulent part (as introduced by Chahed et al. (2003)) and dissipative turbulent part produced by local shear in bubble’s wake as explained schematically in Figure-7.

-29-

Figure 7 : Schematic decomposition of total turbulent fluctuations proposed for FiveEquation turbulence model Form dynamical point of view, the two components of the BIT contribution in homogeneous bubbly turbulence may be interpreted as the energy transfer rates by means of drag force and added mass force. Hence, the total turbulent energy k is split into three contributions: two dissipative Shear-Induced-Turbulence contributions kT and kS produced by shear in the continuous phase respectively by the mean liquid velocity and by the bubbles’ wakes, and a non-dissipative pseudo-turbulent contribution kP , so that we may write:

k  3.1.

k T  kS  k P

(68)

Transport equations modeling

The dissipative turbulent part of kinetic energy k T and its dissipation rate  T as well as the pseudo-turbulent part kP are still modeled following Chahed et al. (2003). We thus keep with the two characteristic time-scales of turbulent fluctuations (produced by

-30-

gradient of mean velocity in liquid phase) and bubbles agitation  t and  b involved in the three-equation turbulence model (equation 19). In addition, we introduce a supplementary time-scale characteristic of the dissipative turbulent fluctuations produced by shear in bubble’s wake. It is expressed by:

s 

kS S

(69)

The transport equations of k T and  T are therefore written as follows:

 1    k T    . 1    v k T     t  . 1    t k T  s k S  b k P   k T   1     PT  T   1    T    . 1    v T     t  . 1    t k T  s k S  b k P   T  

1    t

C1 PT

 C2 T 

(70)

(71)

Where the production term PT of kT is expressed using turbulent viscosity formulation as:

PT 

 

 

t   t  v .   v   v   

(72)

The transport equation of pseudo-turbulent part kP is also written as:

 1    k P    . 1    v k P     t 





 k   . 1    t k T  s k S  b k P   k P    PH   . v k PH   t 

(73)

The dissipative turbulent part of the Bubble Induced turbulence k S assumed to be produced by the drag force and dissipated in the bubbles’ wake. Such an assumption leads us to model its transport equation as follows:

-31-

 1    k S    . 1    v k S     t   . 1    t k T  s k S  b k P   k S   1     PS  S 

(74)

where the production term PS is equal to the turbulent energy production by drag force in the relative motion of the bubbles. It may be written as follows:

PS 

3 3 CD  vG  v 4 dB

The dissipation rate S

(75) is modeled as in single-phase flow following the

phenomenological approach proposed by Launder et al. (1975). Its transport equation is modeled as follows:

 1    S    . 1    v S     t   . 1    t k T  s k S  b k P   S  

1    s

C1S PS

 C2S S 

(76)

Where C1S and C2S are new model constants. Nonetheless, when considering the bubbly homogeneous turbulence situation where a production- dissipation equilibrium is reached in bubbles’ wake, equations (74) and (76) are reduced to the following forms:

0 

PS  S

(77)

0 

C1S PS  C2S S

(78)

We deduce therefore that to adequately take into account this equilibrium observed in experiments, we should impose:

C1S  C2S  CS

(79)

We pose CS  1. The transport equation of S is thus written as follows:

-32-

 1    S    . 1    v S     t   . 1    t k T  s k S  b k P   S  

1    s

 PS

 S 

(80)

The diffusion terms appearing in the five transport equations of model (70), (71), (73), (74) and (76) are modeled using a gradient law (following Bellakhal et al. (2004b)) with diffusion coefficient including besides turbulent diffusion effect (expressed by t kT ) two supplementary diffusion effects corresponding to diffusion induced by bubbles agitation (expressed by bkP ) as well as a diffusion effect induced by the turbulence in the bubbles’ wake (expressed by skS ). 3.2.

Turbulent viscosity

The new time scale s involved in the present study is characteristic of turbulence produced by local shear in the bubbles’ wake. It may be estimated using the following relation deduced from the following production-dissipation balance:

3 CD  u 3R 4 dB

 S 

kS s

(81)

On the other hand, several works (Lopez de Bertodano et al. (1994); Risso and Ellingsen (2002) ; Chahed et al. (2004) ; Riboux et al. (2010) ; Mrabtini et al. (2017) ; …) suggest that kS may be scaled by:

k S   u 2R

(82)

Which, taking account of (80), implies that :

S

 4  dB     3CD  u R

(83)

A similar result was established by Lopez de Bertodano et al. (1994). We deduce therefore that the two characteristic time scales S and b are of the same order of magnitude. Consequently, the Bubble-Induced-Turbulence (BIT) as described in the five-equation turbulence model ( k BI  k S  k P ) may be roughly characterized by the time scale b .

-33-

Based on this analysis, we propose to generalize the formulation of the turbulent viscosity developed for the three-equation model (equation (22)) in order to adapt its expression to the new decomposition used in the five-equation turbulence model in the form:

t



c

2 T

k T

1

cb  k P  k S    c 0  k T   1  t b

(84)

It should be observed that equation (84) is similar to the previous formulation of the turbulent viscosity (24) resulting from the reduction of second order turbulence closure except for the fact that pseudo-turbulence k P (in the ratio k P to k T ) is replaced by the total Bubble Induced turbulent energy k BI  k S  k P 3.3.

Reynolds stress tensor

Different studies carried out on dispersed bubbly both experimental (Amoura et al. (2010) ; Risso and Ellingsen (2002) ; Riboux et al. (2010)) and numerical (Lahey et al. (2005) ; Du Cluzeau et al. (2019)) have pointed out the anisotropic structure of BIT contribution. Studies of bubbles swarm rising in quiescent water or dispersed bubbles rising in homogeneous turbulent flow show that longitudinal turbulent fluctuations of BIT are more important than transversal ones. Analysing rising of millimetric bubbles in water ( d B  2.1 mm ; Vslip  0.32 m / s ), Riboux et al (2010), have reported that the ratio between the longitudinal and transversal velocities fluctuations is equal to 3/2; a very higher amount compared to the turbulence structure of the pseudo-turbulence provided by the theoretical solution of Bisheuvel and Van Vijingaarden (1984). Moreover Risso (2018) evoked much higher values for liquid flow over a random array of fixed spheres. Based on Riboux et al (2010) result, we assume that for this category of bubbles (millimitric bubbles), the anisotropic structure of Reynolds stress tensor associated to BIT may be described, by:

T  T 

BIT Re zz BIT Re xx



T  T 

BIT Re zz BIT Re yy



9 4

(85)

-34-

Recall that the z-axis is the upward vertical direction. The Bubble Induced turbulent energy may thus be expressed as:

k BI

 kS  k P



1  BIT TRe 2 





zz



BIT  2 TRe



yy

 

(86)

We deduce therefore the following relations:

k BI



17 BIT TRe 18





k BI



17 BIT TRe 8





(87)

zz

xx



17 BIT TRe 8





(88)

yy

On the other hand, according to the second decomposition adopted in the present modelling approach, we have: BIT P S  TRe  TRe TRe

(89)

P Similarly to the previous spherical tensor A defined to describe the anisotropy of TRe

(equation (37)), we introduce here the two spherical tensors D and C describing S BIT respectively the anisotropic structures of TRe and TRe . Regarding relations (86), (87)

and (88), we can write:

0   8 17 0   0  C   0 8 17  0 0 18 17  

(90)

Remarking that:

 kP  kS  C

 kP A  kS D

(91)

We thus deduce that the turbulence produced in the wake of millimitric bubbles has an anisotropic structure which may be described by the spherical tensor D given by the following expression:

-35-

k  D  C   P  C  A   kS 

 8 11 k P    17 85 k S  0     0  

0 8 11 k P  17 85 k S 0

    0   18 22 k P    17 85 k S  0

(92)

In more general terms, the Reynolds stress tensor may be calculated using the following anisotropic relation:

 

t 2   T S P    t   v   v   k T I  kS D  k P A + TRe + TRe TRe  TRe 3  

(93)

However, this description is based on a particular outcome relating to the ratio between the longitudinal and transversal velocities fluctuations (equation (85)) obtained for bubbles of diameter d B  2.1 mm rising in water at rest with velocity Vslip  0.32 m / s (Riboux et al., 2010). Further experimental and numerical works are needed to provide a general description the structure of the Bubble Induced Turbulence; especially concerning the turbulence that is produced and dissipated in the bubbles’ wake. 4. Results and interpretations 4.1.

Numerical simulation of homogeneous turbulent bubbly flows

The Five-Equation turbulence model developed here supposes that we have to model the transport equations of dissipative turbulent energies k T (70), k S (74), as well as their dissipation rates T (71), S (76), and the transport equation of the non-dissipative pseudo-turbulent energy k P (73). The turbulent viscosity is calculated using formulation (84). The Reynolds stress tensor is then deduced using relation (93). In order to be applied and validated in homogeneous turbulent bubbly flows, this model has been implemented in the CFD code ANSYS CFX. The numerical implementation is similar to that detailed in Chaibina et al. (2019) for Three-equation turbulence model.

-36-

Geometry and meshing The computational domain conceived in this work is a rectangular vertical channel corresponding to the test section of Lance and Bataille (1991) and Lance et al. (1991) experiences. It was meshed using structured non uniform cartesian mesh refined near the inlet transverse section. Boundary conditions Symmetry boundary conditions are applied to the lateral walls of the domain. On the other hand, an outlet boundary condition is set at the outflow section where an atmospheric pressure is imposed and zero gradient for the remaining flow variables. About inlet conditions, the same profile uniformly sheared with constant shear rate S is set for mean velocity of two phases (including the uniform situation where S  0 s 1 ). This means that dispersed phase is assumed to be injected with no slip velocity. For its volume fraction, it is imposed uniform through the whole inlet section. Inlet turbulent conditions are also imposed uniform through inlet section. Numerical processing We start at first step by simulating single-phase case. Turbulent energy and its dissipation rate are adjusted at the inlet for standard k-ε model to reproduce experimental data of Lance and bataille (1991) and Lance et al. (1991). According to the decomposition adopted in five-equation turbulence model, the Shear Induced turbulence is assimilated tor the equivalent single-phase flows so that we may consider that kT and T at the inlet section take the same values as in the equivalent single-phase flows. The pseudo-turbulent energy kP is imposed equal zero at inlet. The same holds for the turbulent part produced by local shear in bubble’s wake kS and its dissipation rate S . However, regarding the modeled transport equation (76) of S and expression (83) of characteristic time scale S , we affect a very little value to kS (

k Sinlet 

k inlet T ) in order to avoid division by zero at first computing iterations. The whole 103

-37-

of simulations in bubbly cases have been carried out with a bubble characteristic diameter equal to 5 mm as it was measured by Lance and bataille (1991) and Lance et al. (1991). Homogeneous uniform turbulence Figure-8 shows simulation results obtained from Five-Equation turbulence model and measured turbulent intensities in homogeneous uniform turbulent flow in single-phase flow and in three bubbly flows with void fractions 0.5%, 1% and 2%. This figure shows a satisfactory agreement between numerical and experimental data. The Five-Equation turbulence model succeeds in reproducing the adequate decay of the turbulence behind the grid as observed in experiments. In these simulations, the relative velocity is equal to that measured by Lance and Bataille (1991) ( u R  0.27 m / s ), which means that there was no need to overestimate relative velocity by calibrating drag force coefficient CD to get satisfactory agreement with experimental data when using five-equation turbulence model as was the case with the three-equation model ( u R  0.4 m / s ).

Figure 8 : Decay of turbulent intensities given by Five-Equation turbulence model in homogeneous uniform turbulent bubbly flow. Confrontation with experimental data of Lance and Bataille (1991)

-38-

Furthermore, in order to highlight the improvements obtained by adopting the double separating approach , the same homogeneous uniform bubbly turbulence situations have been simulated using three-equation model with identical numerical conditions (with same drag force coefficient CD ). Figure-9 illustrates a comparison between the numerical results in this case and experimental data.

Figure 9 : Decay of turbulent intensities in homogeneous uniform turbulent bubbly flow obtained with three-equation model using the same numerical conditions. Comparison to experimental data of Lance and Bataille (1991) This figure shows a significant gap between numerical results obtained with threeequation model and experimental data. This may be explained by the inability of only the pseudo-turbulent component to produce the increasing of turbulent energy induced by bubble’s presence in bubbly flows.

-39-

Figure 10 : Longitudinal evolution of turbulent energy in homogeneous uniform turbulent bubbly flow case (α = 1%). Comparison of numerical results to experimental data of Lance and Bataille (1991) Figure-10 presents a curve superposing of total turbulent energy k (compared with experimental data) with its three contributions ( k t , k s and k p ) in homogeneous uniform turbulent bubbly flow with void fraction equal to 1%. This figure shows a same decay of

k t as in single-phase case as well constant values of k s and k p quickly reached behind the grid. This means that turbulent component k t in this situation is equivalent to the turbulent energy in single-phase flow while the constant value of dissipative contribution

k s results from production-dissipation balance in the bubbles’ wake. The non-dissipative contribution k p reaches constant values when bubbles reach their limit slip velocity. In summary, we deduce that the three-equation turbulence models yield result as good as those of the five equations turbulence model in homogeneous bubble turbulence. However, this is only possible once the total turbulent energy at the inlet (including the Bubble-Induced-Turbulence) is suitably adjusted. Identifying the Shear-InducedTurbulent energy amount to the corresponding turbulent energy in single-phase case, the Bubble-Induced-Turbulent energy amount at inlet should be set arbitrarily or adequately produced by adjusting the interfacial forces. These adjustments make the three-equation

-40-

model non-predictive of the initial phase of bubble acceleration. By decomposing bubbleinduced turbulence into a dissipative-turbulent contribution and a non-dissipative pseudoturbulent contribution for which specific transport equations are modelled, the fiveequation model succeeds in obtaining the turbulent kinetic energy levels observed in homogeneous bubble turbulence experiments without the need for prior adjustments at the inlet. With this decomposition, the homogeneous turbulent component is identified to single-phase flow turbulence with the same velocity of the liquid; and the bubble-induced turbulence is expressed as the sum of the turbulent and pseudo-turbulent contributions calculated with conventional added mass force coefficient and the same values of the relative velocity observed in the experiments. It should be noted on the other hand that the fluctuating field structure provided by the theoretical potential flow solution of Biesheuvel and Van Wijngaarden (1984) which is adopted in the three-equation turbulence models closure, is inconsistent with the structure of the bubble-induced turbulence produced by Direct Numerical Simulation of the rising bubbles in a homogenous swarm Risso (2018). In comparison to the theoretical solution obtained in potential flow, the level of bubble-induced turbulent energy produced by Direct Numerical Simulation is much higher and the turbulent field is much more anisotropic. Based on Direct Numerical Simulation of bubble-induced turbulence, Du Cluzeau et al. (2019) indicated that it is necessary to proceed with the decomposition of the turbulence-induced turbulence associated with the fluctuations of potential origin and turbulence produced in the bubble’s wakes in order to develop suitable turbulence closure of turbulence in bubbly flows. This recommendation is in line with the conclusions of previous works (Chahed et al. 2003 and 2004) that mentioned the need to better specify the contribution of turbulence in the bubble’s wake which, because of hypothesis of the production-dissipation equilibrium in the bubbles wake, does not appear explicitly in the transport equations of the three-equation model. Other modelling issues relating to the hydrodynamic interactions associated with the deformation of the bubbles or those, which are likely to occur in bubbly flows with relatively high void fractions are also evoked (Chahed et al. 2004, Lopez de Bertodano 2017). 4.2.

Homogeneous uniformly sheared turbulence

Figure-11, Figure-12 and Figure-13 present numerical simulation and experimental data of longitudinal evolution of turbulent energy and turbulent shear respectively in

-41-

single-phase and bubbly cases (   1% and   1.4% ) for homogeneous uniformly sheared turbulence where the same shear rate was kept (S = 2.9 s-1).

Figure 11: Longitudinal evolution of total turbulent energy and turbulent shear rate, in homogeneous uniformly sheared turbulent single-phase flow. Comparison to experimental data of Lance et al. (1991)

-42-

Figure 12: Longitudinal evolution of total turbulent energy and turbulent shear rate, in homogeneous uniformly sheared turbulent bubbly flow (α = 1%). Comparison to experimental data of Lance et al. (1991)

Figure 13: Longitudinal evolution of total turbulent energy and turbulent shear rate, in homogeneous uniformly sheared turbulent bubbly flow (α = 1%). Comparison to experimental data of Lance et al. (1991)

-43-

Figures (11), (12) and (13) show good agreement between simulations and experimental results. These figures indicate that the five-equation turbulence model is in a position to reproduce the decrease of turbulent shear induced by bubble’s presence. Similar results have been previously achieved using three-equation turbulence model (Bellakhel et al. (2004b); Chaibina et al., 2019). This important result obtained by both three-equation and five equation turbulence models was made possible due to the turbulent viscosity formulation deduced from second order turbulence closures. It should however be noted, that with three-equation turbulence model, it was necessary to overestimate relative velocity to agree with experimental data. The implementation of the Five-equation turbulence model allows to obtain similar results with the value of the relative velocity observed in the experiments. Similarly to the above analysis for homogeneous uniform turbulence case, Figure-14 presents for the bubbly uniformly sheared turbulence case with void fraction 1 % a curve superposing of total turbulent energy k (compared to experimental data) with its different contributions : k T , k S and k P .

-44-

Figure 14: Longitudinal evolution of turbulent energy in bubbly homogeneous uniformly sheared turbulent flow case with S = 2.9 s-1 and α = 1% . Comparison of numerical results to experimental data of Lance and Bataille (1991)

Figure (14) shows the longitudinal evolution of k t . An equilibrium for both k s and k p is rapidly reached. The pseudo-turbulent energy amount seems to be so low compared to the one of the turbulent energy produced in the bubbles’ wake. Furthermore, figure-15 illustrates the longitudinal evolution for this case of the two dissipation rates T and S . It shows equally a longitudinal evolution of T similar to that observed in single-phase case, whereas S points out an evolution with constant value resulting from balance between production and dissipation established in bubble’s wake.

-45-

Figure 15: Longitudinal evolution of the two dissipation rates T and S in bubbly homogeneous uniformly sheared turbulent flow case with S = 2.9 s-1 and α = 1%

As regards the turbulent viscosity, we present in Figure-16 a comparison between numerical produced values for this variable in both single-phase case and bubbly case with α = 1%.

-46-

Figure 16: Longitudinal evolution of the turbulent viscosity in both single-phase case and bubbly case with α = 1% for homogeneous uniformly sheared turbulence Figure (16) shows a decreasing of the turbulent viscosity induced by bubbles agitation which makes possible to reproduce the attenuation of the turbulent shear as observed in the experiments of Lance et al. (1991). This highlights the well founded of the turbulent viscosity formulation. In fact, the expression of the turbulent viscosity (equation 84) involves two ratios: (k P  k S ) to k T and  t to  b . This represents competitive effects related to two antagonist aspects induced by presence of bubbles : enhancement of turbulent viscosity by bubbles agitation and random eddies stretching in liquid phase which reduces anisotropic structure of turbulence and consequently the turbulent shear and the turbulent viscosity. Such competition aspect is illustrated by Figure-17 , which shows the different ratios between the involved time scales.

-47-

Figure 17 : Longitudinal evolution of the ratios of times scales

 t t , and s in s b b

bubbly homogeneous uniformly sheared turbulent flow case with S = 2.9 s-1 and α = 1%

Table-2 shows the values of the ratio

 TRe  yz k

obtained from numerical results in

homogeneous uniformly sheared flows where equilibrium production-dissipation is reached. These values are compared to those deduced from experiments, from the analytical expression established in Bellakhal et al (2004a) and from numerical results obtained by three-equation turbulence model in Chaibina et al. (2019).

-48-

Table 2 :  TRe  yz k ration in homogeneous uniformly sheared bubbly flows

Experiment (Lance et al., 1991) Analytical expression (Bellakhal et al., 2004) Three-equation model (Chaibina et al., 2019) Current work

0%

1%

1.4%

-0.3

-0.21

-0.15

-0.3

-0.2

-0.17

-0.3

-0.2

-0.15

-0.3

-0.17

-0.143

5. Conclusion We have proposed in this work a five-equation turbulence model for bubbly flows where total turbulent energy is split into three contributions : two turbulent dissipative quantities ( k t and k s ) produced respectively by gradient of mean velocity in liquid phase and local shear in bubble’s wake, and a non-dissipative component k p attributed to fluctuations of potential flow around bubbles. The proposed model has been validated in homogeneous bubbly turbulent flow uniform and uniformly sheared. Numerical results succeeded to reach satisfactory concordance with experimental data of Lance and Bataille (1991) and Lance et al. (1991) with correct numerical relative velocity as observed in experiments. As perspectives of this work, we suggest that it should be so useful to investigate inhomogeneous turbulent bubbly flows situations using the proposed five-equation model. The bubbly turbulent free shear flows should be so appropriate at this advancing stage to be tackled as the difficulties related to boundaries presence are still avoided in those cases. In previous works where three-equation turbulence model was validated in bubbly shear layer (Chahed et al., 2003 ; Ayed et al., 2007) and bubbly wake flow (Chahed et al., 2003). As turbulence produced by local shear in bubble’s wake is hidden

-49-

in the adopted turbulence modeling approach, it was taken account for by adjusting the pseudo-turbulence energy at the inlet as it was explained by authors. In fact, experimental data exhibit a significant gap between bubbly case measures and the corresponding ones obtained in single-phase case. Such difference was taken into account by adjusting pseudo-turbulence at inlet. We expect that five-equation model will be able to avoid such empirical adjustment when introducing the local wake shear produced turbulence.

Acknowledgements This work was supported by the Modelling Laboratory in Hydraulics and Environment at the National Engineering School of Tunis, Tunisia.

References Amoura, Z., Roig, V., Risso, F., Billet, A. M., 2010. Attenuation of the wake of a sphere in an intense incident turbulence with large length scales. Physics of Fluids, Vol (22) n°5, pp. 055105-055105.

Aouadi, A., Bellakhal, G., Chahed, J., 2019. Effect of the turbulent contribution of the interfacial momentum transfer on the bubbles dynamics and void fraction distribution in vertical bubbly jets. International Journal of Multiphase Flow, vol. 114, pp. 82-97.

Ayed, H., Chahed, J., Roig., V., 2007. Hydrodynamics and Mass Transfer in a Turbulent Buoyant Bubbly Shear Layer. AIChE Journal, Vol (53) N°11, pp 2742-2753

Balachandar S, Eaton JK., 2010. Turbulent dispersed multiphase flow. Annu. Rev. Fluid Mech., vol(42), pp 111–33.

Bellakhal G., Chahed J., Masbernat L., 2004a. Analysis of the turbulence structure in homogeneous shear bubbly Flow using a turbulent viscosity model. Journal of Turbulence, Vol (5), 036

Bellakhal, G., Chahed, J., Masbernat, L., 2004b. K-Omega turbulence model for bubbly flows. 5th International Conference on Multiphase Flow, ICMF’04, Yokohama, Japan, May 30–June 4, 2004, Paper No. 319

-50-

Biesheuvel, A., Van Wijngaarden, L., 1984. Two-phase flow equation for a dilute dispersion of gas bubbles in liquid, J. Fluid Mech., Vol (148), pp. 301-318.

Chahed J., Colin C., Masbernat L., 2002. Turbulence and phase distribution in bubbly pipe flow under micro-gravity condition. Journal of Fluids Engineering, vol. 124(4), pp. 951-956.

Chahed J., Roig V., Masbernat L., 2003. Eulerian-Eulerian two-fluid model for turbulent gaz-liquid bubbly flows. International Journal of Multiphase Flow, vol. 29-1, pp. 23-49.

Chahed, J., Bellakhal, G., Masbernat, L., 2004. Turbulence and pseudo-turbulence modelling in low and high void fraction bubbly flows. 5th International Conference on Multiphase Flow, ICMF’04, Yokohama, Japan, May 30–June 4, 2004, Paper No. 318

Chahed, J., Aouadi, A., Rezig, M., Bellakhal, G., 2016. The phenomenon of bubbles negative relative velocity in vertical bubbly jets. Journal of Fluids Engineering, vol. 138(12), pp. 121301.

Chaibina, F., Bellakhal, G., Chahed, J., 2019. First and second order turbulence closures applied to homogeneous turbulent bubbly flows. Journal of Applied Fluid Mechanics, Vol (12), No. 6, pp. 1813-1823

Colin, C., Fabre, J., Kamp, A., 2012. Turbulent bubbly flow in pipe under gravity and microgravity conditions. Journal of Fluid Mechanics, vol. 711, pp 469-515. ISSN 00221120

Colombet D, Legendre D, Risso F, Cockx A, Guiraud P., 2015. Dynamics and mass transfer of rising bubbles in a homogenous swarm at large gas volume fraction. J. Fluid Mech. Vol (763), pp. 254-85.

Colombo, M., Fairweather, M., 2015. Multiphase turbulence in bubbly flows: RANS simulations. Int. J. Multiph. Flow, Vol (77), pp 222–243.

Drew, D.A., 1983. Mathematical modeling of two-phase flow. Annu. Rev. Fluid Mech. 15, 261–291.

-51-

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., Vol (866), pp 132–168.

Guan, X., Li, Z., Wang, L., Li, X., Cheng, Y., 2015. A dual-scale turbulence model for gas–liquid bubbly flows. Chinese J. Chem. Eng. 23, 1737–1745.

Hassan, Y. A., Blanchat, T. K., Seely, J.R., Canaan, R. E., 1992. Simultaneous velocity measurements of both components of a two-phase flow using particle image velocimetry. Int. J. Multiphase Flow, Vol (18), No. 3, pp. 371-395.

Ishii, M., 1975. Thermo-fluid dynamic theory of two-phase flow. Eyrolls, Collection de la Direction des Etudes et Recherche d'Electricité de France.

Ishii, M., Hibiki, T., 2011. Thermo-fluid Dynamics of Two-phase Flow. Second ed. Springer.

Kataoka, I., Serizawa, A., 1989. Basic equations of turbulence in gas-liquid two-phase flow. Int. J. Multiph. Flow, Vol (15), pp. 843-855.

Kato, H., Iwashina, T., Miyanaga, M., Yamaguchi, H., 1999. Effect of microbubbles on the structure of turbulence in a turbulent boundary layer. J Mar Sci Technol, Vol (4), pp 155-162.

Kawamura, T., Kodama, K., 2002. Numerical simulation method to resolve interactions between bubbles and turbulence. International Journal of Heat and Fluid Flow, Vol (23) 5, pp 627-638.

Krepper, E., Lucas, D., Prasser, H.-M., 2005. On the modelling of bubbly flow in vertical pipes. Nuclear Engineering Design, Vol (235), pp 597–611.

Krepper, E., Rzehak, R., 2011. CFD for subcooled flow boiling: Simulation of DEBORA experiments. Nuclear Engineering and Design, Vol (241), pp 3851–3866.

-52-

Lahey, R. T. Jr., 2005. The simulation of multidimensional multiphase flows. Nuclear Engineering and Design, Vol (235), pp 1043–1060.

Lance, M., Bataille, J., 1991. Turbulence in the liquid phase of a uniform bubbly air water flow. J. Fluid Mech., Vol (222), pp. 95-118. Lance, M., Marié, J.L., Bataille, J., 1991. Homogeneous turbulence in bubbly flows. Journal of Fluids Engineering, Vol (113), pp. 295-300.

Launder, B.E., Reece, G.J., Rodi, W., 1975. Progress in the development of a Reynolds stress turbulence closure. J. Fluid mech., Vol (68), part 3, pp. 537-566.

Lee, S., Lahey Jr., R.T., Jones, O.C., 1989. The prediction of two-phase turbulence and phase distribution using a k-ε model. Jpn. J. Multiphase Flow, Vol (3), pp 335–368.

Liu T.J., Bankoff S.G., 1990a. Structure of air-water bubbly flow in a vertical pipe: ILiquid mean velocity and turbulence measurements. Int. J. Heat and Mass Transfer., vol. 36 (4), pp 1049-1060.

Liu T.J., Bankoff S.G., 1990b. Structure of air-water bubbly flow in a vertical pipe: IIVoid fraction, bubble velocity and bubble size distribution. Int. J. Heat and Mass Transfer., vol. 36 (4), pp 1061-1072.

Lopez de Bertodano M., Lee S.J., Lahey R.T., Drew D.A., 1990. The prediction of twophase turbulence and phase distribution using a Reynolds stress model. Journal of Fluids Engineering, vol. 112, pp. 107-113

Lopez de Bertonado M., Lee S.J., Lahey R.T., Jones. O. C., 1994. Development of a k-ε model for bubbly two-phase flow. Journal of Fluids Engineering, vol. 116, pp. 128-134.

Lopez de Bertonado M., Fullmer, W., Clausse, A., Ransom, V. H., 2016. Two-Fluid Model Stability, Simulation and Chaos, Springer, Book, ISBN 978-3-319-44967-8, ISBN (eBook) 978-3-319-44968-5

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.

-53-

Magolan, B., Lubchenko, N., Baglietto, E., 2019. A quantitative and generalized assessment of bubble induced turbulence models for gas-liquid systems. Chemical Engineering Science: X 2, 100009.

Moursali, E., Marié, J.,L., Bataille, J., 1995. An upward turbulent bubbly layer along a vertical flat plate. Int. J. Multiphase Flow, Vol (21), pp. 107-117.

Morel, C., 1995. An Order of Magnitude Analysis of the Two-Phase K-ε Model. International Journal of Fluid Mechanics Research, Vol (22), issue 3-4, pp 21-44.

Morel, C., 1997. Turbulence modeling and first numerical simulation in turbulent twophase flows. Technical Report. CEA Grenoble, France.

Mudde, R. F., Gravity-driven bubbly flows. Annu. Rev. Fluid Mech. Vol (37), pp. 393423.

Mrabtini, H. A., Bellakhal, G., Chahed, J., 2017. Analysis of the homogeneous turbulence structure in uniformly sheared bubbly flow. Nuclear Engineering and Design, Vol (320), pp 112-122.

Politano, M., Carrica, P., Converti, J., 2003. A model for turbulent polydisperse twophase flow in vertical channels. Int. J. Multiph. Flow, Vol (29), pp 1153–1182.

Rezig, M., Bellakhal, G., Chahed, J., 2017. Phase distribution in dispersed liquid-liquid flow in vertical pipe: mean and turbulent contributions of interfacial force. AIChE J., Vol (63), pp 4214–4223.

Riboux G, Risso F, Legendre D. 2010. Experimental characterization of the agitation generated by bubbles rising at high Reynolds number. J. Fluid Mech., Vol (643), pp 509–39.

Risso, F., Ellingsen, K., 2002. Velocity fluctuations in a homogeneous dilute dispersion of high-Reynolds-number rising bubbles. J. Fluid Mech., Vol (453), pp 395-410.

-54-

Risso, F., 2018. Agitation, mixing, and transfers induced by bubbles. Annu. Rev. Fluid Mech., Vol (50), pp 25–48.

Rzehak, R., Krepper, E., 2013. Cfd modeling of bubble-induced turbulence. Int. J. Multiph. Flow, Vol (55), pp 138–155.

Rzehak, R., Kriebitzsch, S., 2015. Multiphase CFD-simulation of bubbly pipe flow: A code comparison. International Journal of Multiphase Flow, 68, 135–152.

Rzehak, R., Ziegenhein, T., Kriebitzsch, S., Krepper, E., Lucas, D., 2017. Unifed modeling of bubbly flows in pipes, bubble columns, and airlift columns. Chem. Eng. Sci. 157, 147–158.

Sato, Y., Sadatomi, M., Sekoguchi, K., 1981. Momentum and heat transfer in twophase bubble flow-i. Theory. Int. J. Multiph. Flow, Vol(7), pp 167–177.

Serizawa A., Kataoka I., Michiyoshi I., 1975a. Turbulence structure of air-water bubbly flow : I Measuring techniques, Int. J. Multiphase Flow, vol. 2, pp. 221-233

Serizawa A., Kataoka I., Michiyoshi I., 1975b. Turbulence structure of air-water bubbly flow : II- Local properties", Int. J. Multiphase Flow, vol. 2, pp. 235-246

Serizawa A., Kataoka I., Michiyoshi I., 1975c. Turbulence structure of air-water bubbly flow : III- Transport properties. Int. J. Multiphase Flow, vol. 2, pp. 247-259

Simiano, M., Lakehal, D., Lance, M., Yadigaroglu, G., 2009. Turbulent transport mechanisms in oscillating bubble plumes. J. Fluid Mech. Vol (633), pp 191-231.

Simiano, M., Lakehal, D., 2012. Turbulent exchange mechanisms in bubble plumes. Int. Journal of Multiphase Flow. Vol (47), pp 141-149.

Simonin, O., 1990. Eulerian formulation for particle dispersion in turbulent two-phase flows. Sommerfeld M, Wennerberg D, editors. Proceedings of the Fifth Workshop on Two-phase Flow Predictions. Erlangen, Germany, pp 156-166.

-55-

Simonin, O., Viollet, P.L., 1990a. Modelling of turbulent two-phase jets loaded with discrete particles. Hewitt FG, et al., editors. Phase-Interface Phenomena in Multiphase Flows. Washington, DC: Hemisphere Publishing Corporation, pp 259–269.

Simonin, O., Viollet, P.L., 1990b. Prediction of an oxygen droplet pulverization in a compressible subsonic co-flowing hydrogen flow. Symposium on Numerical Methods for Multiphase Flows, Toronto, Canada, Jan. 4–7.

Tabib, M. V., Roy, S. A., Joshi, J. B., 2008. CFD simulation of bubble column—An analysis of interphase forces and turbulence models. Chemical Engineering Journal, Vol (139), pp 589–614

Tomiyama, A., Kataoka, I., I. Zun, I., Sakaguchi, T., 1998. Drag coefficients of single bubbles under normal and micro gravity conditions. JSME international journal. Series B, Fluids and thermal engineering, Vol (41) 2, pp. 472-479.

Vaidheeswaran, A., Hibiki, T., 2017. Bubble-induced turbulence modeling for vertical bubbly flows. Int. J. Heat Mass Transfer, Vol (115), pp 741–752.

Wang S.K, Lahey Jr R.T, Jones Jr O.C., 1987. Three-dimensional turbulence structure and phase distribution measurements in bubbly two phase flows. Int. J. Multiphase Flow., vol. 13, pp. 327-343.

Ziegenhein T, Rzehak R., Lucas, D., 2015. Transient simulation for large scale flow in bubble columns Chemical Engineering Science, Vol (122), pp 1-13.

Ziegenhein T, Rzehak R, Ma T, Lucas D. 2017. Towards a unified approach for modelling uniform and non-uniform bubbly flows. Can. J. Chem. Eng., Vol (95), pp 170 – 79.

-56-

Highlights : • Turbulence in bubbly flows is governed by different but strongly coupled dynamics • Total turbulent in bubbly flow is first split into SIT and BIT components • BIT part is secondly split in pseudo-turbulence and shear wake produced turbulence • Five-equation turbulence model for bubbly flow is developed

-57-

Authors Contribution Statement

Ghazi Bellakhal Conceptualization, Methodology, Writing Original Draft Fethia Chaibina Conceptualization, Methodology, Software Jamel Chahed Supervision, Conceptualization, Methodology, Writing- Reviewing and Editing

-58-

-59-