Wall function treatment for bubbly boundary layers at low void fractions

Abstract

The present work investigates the role of different treatments of the lower boundary condition on the numerical prediction of bubbly flows. Two different wall function formulations are tested against experimental data obtained for bubbly boundary layers: (i) a new analytical solution derived through asymptotic techniques and (ii) the previous formulation of Troshko and Hassan (IJHMT, 44, 871-875, 2001a). A modified k-e model is used to close the averaged Navier-Stokes equations together with the hypothesis that turbulence can be modelled by a linear superposition of bubble and shear induced eddy viscosities. The work shows, in particular, how four corrections must the implemented in the standard single-phase k-e model to account for the effects of bubbles. The numerical implementation of the near wall functions is made through a finite elements code.

Keywords
bubbly flow; law of the wall; k-e; pipe flow; turbulence

INTRODUCTION

The significant advances experienced over the last forty years on turbulence modeling have meant that the prospect of achieving a three-dimensional representation of multiphase flows has evolved into procedures of practical application. The large increase on computer power and improvement on numerical techniques have, in particular, succeeded in providing a framework onto which sophisticated predictive codes can be constructed.

For disperse flows, a common solution strategy is to consider the phases separately using the inter-penetrating continua concept. In this case, the macroscopic effects of the interaction between phases must be correctly addressed through constitutive equations. Many examples can be found in literature, notably the works of ^{Drew and Lahey (1982)}, Lopez de Bertodano et al. (¹⁹⁹⁰) and Troshko and Hassan (^2001a).

Irrespective of the chosen solution strategy, the treatment of the wall boundary condition always presents great difficulties. The very steep changes in flow properties in the near wall region result in an extremely thin characteristic layer that, as expected, demands the use of exceptionally fine meshes in the numerical computation of flows.

Some early studies have proposed to separate this thin layer into two typical regions: a viscous (inner) sublayer where turbulent and laminar stresses are of comparable magnitude and a defect (outer) layer where the turbulent stresses provoke a small perturbation to the external flow solution. This natural division of the flow suggests the use of local analytical solutions to bridge the viscous dominated region, in what is normally referred to as the wall function method.

Unfortunately, the structure of a near wall turbulent bubbly boundary layer has long been recognized as depending strongly on the complex interactions that occur between the discrete and continuous phases. In addition to turbulence caused by shear, bubble agitation greatly contributes to the transfer of momentum in the liquid phase. Pipe flow experiments on the distribution of void fraction have clearly established the existence of profiles changing from a saddle to a parabolic shape as the void fraction increases at constant water velocity (^{Serizawa et al. 1975}, Sato et al. 1981a, b). The saddle distribution is typical of bubbly flow whereas the parabolic distribution corresponds to slug flow.

Hence, in bubbly flows there exists a peak in void fraction, normally at a distance from the wall of about one mean bubble radius. The effect of the bubbles in the liquid phase in the immediate vicinity of the wall is to provoke an increase in the mean velocity gradient and in the turbulence intensity. The overall asymptotic structure of a two-phase turbulent boundary layer is, however, found to be in accordance with single-phase flow theory (^{Marié et al. 1997}).

Early numerical simulations of two-phase flows resorted to procedures developed specifically for single-phase flows. To account for void fraction effects on flow solution, more realistic wall functions were proposed by Marié et al. (1997) and Troshko and Hassan (^2001b). Wall functions for the mean liquid velocity profile (U), the turbulent kinetic energy (κ) and the dissipation rate (ϵ) were developed in terms of the local void fraction (α), wall shear stress (${\tau }_w$ ), slip velocity ($U_R$) and two empirical parameters (𝜘, ${𝜘}_l$) for horizontal flow. Results were then validated against the vertical flow data of Marié et al. (1997), Sato et al. (^1981b) and Nakoryakov et al. (¹⁹⁸¹).

The present work revisits some of the previous analysis on bubbly boundary layers at low void fractions, proposing in its own term new local solutions for U, κ and ϵ. Here, we follow a different route from other authors (Marié et al. 1997, Troshko and Hassan 2001a), by using intermediate limits to derive equations that are subsequently integrated to yield local solutions. In the process, a new equation for the transport of ϵ is also derived. In fact, one of the standard single-phase constants in the ϵ-equation, $C_{{ϵ}_1}$, is shown to be a function of α, ${\tau }_w$, $U_R$, 𝜘 and ${𝜘}_l$ for two-phase flows. The new solutions for U, κ and ϵ correct the expressions of Troshko and Hassan (2001a).

The numerical implementation of wall functions in bubbly flows is also discussed here. An investigation of the parametric sensitivity of the proposed wall function on the void fraction is of special concern. The new boundary conditions are implemented in the computational fluid dynamics program TURBO2D. Predictions are compared with the bubbly flow data of Marié et al. (1997) and with results obtained with the model of Troshko and Hassan (2001b).

SHORT REVIEW ON SOME PREVIOUS CONTRIBUTIONS TO THE PROBLEM

For two-phase bubbly flows, early studies tended to consider that very close to the wall the classical single-phase logarithmic law persisted. This belief was partially motivated by the small diameter of the pipes that were used in the experiments. In fact, the placement of standard local probes near to the wall was a problem that seriously compromised detailed studies.

The existence of peak values of void fraction in the wall region of a bubbly flow was reported by Zun (¹⁹⁸⁰). This phenomenon was then associated to a transverse bubble migration from the core region to the wall for certain flow regimes. A possible physical mechanism was proposed by considering the effects of Magnus and Zhukovski forces. Thus, the non-equilibrium bubble transverse migration was treated by combining the bubble dispersion and the circulation of liquid around the bubble provoked by the liquid velocity gradient.

A discrepancy between the velocity logarithmic distribution occurring in a single-phase flow and the velocity profile in a bubbly flow was observed by Sato et al. (1981a, b). By considering the turbulent structure of the liquid phase to consist of two components, one dependent only on the shear stress of the liquid phase and the other on the additional turbulence caused by the bubbles, an eddy diffusivity model was advanced capable of accounting for the transfer of momentum and heat in a bubbly flow. The theory permitted the prediction of mean liquid velocity profiles and frictional pressure gradients for a given void fraction distribution. A comparison with experiments was performed for model validation. The overall agreement between data and theory was found to be good.

An empirical correlation for the turbulent viscosity in a two-phase flow was also advanced by van der Welle (¹⁹⁸¹). The theoretical arguments were again based on the assumption that the turbulent field could be divided into two independent components: one due to the momentum exchange of the liquid phase, the other due to the movement of the dispersed phase. Tests of the correlation against data from various authors were performed, showing a standard deviation of 22%.

^{Beyerlein et al. (1985)} reported that the bubble void fraction profile in a vertical upward flow is wall-skewed. To account for this flow behaviour, the authors incorporated into the equations of motion a lateral force due to the relative velocity of the two phases and the eddy diffusivity of the liquid. A good agreement was noted between the theoretical predictions and the experimental data.

The three dimensional turbulence structure and the phase distribution in bubbly two-phase flows were investigated by Wang et al. (¹⁹⁸⁷). Using both single- and three-sensor hot-film anemometer probes, measurements of local void fraction, liquid velocity and Reynolds stresses were made. For upward flows, authors found that bubbles tend to migrate to the wall whereas for downward flows bubbles migrate to the pipe centerline. These two distinct tendencies result respectively in a void fraction peaking at the wall and in a “coring” phenomena that can be predicted by considering the turbulence structure of the continuous phase and the lateral lift force acting on the bubbles. Measurements of the Reynolds stress components showed nearly flat profiles in the core region, but an anisotropic structure near the wall. The presence of bubbles in a liquid flow is normally observed to increase the level of turbulence. Wang et al. (1987), however, observed that for the higher flow rates bubbles suppressed turbulence.

The modelling of the skin-friction and of the heat transfer in bubbly flows in pipes was considered by Marié (¹⁹⁸⁷). Two main assumptions were used by the author: the persistence of the logarithmic region and the existence of similarity between the modifications caused by the bubbles and those that would be caused by a grid in a single phase flow. Therefore, the resulting expressions were supposed to be valid just for low gas concentrations. The model was shown to work well for up to 0.2-0.3 void fractions.

To overcome some of the experimental difficulties previously found in other works, Moursali et al. (¹⁹⁹⁵) and Marié et al. (1997) turned their attention to flows developing over a vertical, smooth, flat plate in the presence of small bubbles. In Moursali et al. (1995) graphs of void fraction distribution, wall shear stress and liquid mean velocity profiles were presented for different mean bubble diameter. An important result was the recognition that a significant fraction of the bubbles was deflected toward the wall depending on their size. This migration together with a marked deceleration of the bubbles in the near wall region proved to be the two main mechanisms responsible for the void peaking phenomenon. The presence of the dispersed phase was found to increase the skin-friction coefficient. This increase is reflected on a modification of the classical law of the wall, and on a depression of the wake. Marié et al. (1997) studied the changes in the logarithmic law of the wall resulting from the presence of small bubbles. Then, through simple analytical considerations and dimensional analysis, a modified law was proposed. The wall friction calculated on the basis of the new law was shown to fit well the experimental data. The authors also presented longitudinal turbulence intensity profiles and showed that turbulence is increased by two main mechanisms: a modification of the wall production and the creation of pseudo-turbulence in the external layer. The mixing length calculated from the data was compared with some other models proposed in literature. Their conclusion was that the single-phase logarithmic velocity profile is definitively altered when subjected to the presence of the bubbles.

The effects of bubble size and of two-phase flow rates on the wall shear stress were investigated experimentally by Liu (¹⁹⁹⁷). Using a flush-mounted hot film sensor, the time varying fluctuations of the wall shear stress were measured in a air-water bubbly flow in a vertical channel. The reported experiments were unique in the sense that a special bubble generator was capable of decoupling the bubble size effect from the inlet conditions. Thus, the experiments were carried out under various fixed gas and liquid fluxes, with only the bubble size being a variable. The data show that the wall shear stress is strongly influenced by the wall structure of the flow, while both the liquid phase velocity and the wall concentrated bubbles are the dominant parameters on both the magnitude and the fluctuation intensity of the wall shear stress in the regime of bubbly flow. The findings were compared with the data of other authors as well as with other models for the prediction of the wall shear stress.

Troshko and Hassan (2001a) developed a new formulation for the law of the wall on horizontal flows by considering the total liquid turbulent stress to result from the summation of the bubble induced local stress and the shear induced stress. Both stress components were estimated through the Boussinesq turbulent viscosity approximation. The non-linear interaction between the shear and the bubble induced turbulence fields was accounted for by a proportionality coefficient. The authors conclude through a numerical simulation that the new law performs better than the classical single-phase law.

^{Colombo and Fairweather (2015)} stated that, at the present time, no generally accepted formulation that accounts for bubble induced turbulence has yet emerged. Some works considered similarity with single-phase shear-induced turbulence, while others (^{Rzehak and Krepper, 2013}) proposed the use of a mixed time scale with the velocity scale derived from the liquid turbulent kinetic energy and the length scale equal to the bubble diameter, which should be known or estimated in advance. Bubble induced turbulence modelling is a complex subject and, despite the number of works dedicated to the subject, the literature still lacks a thorough understanding of the governing mechanics.

BASIC EQUATIONS

Transport equations for bubbly two-phase flows have been introduced by Lopez de Bertodano et al. (¹⁹⁹⁴), and, more recently, by Troshko and Hassan (2001b). For an account on the full set of equations the reader is referred to the original references.

Here, only the x-component of the liquid momentum equation is considered. For an incompressible, isothermal, two-phase, turbulent boundary layer in a Cartesian coordinate system, Troshko and Hassan (2001a) write

$\begin{array}{cc}\hfill \frac{\partial \left({\rho }_l{\alpha }_l{U}_l{U}_l\right)}{\partial x}+\frac{\partial \left({\rho }_l{\alpha }_l{U}_l{V}_l\right)}{\partial y}=& -{\alpha }_l\frac{\partial P}{\partial x}+F_x+{\rho }_l{\alpha }_l{g}_x+\frac{\partial }{\partial x}\left({\rho }_l{\alpha }_l\left(2{\nu }_l\frac{\partial U_l}{\partial x}-⟨u^2⟩\right)\right)\hfill \\ & +\frac{\partial }{\partial y}\left[{\rho }_l{\alpha }_l\left({\nu }_l\left(\frac{\partial U_l}{\partial y}+\frac{\partial V_l}{\partial x}\right)-⟨uv⟩\right)\right],\hfill \end{array}$

where $U_l$ and $V_l$ are the longitudinal and transversal components of the mean liquid velocity, $⟨u^2⟩$ and $⟨uv⟩$ are the Reynolds stress components, $F_x$ and $g_x$ are the interfacial force density and gravity projections and ${\alpha }_l$ is the local liquid fraction.

Equation (1) can only be solved provided the interfacial forces and the Reynolds stress components are modeled and appropriate boundary conditions are furnished.

Concerning the interfacial forces, the simplest possible approach is to consider the gas bubbles as mere voidages, so that no transfer of momentum occurs in the gas phase and, therefore, the flow dynamics is entirely determined by the liquid phase. Turbulence in the liquid phase is decomposed into contributions due to shear and to bubble agitation. This latter assumption is considered valid for void fraction levels below 10% (^{Lance and Bataille, 1991}).

The linear superposition of shear and bubble induced turbulence effects means that in the boundary layer the turbulent shear stress can be written as

$-⟨uv⟩={\nu }_t\left(\frac{\partial U_l}{\partial y}\right),$

where, according to the previous remarks,

${\nu }_t={\nu }_t^s+{\nu }_t^b,$

with ${\nu }_t^s$ = eddy viscosity due to shear; ${\nu }_t^b$ = eddy viscosity due to bubble agitation.

Different approaches can be used to specify the eddy viscosities ${\nu }_t^s$ and ${\nu }_t^b$. The shear induced viscosity can be modeled through the mixing-length hypothesis, that is,

${\nu }_t^s=l_m^2\frac{d{U}_l}{dy},$

where $l_m$($=𝜘y$) is the mixing length and $𝜘$ ($=$ 0.4) is the von Kármán’s constant.

An alternative is to consider the $\kappa$-$ϵ$ model, so that

${\nu }_t^s=c_{\nu }\frac{{\kappa }^2}{ϵ},$

and the turbulent kinetic energy, $\kappa$, and the dissipation rate of turbulent kinetic energy, $ϵ$, are given by the two transport equations

$\frac{D\kappa }{Dt}=P-ϵ+\frac{\partial }{\partial y}\left(\frac{{\nu }_t}{{\sigma }_{\kappa }}\frac{\partial \kappa }{\partial y}\right),$

$\frac{Dϵ}{Dt}=c_{{ϵ}_1}\frac{ϵ}{\kappa }P-c_{{ϵ}_2}\frac{{ϵ}^2}{\kappa }+\frac{\partial }{\partial y}\left(\frac{{\nu }_t}{{\sigma }_{ϵ}}\frac{\partial ϵ}{\partial y}\right),$

$P=-⟨uv⟩\left(\frac{\partial U_l}{\partial y}\right),$

where all the $c\prime$s and $\sigma \prime$s are model constants. Typical values of the empirical constants for single phase flows are shown in Table I.

Model constants for single phase flows.

$c_{\nu }$	$c_{{ϵ}_1}$	$c_{{ϵ}_2}$	${\sigma }_{\kappa }$	${\sigma }_{ϵ}$
0.09	1.44	1.92	1.0	1.30

The bubble induced turbulence is suggested by Sato et al. (^1981a) to be modelled by accounting for the drift phenomena of liquid particles due to the relative motion of gas bubbles. The result is

${\nu }_t^b={𝜘}_l{\alpha }_{g_{max}}y{U}_R,$

where ${𝜘}_l$ is an empirical constant (= 1.2 to 1.4), ${\alpha }_{g_{max}}$ the peak of gas void fraction and $U_R$ the mean relative velocity of the bubbles.

Sato et al. (1981a) mention that Eq. (9) bears similarity to the kinematic viscosity of a free turbulent flow such as a wake behind a solid body. Parameter ${\alpha }_{g_{max}}$ corresponds to the probability of vortex generation by the transit of bubbles.

The system of equations (1) to (9) needs to be complemented by appropriate boundary conditions. This is normally achieved with the specification of wall functions. We shall see this next.

ASYMPTOTIC ANALYSIS

The purpose here is to find a local solution to the system defined by Eqs. (1) to (9) in regions adjoining a solid wall. In many works this has been made through ad hoc arguments. For example, a common approach is to consider the existence of a flow region where turbulent effects dominate and an equilibrium condition is established between turbulence production and dissipation. While some useful results may be inferred sometimes, procedures with a lack of rigorous foundation very often yield misleading conclusions.

Perturbation methods have proved to be a powerful tool in the determination of global solutions. They establish a mathematical formalism, with simple rules and recipes, that can be used to find global solutions from matched local solutions. For most of these simple rules, the concepts of matching and overlap do not appear explicitly. However, matching is, by its nature, a comparison of two approximations in their domain of overlap.

Here, an application of Kaplun’s intermediate variable theory to Eqs. (1) to (9) is made to find local equations defined in terms of a formal domain of validity. In fact, Kaplun limits (^{Kaplun, 1967}) have been applied to a single-phase boundary layer by ^{Cruz and Silva Freire (1998)}. The paper studies the asymptotic structure of both the velocity and the temperature turbulent boundary layers for attached and separating flows. A relevant result is the determination of two distinguished limits according with the scales $ord\eta =ord{u}_{*}^2$ and $ord\eta =ord\left(1/\right(u_{*}R\left)\right)$, where $u_{*}$ and $R$ denote the friction velocity and the Reynolds number, respectively. This result sheds new light onto the asymptotic structure derived by Sychev and Sychev (¹⁹⁸⁷), giving a more complete interpretation to the stretching scales derived there.

The details of the analysis of Cruz and Silva Freire (1998) are not repeated here. They can be obtained directly from the original reference. We just emphasize that all results are derived without appealing to any particular closure model; only asymptotic arguments need to be used. On this note, it must be said that the asymptotic structure remains the same as considered in a $\kappa$-$ϵ$ model. In fact, if the $\kappa$ and $$ equations are considered in the analysis of Cruz and Silva Freire (1998), it can be shown that the overlap domain defined by the region where turbulent effects dominate is given by the set

$D=\left\{\eta /ord\right(1/\left(u_{*}R\right))<ord(\eta )<ord({u_{*}}^2\left)\right\},$

where $\eta$ denotes the region of validity of the approximated equations.

The fully turbulent region, determined by Eq.(10), must be interpreted as the overlap domain of the inner and outer solutions.

Passing the $\eta$-limit (see Cruz and Silva Freire 1998) with $ord\left(\eta \right)$ = $ord\left({u_{*}}^2\right)$ onto Eqs. (1) to (9), we get

$0={\alpha }_l\frac{\partial }{\partial y}\left[{\nu }_t\frac{\partial U_l}{\partial y}\right],$

$0=c_{{ϵ}_1}\frac{ϵ}{\kappa }{\nu }_t{\left(\frac{\partial U_l}{\partial y}\right)}^2-c_{{ϵ}_2}\frac{{ϵ}^2}{\kappa }+\frac{\partial }{\partial y}\left(\frac{{\nu }_t}{{\sigma }_{ϵ}}\frac{\partial ϵ}{\partial y}\right).$

These are the equations of motion, in the sense of Kaplun, that hold in the fully turbulent region. Wall functions for $U$, $\kappa$ and $ϵ$ must be constructed on their basis.

Note that according to the above analysis, the local motion equations do not contain a gravity term. In fact, the formal passage of the $\eta$-limit process onto Eqs. (1) to (9) shows that gravity effects are balanced by the inertial and pressure terms in the outer layer. These effects are then transmitted to the inner layer through the turbulent stresses in the region defined by domain (10).

Thus, any solution to Eqs. (11) to (13) should apply to horizontal as well as vertical flows.

LAW OF THE WALL FOR BUBBLY FLOW

ALGEBRAIC MODEL

The simplest way to find a local solution for the velocity profile is to solve Eq. (11) together with Eq. (3), that is,

$0={\alpha }_l\frac{\partial }{\partial y}\left[({𝜘}^2{y}^2\frac{\partial U_l}{\partial y}+{𝜘}_l{\alpha }_{g_{max}}{U}_{R}y)\frac{\partial U_l}{\partial y}\right].$

To integrate Eq. (14), the void fraction profile is considered to have a rectangular step shape with a peak value in the fully turbulent region. The following approximation is then considered ${\alpha }_l\cong 1-{\alpha }_{g_{max}}$.

A first integration of Eq. (14) yields

$({𝜘}^2{y}^2\frac{\partial U_l}{\partial y}+{𝜘}_l{\alpha }_{g_{max}}{U}_{R}y)\frac{\partial U_l}{\partial y}=\frac{u_{*}^2}{1-{\alpha }_{g_{max}}},$

where $u_{*}$ (=$\sqrt{{\tau }_w/{\rho }_l}$) is the friction velocity.

A second integration furnishes

$U_l^{+}=\frac{\beta }{𝜘}ln\left(y^{+}\right)+B^{+},$

where the wall variables are defined as $U_l^{+}$ = $U_l/u_{*}$ and $y^{+}$ = $y{u}_{*}/\nu$, and

$\beta =\frac{{\varkappa }_l{\alpha }_{g_{max}}{U}_R}{2\varkappa u_{*}}(\sqrt{1+\frac{{\left(2\varkappa u_{*}\right)}^2}{{\left({\varkappa }_l{\alpha }_{g_{max}}{U}_R\right)}^2\left(1-{\alpha }_{g_{max}}\right)}}-1).$

Note that as ${\alpha }_{g_{max}}$ tends to zero Eq. (16) reduces to the classical single-phase law of the wall. Troshko and Hassan (2001a) have also derived a two-phase law of wall through similar reasonings. However, instead of considering Eq. (4) for closure of the shear induced turbulence, they have taken ${\nu }_t^s$ = $𝜘y{u}_{*}$. This latter choice comes clearly from a less primitive assumption; it is not equivalent to Prandtl’s mixing length hypothesis, Eq. (4). The objective consequence is that Troshko and Hassan (2001a) also arrive at a logarithmic equation but with a $\beta$ different from that defined by Eq. (17).

DIFFERENTIAL MODEL

Solutions for $\kappa$ and $ϵ$ may be obtained by considering Eq. (5) as the closure relation for the turbulence generated by shear. The set of equations to be solved – Eqs. (11), (12) and (13)– must necessarily incorporate Eq. (3) with Eqs. (5) and (9) in its definition so that the effects of bubble induced turbulence are accounted for.

The resulting solution is given by

$U_l=\frac{\beta u_{*}}{𝜘}ln\left(y\right)+B,$

$\kappa =\frac{\beta u_{*}^2}{\sqrt{c_{\nu }}},$

$ϵ=\frac{\beta u_{*}^3}{𝜘y},$

where $\beta$ is defined by Eq. (17)

These equations have been written in dimensional form so that they can be easily compared with the results of Troshko and Hassan (2001a). In addition to the differences provoked by the definition of $\beta$, Eq. (19) is much distinct from Troshko and Hassan (2001a), who write $\kappa$ = $(u_{*}{)}^2/\sqrt{c_{\nu }}$.

In non-dimensional form, Eqs. (18) to (20) can be cast as

$U_l^{+}=\frac{\beta }{𝜘}ln\left(y^{+}\right)+B^{+},$

${\kappa }^{+}=\frac{\beta }{\sqrt{c_{\nu }}},$

${ϵ}^{+}=\frac{\beta }{𝜘y^{+}},$

where ${\kappa }^{+}$ = $\kappa /(u_{*}{)}^2$, ${ϵ}^{+}$ = $ϵ/\left(\right(u_{*}{)}^4/\nu )$ and clearly Eq. (21) is identical to Eq. (16).

The solution of Eqs. (11), (12) and (13) also implies that

$c_{{ϵ}_1}=c_{{ϵ}_2}-\frac{1}{{\sigma }_{ϵ}}\frac{{𝜘}^2}{\sqrt{c_{\nu }}}\frac{1}{\beta }.$

In the limit as $\beta$ tends to unit, Eq. (24) reduces to the single-flow equation for $c_{{ϵ}_1}$.

Previous studies have show that free turbulent shear flows are very sensitive to changes in $c_{{ϵ}_1}$ and $c_{{ϵ}_2}$. Variations of 10% in their values might result in changes of about 40% in the growth rate of a shear layer. In general, ${\sigma }_{ϵ}$ is fixed so that adequate adjustments on $c_{{ϵ}_1}$ are made to the flows of interest.

In a two-phase flow, Eq. (24) shows that $c_{{ϵ}_1}$ must be corrected to account for the action of the bubbles. Furthermore, Eq. (24) shows that this correction must be made in terms of ${\beta }^{-1}$.

NUMERICAL SIMULATIONS

The present near wall formulation and the theory of Troshko and Hassan (2001b) will be validated against the data of Marié et al. (1997). These authors have studied experimentally the wall region of a turbulent boundary layer developing over a vertical flat plate. In addition to mean velocity data, void fraction, wall shear stress and longitudinal turbulent intensity profiles were reported. To the best of our knowledge, the data set of Marié et al. (1997) constitutes the best account of bubbly boundary layer flow that can be found in literature.

Regarding the theory of Troshko and Hassan (2001b), $\beta$ has been defined according to

${\beta }_{TH}={\left[\right(1-{\alpha }_{g_{max}}\left)(1+{\alpha }_{g_{max}}\frac{{𝜘}_l}{𝜘}\frac{U_R}{u_{*}})\right]}^{-1}$

The bubbles slip velocity is evaluated from (^{Ishii and Zuber, 1979}),

$U_R=\left[4g\sigma \Delta \rho /{\rho }_l^2{]}^{1/4}\right(1-{\alpha }_{g_{max}}{)}^{3/4},$

where $\sigma$ is the surface tension and $\Delta \rho$ is the density difference of the phases.

The coefficient ${𝜘}_l$ was determined by Troshko and Hassan (2001b) by a direct fitting of Eq. (18) to the experimental data of Marié et al. (1997). Troshko and Hassan (2001b) also determined $B^{+}$ from a fitting procedure to make sure that ${\delta }^{+}=11$, where ${\delta }^{+}$ denotes the thickness of the viscous sublayer.

The distinct behaviour of $\beta$ (Eq. 17) and ${\beta }_{TH}$ (Eq. 15) is shown in Fig. 1.

Figure 1
Behaviour of β according to Eqs. (17) and (25).

COMPUTATIONS

The numerical simulations were performed with code TURBO2D, which is a two-dimensional code based on the finite elements method. The governing equations are discretized in space through triangular finite elements, defined by linear interpolation functions. The compatibility conditions between pressure and velocity is preserved using two calculation meshes. The pressure field is calculated with a mesh with elements of types P1. The velocity and all other variables are calculated using a P1-isoP2 mesh, defined from the P1 mesh by dividing one segment into two. This procedure generates four P1-isoP2 elements from one P1 element.

Temporal discretization of the governing equations is made through a sequential semi-implicit finite difference algorithm. A time iteration process was used to remove the influence of the initial conditions on the final calculations, so that the simulation ends only when statistically steady results have been reached. To optimize the convergence, a temporal integration procedure used increasingly successive time steps, with typical values ranging from ${10}^{-6}$ seconds at the beginning of a run to 5x${10}^{-2}$ seconds at the end of the simulation. For the velocity and pressure fields and for the velocity boundary conditions, the convergence criteria was set to ${10}^{-8}$ and ${10}^{-7}$ respectively. For more details on the optimization algorithm and numerical schemes implemented in TURBO2D the reader should refer to Loureiro et al. (²⁰⁰⁷) and Fontoura Rodrigues et al. (²⁰¹³).

The final grid was obtained from successive mesh refinements with the condition that the output be independent on the number of nodes by a margin of less than 2%. Data post-processing was done with the aid of GnuPlot and GMSH (see ^{Geuzaine and Remacle, 2009}).

Tables II and III show the physical properties and the flow conditions used in the present validation, where $U_{\delta }$ stands for the external boundary layer velocity. Figure 2a shows a sketch of the flow configuration, adapted from Marié et al. (1997). Figs. 2b-c present the details of the computational domain as well as the mesh refinement considered.

Figure 2
Sketch of the flow configuration, simulation domain and detail of the isoP2 mesh distribution.

TABLE I Physical properties of the fluids.

${\nu }_{water}$ $\left[m^2{s}^{-1}\right]$	$g$ $\left[m{s}^{-2}\right]$	$\sigma$ $\left[N/m\right]$
1.0$\times$ ${10}^{-6}$	9.80	0.073

Reference parameters from experiments.

$U_{\delta }$ $\left[m{s}^{-1}\right]$	${\alpha }_{g_{max}}$	$u_{*}$ $\left[m{s}^{-1}\right]$
0.75	0.02	0.037
0.75	0.035	0.039
0.75	0.06	0.044

In contrast to the model of Troshko and Hassan (2001b), the present model has not been calibrated against a particular set of flow conditions. To establish a common ground for comparison, the present computations have been carried out with $B^{+}$ = 7.6 and ${𝜘}_l$ = 1.4.

The boundary conditions were taken directly from the experiments of Marié et al. (1997). The inlet conditions were defined by a uniform vertical velocity profile with $U_{\delta }$ = $0.75$ ms${-1}^{}$ and turbulence intensity of $0.33$%. Non-wall boundary conditions were specified in terms of zero gradients, that is, boundary conditions of Neumann type.

The wall boundary conditions were specified with a vertical displacement of $0.4$ mm from the wall. This resulted in typical values for ${\delta }^{+}$ of about 20 in the computational domain. This procedure bridges the turbulent region directly to the wall so that no computation of the buffer and viscous sublayers is required.

MODEL SENSITIVITY TO ${\alpha }_{g_{max}}$

The sensitivity of the present formulation to changes in ${\alpha }_{g_{max}}$ is expressed in terms of four effects, the inclusion of parameter $\beta$ into the following equations:

• the slope of the law of the wall, Eq. (18),

• the lower boundary condition for $\kappa$, Eq. (19),

• the lower boundary condition for $ϵ$, Eq. (20),

• the equation for $c_{{ϵ}_1}$, Eq.(24).

Figure 3 shows the effects of the peak void fraction, ${\alpha }_{g_{max}}$, on predictions of the friction velocity. The models were tested with ${\alpha }_{g_{max}}$ as high as 20%.

An increase in void fraction results in an increase of $u_{*}$ that clearly cannot be captured by the single-phase law of the wall. Both theories, present and Troshko and Hassan’s, exhibit the same general behaviour, agreeing with the general behavior of the data of Marié et al. (1997). However, the theory of Troshko and Hassan (2001b) consistently furnishes higher values of $u_{*}$. At the highest experimental value of ${\alpha }_{g_{max}}$, 0.06, our prediction of $u_{*}$ exceeds the measurements by 10%. The theory of Troshko and Hassan (2001b), on the other hand, exceeds that value by 20%.

Figure 3
Model sensitivity to ${\alpha }_{g_{max}}$.

EFFECTS ON COMPUTATIONAL TIME

Figure 4 compares the required computational times between the present formulation and the theory of Troshko and Hassan (2001b) for various values of ${\alpha }_{g_{max}}$. The theory of Troshko and Hassan (2001b) is more time consuming for all void fractions. The differences, however, are negligible considering that the average time difference, about 100 seconds, represents $8$% of machine usage time.

Figure 4
Computational time.

MEAN VELOCITY PROFILES

Predictions of the mean velocity profiles for three different void fractions, ${\alpha }_{g_{max}}$ = 0.02, 0.035 and 0.06, are shown in Figs. 5, 6 and 7, respectively. The three typical single-phase flow regions are also reproduced in a bubbly flow: the viscous region, the fully turbulent region and the wake region. The profiles were taken at a position $y$ = 1 m, according to the coordinate system shown in Fig. 2. As the void fraction increases, the logarithmic region becomes less step.

For the lowest void fraction, the agreement between the present formulation and the experimental data in the logarithmic region, 100 $<$ $y^{+}$ $<$ 800, is very good. The absence in the present formulation of wall damping functions imply in a poor agreement in the viscous region. The theory of Troshko and Hassan (2001b) overestimates $u_{*}$, and for this reason shifts the velocity profile to a lower level. For the higher ${\alpha }_{g_{max}}$, 0.06, we still get a better agreement from the present formulation. For this condition, the friction velocity calculated by the present formulation presents an agreement within 11.21% of the experimental data. The formulation of Troshko and Hassan (2001b), on the other hand, furnished values of $u_{*}$ that differed by 24% from the measured values.

Figure 5
Mean velocity profiles, ${\alpha }_{g_{max}}=0.02$.

Figure 6
Mean velocity profiles, ${\alpha }_{g_{max}}=0.035$.

Figure 7
Mean velocity profiles, ${\alpha }_{g_{max}}=0.06$.

MODEL SENSITIVITY TO ${𝜘}_l$

Table IV shows the predicted values of friction velocity for simulations where the value of ${𝜘}_l$ is varied from $1.2$ to $1.4$ while all other parameters are kept constant. Only cases related to the peak local void fraction of $0.06$ are shown. For both models, lower values of ${𝜘}_l$ result in lower values of friction velocity. For ${𝜘}_l$ = 1.2, the agreement between the present model and the data of Marié et al. (1997) is within $10.32$%.

Troshko and Hassan (2001b), as previously explained, used the data of Marié et al. (1997) to obtain a direct fitting for ${𝜘}_l$ and $B^{+}$. When these fittings are used in the computations, the error in estimated value of $u_{*}$ decreases to 1.46%. This comparison shows the dependence of this formulation on adjusted values of ${𝜘}_l$ and $B^{+}$.

TABLE IV Model sensitivity to ${𝜘}_{𝐥}$ for ${𝛂}_{{𝐠}_{𝐦𝐚𝐱}}=0.06$. Reference data from Marié et al. (1997): ${𝐮}_{*}=0.044\left[𝐦{𝐬}^{-1}\right]$.

${𝜘}_l$	Present work	Error(%)	Troshko and Hassan (2001b)	Error(%)
	$u_{*}$ $\left[m{s}^{-1}\right]$	$u_{{*}_{PW}}-u_{{*}_{exp}}/u_{{*}_{exp}}$	$u_{*}$ $\left[m{s}^{-1}\right]$	$u_{{*}_{TH}}-u_{{*}_{exp}}/u_{{*}_{exp}}$
1.2	0.04854	10.32%	0.0536	21.87%
1.3	0.04874	10.77%	0.0541	23.04%
1.4	0.04893	11.21%	0.0546	24.17%
Variable	–	–	0.0446	1.46%

EFFECTS OF ${\kappa }^{+}$ AND $c_{{ϵ}_1}$ ON FLOW SOLUTION

The combined effects of Eqs. (19) and (24) on predictions for the friction velocity have been investigated. In fact, corrections due to both equations have been implemented in both models, the present formulation and the model of Troshko and Hassan (2001b). All other parameters are held constant.

Estimation of $u_{*}$ through the present model with corrections given by Eqs. (19) and (24) are improved by 5.5%, as compared with the data of Marié et al. (1997). Results obtained with the model of Troshko and Hassan (2001b) are improved by 8.4%. Because these percentages lie within the error margin of the experimental data, the influence of the present modifications in ${\kappa }^{+}$ and $c_{{ϵ}_1}$ might seem small. However, the generality introduced by Eqs.(19) and (24) is important. These equations improve estimations for flows subjected to adverse pressure gradients and separation, when the calculation of $u_{*}$ becomes critical to the numerical solution.

CONCLUSIONS

The present work has performed an analytical study of the fully turbulent region of bubbly flows. In the course of the research, new expressions have been derived to provide the lower boundary condition for a $\kappa$-$ϵ$ modeling of the flow. These expressions include profiles for the mean velocity, the turbulent kinetic energy and the rate of dissipation. The solutions are observed to satisfy the set of local equations, Eqs (11) to (13), and imply that the standard single-phase constants in the $ϵ$-equation, $C_{{ϵ}_1}$, is a function of $\beta$.

TABLE V Summary of the corrections proposed to the single-phase $𝛋-𝛜$ to account for the presence of bubbles and comparison with the formulation of Troshko and Hassan (2001b).

Correction	Present work	Troshko and Hassan (2001b)
Slope of the law of the wall	$\beta$ = Eq. (17)	${\beta }_{TH}$ = Eq. (25)
Lower boundary condition for $\kappa$	$\kappa$ = Eq. (19)	No correction
Lower boundary condition for $ϵ$	Eq. (20)	Eq. (20), $\beta$ = ${\beta }_{TH}$
Equation for $c_{{ϵ}_1}$	$c_{{ϵ}_1}$ = Eq.(24)	No correction
Empirical constant ${\kappa }_l$	${\kappa }_l$ = 1.4	${\kappa }_l$ fitted from experimental data

Several differences exist between the presently proposed equations and those of Troshko and Hassan (2001a). In that reference, the eddy viscosity due to the shear effects is taken as ${\nu }_t^s$ = $𝜘u_{*}y$, implying that the local characteristic velocity is $u_{*}$. Their mean velocity solution, on the other hand, implies that a characteristic velocity of the form ${\beta }_{TH}{u}_{*}$ should be used instead. These two facts cannot be reconciled, yielding an erroneous procedure that among other things suggests a boundary condition for $\kappa$ that is not correct. Here, the correct boundary condition for $\kappa$ is introduced through Eq. (19).

In particular, the present work has shown how four corrections must the implemented in the standard single-phase $\kappa$-$ϵ$ model to account for the effects of bubbles. These corrections are summarized in Table V, which also identifies the differences with the formulation of Troshko and Hassan (2001b).

ACKNOWLEDGMENTS

In the course of the research, JBRL benefited from a CNPq Research Fellowship (Grant No 309455/2016-2) and from further financial support through Grants CNPq 458249/2014-9, FAPERJ E-26/102.212/2013 and FAPERJ E-26/010.002857/2014. APSF is grateful to the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for the award of a Research Fellowship (Grant No 305338/2014-5). The work was financially supported by CNPq through Grant No 477293/2011-5 and by the Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) through Grant E-26/102.937/2011.

References

• ¹
  BEYERLEIN SW, COSSMAN RK AND RICHTER HJ. 1985. Prediction of bubble concentration profiles in vertical turbulent two-phase flow. Int J Multiphase Flow 11: 629-641.
• ²
  COLOMBO M AND FAIRWEATHER M. 2015. Multiphase turbulence in bubbly flows: RANS simulations. Int J Multiph Flow 77: 222-243.
• ³
  CRUZ DOA AND SILVA FREIRE AP. 1998. On single limits and the asymptotic behaviour of separating turbulent boundary layers. Int J Heat Mass Transfer 41: 2097-2111.
• ⁴
  DREW D AND LAHEY RT. 1982. Phase distribution mechanism in turbulent two-phase flow in a circular pipe. J Fluid Mech 117: 91-106.
• ⁵
  FONTOURA RODRIGUES JLA, GONTIJO RG AND SOARES DV. 2013. A new algorithm for the implementation of wallfunctions in high Reynolds number simulations. J Braz Soc Mech Sci Eng 35: 391-406.
• ⁶
  GEUZAINE C AND REMACLE JF. 2009. GMSH: a three-dimensional finite element mesh generator with built-in pre- and postprocessing facilities. Int J Numerical Methods Eng 79(11): 1309-1331.
• ⁷
  ISHII M AND ZUBER N. 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChEJ 25: 843-856.
• ⁸
  KAPLUN S. 1967. Fluid mechanics and singular perturbations. Academic Press, New York, 369 p.
• ⁹
  LANCE M AND BATAILLE J. 1991. Turbulence in the liquid phase of a uniform bubbly air-water flow. J Fluid Mech 222: 95-118.
• ¹⁰
  LIU TJ. 1997. Investigation of the wall shear stress in vertical bubbly flow under different bubble size conditions. Int J Multiph Flow 23: 1085-1109.
• ¹¹
  LOPEZ DE BERTODANO M, LEE S-J, LAHEY RT AND DREW DA. 1990. The prediction of two-phase turbulence and phase distribution phenomena using a Reynolds stress model. J Fluids Eng 112: 107-113.
• ¹²
  LOPEZ DE BERTODANO M, LAHEY RT AND JONES OC. 1994. Development of a κ-ϵ model for bubbly two-phase flow. J Fluids Eng 116: 128-134.
• ¹³
  LOUREIRO JBR, SOARES DV, FONTOURA RODRIGUES JLA, PINHO FT AND SILVA FREIRE AP. 2007. Water tank and numerical model studies of flow over steep 2-D hills. Boundary-Layer Meteorol 122: 343-365.
• ¹⁴
  MARIÉ JL. 1987. Modelling of the skin friction coefficient and heat transfer in turbulent two-component bubbly flows in pipes. Int J Multiphase Flow 13: 309-325.
• ¹⁵
  MARIÉ JL, MOURSALI E AND TRAN-CONG S. 1997. Similarity law and turbulence intensity profiles in a bubbly boundary layer at low void fractions. Int J Multiphase Flow 23: 227-247.
• ¹⁶
  MOURSALI E, MARIÉ JL AND BATAILLE J. 1995. An upward turbulent bubbly boundary layer along a vertical flat plate. Int J Multiphase Flow 21: 107-117.
• ¹⁷
  NAKORYAKOV VE, KASHINSKY ON, BURDUKOV AP AND ODNORAL VP. 1981. Local characteristics of upward gas-liquid flows. Int J Multiphase Flow 7: 63-81.
• ¹⁸
  RZEHAK R AND KREPPER E. 2013. CFD modeling of bubble-induced turbulence. Int J Multiphas Flow 55: 138-155.
• ¹⁹
  SATO Y, SADATOMI M AND SEKOGUSHI K. 1981a. Momentum and heat transfer in two-phase bubble flow I: Theory. Int J Multiphase Flow 7: 167-177.
• ²⁰
  SATO Y, SADATOMI M AND SEKOGUSHI K. 1981b. Momentum and heat transfer in two-phase bubble flow II: A comparison between experimental and theoretical calculations. Int J Multiphase Flow 7: 179-190.
• ²¹
  SERIZAWA A, KATAOKA I AND MICHIYOSHI I. 1975. Turbulence structure of air-water bubbly, Parts I-III. Int J Multiphase Flow 2: 221-259.
• ²²
  SYCHEV VV. 1987. On turbulent boundary layer structure. PMM USSR 51: 462-467.
• ²³
  TROSHKO AA AND HASSAN YA. 2001a. Law of the wall for two-phase turbulent boundary layers. Int J Heat Mass Transfer 44: 871-875.
• ²⁴
  TROSHKO AA AND HASSAN YA. 2001b. A two-equation turbulence model of bubbly flows. Int J Multiphase Flow 27: 1965-2000.
• ²⁵
  VAN DER WELLE R. 1981. Turbulence viscosity in vertical adiabatic gas-liquid flow. Int J Multiphase Flow 7: 461-473.
• ²⁶
  WANG SK, LEE SJ, JONES JR OC AND LAHEY JR RT. 1987. 3-D Turbulence structure and phase distribution measurements in bubbly two-phase flows. Int J Multiphase Flow 13: 327-343.
• ²⁷
  ZUN I. 1980. The transverse migration of bubbles influenced by wall in vertical bubbly flow. Int J Heat Mass Transfer 6: 583-588.

Publication Dates

• Publication in this collection
  Jan-Mar 2018

History

• Received
  10 Feb 2017
• Accepted
  03 Oct 2017