Sensitivity of the Coefficients of the Closure Model of the Interfacial Transfer in Bubble Columns

Hamdi AYED***, Abir MOULDI*****, Khaled KHEDHER***** *Department of Civil Engineering, College of Engineering, King Khalid University, Abha – 61421, KSA, E-mail: hayed@kku.edu.sa (Corresponding author) **Higher Institute of Transport and Logistics of Sousse, University Sousse, Erriadh city BP 247 4023 Sousse, Tunisia ***Department of Industrial Engineering, College of Engineering, King Khalid University, Abha – 61421, KSA, E-mail: amouldi@kku.edu.sa ****Department of Civil Engineering, ISET, Nabeul, DGISET, Tunisia, E-mail: kkhedher@kku.edu.sa


Introduction
Bubble liquid gas reactors are used in several industrial fields (chemical engineering, process engineering, energy engineering, environmental engineering, waste-water treatment, etc.). These reactors are often designed to ensure transfer phenomena or chemical transformations between the phases in contact. The transformation of different forms of energy takes place, often with phase changes (boiling, condensation). Efficiency in bubble reactors (bubble column or airlift reactors) depends on mixing efficiency and interfacial transfer rates. The energy transfers depend essentially, on the diameters of the bubbles and their velocities, as well as the hydrodynamics of the liquid phase which ensures the dispersion of the bubbles.
The significant progress made in recent two decades in the field of Computational Fluid Mechanics (CFD) allows today a more phenomenological approach to transfer problems in dispersed multiphasic fluid systems. However, these development perspectives require a better understanding of exchanges at the interfaces to describe the transfer mechanisms and then better represent global exchanges. In general, the presence of bubble-flow interfaces introduces new scales characteristic of the movement of interfaces that can profoundly modify the structure of the fluctuating field. The modulation of the turbulence of the liquid by the bubbles also results in the modification of the turbulent transport properties at the interfaces.
The development of general computational codes for turbulent flows, based on a local approach, immediately aroused great interest among researchers in the various disciplines concerned with the transfer phenomena in these complex systems. These codes, intended to describe in a predictive way the turbulent flows, allow at the same time to access essential local information for turbulent transport and transfer phenomena (fields of average speed, turbulence, concentration etc.). CFD codes have been used by researchers to study different configurations of gas-liquid flows. Nevertheless, the local modeling of multiphase flows for the study of gas-liquid reactors is today a possible way of progress both in the understanding of fundamental mechanisms and in the construction of more general operational models. Some studies based on local two-phase modeling of reactors have thus been able to show that local approaches can feed one-dimensional models useful to industrialists [1,2].
If the utility of these codes is now recognized for the predetermination of turbulent flows in the presence of more or less homogeneous mixtures, at the current stage of their development, their capacity to describe the whole of the two-phase flows remains much less obvious. Difficulties persist in particular in industrial gas-liquid systems, which involve turbulent flows at high vacuum rates, in the presence of bubbles exhibiting relative movements with large Reynolds numbers, in very highly agitated configurations (agitated reactors), or in flows completely governed by gravity (bubble column or air lift) in which the very structure of the flow depends on the distribution of the phases. These difficulties relate to the modeling of turbulence and its effect on the phase distribution [3,4]. A large part of the problem lies in the still incomplete formulation of the closures used to describe interfacial transfers.
Several experiments carried out in bubbly column [5][6][7] show that the structure of the two-phase flow is strongly changed. Freedman and Davidson [8] studied the columns experimentally in terms of gas retention and liquid circulation. Thorat et al. [9] have studied the effect of the sparger design and gas retention. Buwa and Ranade [6][7] and Pfleger et al. [5] studied the hydrodynamics inside a bubble suspension column and the oscillation of the plume.
This reveals one of the difficulties of modeling the interfacial exchange of bubble flows. Indeed, the expression of the instantaneous force on which the modeling is based, is valid only for isolated bubbles without significant deformation. If one can admit that the validity of this hypothesis is realistic for flows with inclusions of small size and for diluted diets; on the other hand, the validity of this formulation of forces can be raised in the general case for bubble column with deformable bubbles. In most cases, these flows are the seat of strong hydrodynamic interactions that disturb the flow around the bubbles and thus modify the field of forces that the liquid exerts on them. These interactions are all the stronger as the void rates are high and some authors have proposed formulations that attempt to account for these effects.
In this paper, numerical simulations were performed using Eulerian two-fluid modeling developed by Chahed et al. [10] and Bellakhal et al. [11]. In closing the equations of the two-fluid model, the modulation of turbulence of the liquid phase due to the presence of bubbles is still a major problem. Indeed, the closure of the turbulence is based on a separation of the turbulence induced by the shearing and the turbulence induced by the bubbles for which specific transport equations are modeled. Two time scales are involved in turbulent transport in two-phase flows: a first time scale linked to the stretching of the vortices, as for the single-phase flow, and, a second time scale linked to the relative movements of the bubbles. This turbulence closure leads to the development of second and firstorder turbulence models that have succeeded in reproducing the experimental data in many turbulent basic bubble flows with low and moderate vacuum fractions (homogeneous turbulence, shear layers [10] and bounded bubble flows [12]).
In this work, the interfacial momentum transfer modelling is also presented and compared with the experimental data of Pfleger et al. [5] and Buwa et al. [6,7]. This configuration is chosen because comprehensive experimental data exists concerning gas hold-up and liquid velocity.

Eulerian two-fluid model for gas-liquid bubbly flows
The average balance equations, for stationary incompressible bubbly flows, are formulated in terms of weighted average by the characteristic function of each phase (symbol  ) with the following notations: the variables related to the gas are distinguished by the subscript G, the variables related to the liquid have no phase index, ̿ and ̿ are the average velocity and pressure fields,  is the density, gi is the gravity acceleration and indicates the gas presence rate. We consider turbulent gas-liquid bubbly flow with low solubility of the gas in the liquid. We assume that no coalescence or break-up occur, in these conditions we consider that the bubble diameter is still roughly constant.
The mass and momentum balances averaged in the liquid and the gas phase are thus given by the following relations.
In the liquid: In the gas: where: For gas-liquid bubbly flows, the acceleration and the weight of the bubbles are low in comparison with the force exerted by the liquid on the dispersed phase (ρG <<ρL). Thus, the equation of the momentum balance (4) reduces to the Eq. (6) which indicates that the interfacial force exerted by the liquid on the bubbles is zero: this force includes the pressure effect of the flow (Tchen force) and the interfacial term that represents the forces due to the disturbed flow.
where: 0 () Lij  is the liquid stress field of the "non-perturbed" flow and the interfacial momentum exchange MGi represents the contribution of the "disturbed" flow. The "non-perturbed" contribution would be present even in the absence of any material interface, it is identified here to the action of the local shear stress of the continuous phase ( (0)
The interfacial moment exchange represents the contribution of the "disturbed" flow around the bubbles to the density of the total force exerted by the liquid on the bubble. The formulation of the average density (per unit volume) is then modeled as follows: where d is the bubble diameter. The velocity of the ''nonperturbed" flow is identified to that of the liquid velocity Li u . The relative velocity is thus calculated as This model contains respectively the well-known drag (coefficient CD), lift (coefficient DL) averaged forces and added mass (coefficient DA) and last term represent a specific turbulent term which comes from the averaging of the added mass term [12].
Several correlations for drag CD, virtual mass CA and lift coefficient CL are proposed for gas-liquid systems. Table 1 gives some correlations for those coefficients. One of the objectives of this work is to compare these values with those proposed by Chahed et al. [12] and used by Ayed et al. [13].
The lift is calculated with a coefficient CL = 0.25, the added mass coefficient is CA= 0.5, and the drag force is calculated with the following drag coefficient: where: Eo is the Eötvös number  is the surface tension (0.072 Nm -1 for gas/water).  In bubbly flow, the Reynolds stress tensors is not only produced by shear stress, but it is necessary to add the bubble effect in their relative displacement which produces an additional turbulence in the liquid phase called "pseudo- where: 0 C  and b C  are coefficients depending on the turbulence and pseudo-turbulence anisotropy. The ratio of these coefficients is of order of unity.
Based on this turbulent viscosity formulation, Bellakhal et al. [11] developed a first-order closure and proposes a three-equation turbulence model developed for turbulent bubble flows, containing modeled transport equations of the turbulent energy k0, the viscous dissipation rate ɛ0 and the pseudo-turbulent energy kS are written: The constants of the turbulence model have the values currently adopted for single-phase turbulence closure. The various coefficients introduced in the turbulent closure laws are reported in Table 2.

Numerical method and boundary conditions
The numerical method used is a combination of finite difference methods -finite volumes based on a nonfractional process [14]. This method consists of solving the evolution equations in several separate steps, each corresponding to an elementary equation. The discretization of the equations is of the finite volume type, it defines a socalled velocity grid where the components of the average velocities and the turbulent quantities are calculated. In order to ensure the validity of the discretization of the mass balance equations in the two phases, the pressure and the void ratio are calculated at the center of each volume element. A new grid called "pressure grid" offset from the first is thus defined.
The calculations are carried out for a rectangular bubble column similar to that used by Pfleger et al. [5], and used by Buwa and Ranade [6,7]. It is of width W= 0.20 m, height H=1.2 m and depth 0.05 m. The water level is 0.45m (H/W = 2.25). The sparger consists of a set of 8 holes placed in a rectangular configuration for the gas injection, each hole having a diameter of 0.8 mm in a square pitch of 6 mm.
In this work, the vertical length of the computational domain is 1.2 m long; the transverse length is 0.2 m. The discrete grid has 101 longitudinal nodes and 91 transverse nodes. We made preliminary tests to ensure that the solution was independent of the grid. For simplicity, the sparger in our numerical simulations is represented by a simple rectangular area placed at the bottom of the column with dimension (18 mm x 6 mm). The geometric configuration and mesh a dimensional rectangular domain as shown in Fig. 1. All walls are treated as non-slip boundaries with standard wall function. The gas flow rate at the sparger is defined via inlet velocity type boundary condition with the gas volume fraction equal to unity. The bubble size at the gas inlet is choosing a uniform bubble size at the inlet of 5mm was taken here based on the study of Buwa et al. [15].
In the experiment, no liquid was injected with the gas (αG=1, ˛αL=0), so the inlet liquid velocity is set to zero in the calculation. The boundary condition for turbulent conditions k0, kS and ɛ0 at inlet are as follows: where: 150 v n  which proved to lead to a good prediction of the single-phase flow, I is the turbulence intensity of liquid phase assumed here to be standard 5% (it was not measured).
is the velocity of the inlet condition.

Results and discussion
According to the experiences of Buwa and Ranade [6,7], when a gas is injected into a column filled with the liquid, gas bubbles formed at purge holes rise in the vertical direction in the column having scales of different length and time. According to the experiences of Buwa and Ranade [6,7], when a gas is injected into a column filled with the liquid, gas bubbles formed at purge holes rise in the vertical direction in the column having scales of different length and time. Fig. 2 illustrate two snapshots of the oscillating bubble plume published by Buwa and Ranade [6] for two gas velocities of 0.13 m/s and 0.73 m/s. Using the 3 equation s model presented above, the time-averaged flow properties (vertical velocity of the liquid and gas void ratio) are compared with the experimental measurements for different measurement sections. The closure model is used with a coefficient of added mass equal to CA = 0.5, the coefficient of the lift is CL = 0.14 and a drag coefficient according to the coefficient of Otvos. A mean bubble diameter (5 mm) is used, as adopted by Buwa et al. [6,15].
In Figs. 3 and 4, it can be seen from the figures that the sinuous movement of the bubble plume is captured in satisfactory qualitative agreement with the experiments and the instantaneous distribution of the gas volume fraction obtained using the simulation respectively for a superficial gas velocity of 0.13 cm/s and 0.73 cm/s. a b Fig. 3 Gas volume fraction (a), and liquid velocity snapshots of meandering bubble plume (b) at a superficial gas velocity of 0.13 cm/s Figs. 5 and 6 respectively show the comparison of average vertical liquid velocity and gas retention profiles measured experimentally and predicted for a superficial gas velocity of 0.13 cm/s (for an H/W ratio of 2.25). Fig. 5 shows a comparison of the experimentally measured time-averaged liquid velocity of Buwa and Ranade [6] and the predicted results at the two liquid heights of 25 cm and 37cm from the bottom of the column. Detailed comparison between the experimental results and the numerical predictions is many perfect and the accelerations outside of the injection zone are well reproduced. Fig. 6 shows a comparison of the experimentally measured time-averaged gas fraction of Buwa and Ranade, It can indeed be verified that the distribution of the void fraction is fairly well reproduced by the model (Fig. 6). The simulation reproduces the transversal migration of the bubbles and the attenuation of the initial peak located at y=0c m, which is the memory of the accumulation of bubbles in the input. Especially for low surface gas velocities, where the plume propagation is narrow, the plume oscillation period is sufficiently large and it is difficult to acquire the experimental data for a very long time to obtain symmetrical gas retention profiles.
The quality of the prediction of the void fraction is strongly related to the model of the interfacial momentum transfer. In fact, the phase distribution in a bubbly flow is governed by the of momentum transfer between the two phases as discussed below. With the closure of the moment transfer given by equation (7), the moment equation (6) of the gas written: According to Eq. (17) the force exerted by the liquid on the bubbles comprises the turbulent correlations issued from the lift force on the bubbles. The lift coefficient CL is tested in various values 0, 0.14 and 0.25. Figs. 7 and 8 presents respectively the mean liquid velocity and the void fraction prediction to the closure laws for the interfacial momentum transfer and the measurement by Buwa and Ranade [6] at an H/W ratio of 2.25 and a superficial gas velocity of 0.13 cm/s at different positions in column. With the simulation without lift, the velocity profile of the liquid showen in  According to Eq. (17) the force exerted by the liquid on the bubbles comprises the turbulent correlations issued from the added mass force on the bubbles. The turbulent contribution of the added mass force contains the turbulent correlations of the liquid and of the gas. So we have to provide closures not only for the turbulence in the liquid but also for the turbulence of the gas. In the present contribution we have a crude model for the turbulence of the gas phase. In a homogeneous bubble-induced turbulence even at low void fraction (2%), Tchen's theory of turbulent dispersion does not hold. This is because the relation between the fluctuating motions of both phases is not only a problem of "one-way coupling" response of the bubbles to the liquid fluctuating excitation. The Lagrangian properties of the bubble-induced turbulence are moreover not known, and are undoubtedly strongly modified. Ayed et al. [13].
In this simulation, we have thus just assumed that the kinetic tensor of the dispersed phase is proportional to the Reynolds tensor, which is reasonable. But we take a constant value for the proportionality ratio, which is a crude approximation. We take for each diagonal component  ,   9 presents the void fraction distribution predicted for various turbulent contribution. The introduction of the averaged added mass force and its turbulent contribution with r=1 to 9. The effect of the turbulent contribution of the added mass is very noticeable on the prediction of the void rate. By increasing this ratio, we notice that the peak of the void rate moves from left to right and gives good results in the section x = 37 cm with a value r = 9 which presents the maximum value given by Tchen's theory.

Conclusions
A two-fluid Eulerian model for predicting the twophase flow dispersed in the bubble columns was presented. Relationships for the interfacial transfer closure model have been proposed. These relationships were obtained using recent data collected in the literature. In addition, a turbulence model based on the three k0-ks-ε0 equations for mixing the two phases has been proposed. This model is based on the decomposition of the turbulence induced by the movement of the bubbles and the turbulence in the liquid.
The effect of interfacial moment transfer modeling in the phase distribution phenomenon is discussed and the analysis is focused on the role played by turbulent contributions. Indeed, different numerical simulations are performed with different models of interfacial moment transfer and the distribution of the void fraction is analyzed as a function of the effect of turbulent terms in the interfacial moment transfer.
The objective of this study is to analyze the sensitivity of the model to the turbulent contribution of the added mass force and therefore a crude model for the turbulence of the gas phase was used. In this model, the turbulent correlations of the gas are linked to that of the liquid using an adjustable coefficient. A constant value of this coefficient fixed to the unit removes the effect of the turbulent contribution of the added mass force and strongly underestimates the gas dispersion. A constant value of 9 overestimates the gas dispersion.
These numerical simulations and the analysis of the transverse displacements of the bubbles in the ascent of air prove that turbulence plays an important role on the distribution of the bubbles not only by the term pressure, but also by the turbulent contribution of the interfacial force. Consequently, it is crucial to take into account the contribution of turbulence in the transfer of interfacial moment, in particular in the term of added mass.
However, new theoretical and numerical developments as well as new specific experiences are necessary to explore the role of other turbulent contributions in the balance of interfacial dynamics.