Wall function treatment for bubbly boundary layers at low void fractions

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 κ-ε 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 κ-ε model to account for the effects of bubbles. The numerical implementation of the near wall functions is made through a finite elements code.


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 interpenetrating 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. (1990) 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 (τ 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. (1981).
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 andHassan 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 α, τ 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).Cienc (2018) 90 (1) WALL FUNCTION TREATMENT FOR BUBBLY BOUNDARY LAYERS AT LOW VOID FRACTIONS

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 singlephase 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 (1980).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 (1981).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 wallskewed.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. (1987).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é (1987).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.3void fractions.
To overcome some of the experimental difficulties previously found in other works, Moursali et al. (1995) 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 (1997).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 singlephase 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.Cienc (2018) 90 (1) WALL FUNCTION TREATMENT FOR BUBBLY BOUNDARY LAYERS AT LOW VOID FRACTIONS

BASIC EQUATIONS
Transport equations for bubbly two-phase flows have been introduced by Lopez de Bertodano et al. (1994), 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 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 α 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 where, according to the previous remarks, with ν s t = eddy viscosity due to shear; ν b t = eddy viscosity due to bubble agitation.Different approaches can be used to specify the eddy viscosities ν s t and ν b t .The shear induced viscosity can be modeled through the mixing-length hypothesis, that is, 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 κ-ε model, so that and the turbulent kinetic energy, κ, and the dissipation rate of turbulent kinetic energy, ε, are given by the two transport equations where all the c s and σ s are model constants.Typical values of the empirical constants for single phase flows are shown in Table I.
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 where κ l is an empirical constant (= 1.2 to 1.4), α 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 α 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 An Acad Bras Cienc (2018) 90 (1) WALL FUNCTION TREATMENT FOR BUBBLY BOUNDARY LAYERS AT LOW VOID FRACTIONS 579 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 η = ord u 2 * and ord η = ord (1/(u * R)), 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 (1987), 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 κ-ε model.In fact, if the κ 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 where η 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 η-limit (see Cruz and Silva Freire 1998) with ord(η) = ord(u * 2 ) onto Eqs.(1) to (9), we get These are the equations of motion, in the sense of Kaplun, that hold in the fully turbulent region.Wall functions for U, κ 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 η-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.

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, 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 α l ∼ = 1 − α g max .
A first integration of Eq. ( 14) yields where u * (= τ w /ρ l ) is the friction velocity.
A second integration furnishes where the wall variables are defined as U + l = U l /u * and y + = yu * /ν, and Note that as α 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 ν s t = κyu * .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 β different from that defined by Eq. ( 17).

DIFFERENTIAL MODEL
Solutions for κ 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 solution of Eqs. ( 11), ( 12) and ( 13) also implies that In the limit as β 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, σ ε 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 β −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), β has been defined according to The bubbles slip velocity is evaluated from (Ishii and Zuber, 1979), where σ is the surface tension and ∆ρ 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 δ + = 11, where δ + denotes the thickness of the viscous sublayer.

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 5x10 −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. (2007) and Fontoura Rodrigues et al. (2013).
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 δ stands for the external boundary layer velocity.Figure 2a     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 δ = 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 δ + 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 α G MAX
The sensitivity of the present formulation to changes in α g max is expressed in terms of four effects, the inclusion of parameter β into the following equations:  20), 4. the equation for c ε 1 , Eq.( 24).
Figure 3 shows the effects of the peak void fraction, α g max , on predictions of the friction velocity.The models were tested with α g max as high as 20%.
An increase in void fraction results in an increase of u * that clearly cannot be captured by the singlephase 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 α 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%.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 α 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.

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 + .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 κ + 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 κ-ε 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 β .
shows a sketch of the flow configuration, An Acad Bras Cienc (2018) 90 (1) WALL FUNCTION TREATMENT FOR BUBBLY BOUNDARY LAYERS AT LOW VOID FRACTIONS 583 adapted fromMarié et al. (1997).Figs.2b-cpresent the details of the computational domain as well as the mesh refinement considered.

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