Large-eddy simulation of gas–liquid two-phase flow in a bubble column reactor using a modified sub-grid scale model with the consideration of bubble-eddy interaction

The Eulerian–Eulerian Large-eddy simulations (LES) of gas–liquid two-phase flow in a cylindrical bubble column reactor have been conducted. When considering the turbulent eddy viscosity in LES, apart from the well-accepted contributions from shear turbulence and bubble induced turbulence (BIT), the effect of the interaction between entrained bubbles and eddies with a similar turbulence length scale to the sub-grid scale (SGS) cannot be neglected. With the consideration of the bubble response to the eddies on the induced sub-grid stresses, a modified SGS model, which incorporates the Stokes number, St, was proposed. The results of LES clearly indicate that the use of the modified SGS model can effectively capture the transient bubbly flow in the cylindrical bubble column. The power turbulent kinetic energy spectrum obtained in LES indicates that a slope similar to Komogorov -5/3 scaling law and the -3 scaling law can still be identified for a critical frequency f=10.70 Hz.


Introduction
Large eddy simulation (LES) of bubbly flow in bubble column reactors adopts two approaches, which are Eulerian-Eulerian (E-E) and Eulerian-Lagrangian (E-L). The E-L approach adopts the way that the liquid phase momentum equations are solved under the Eulerian frame while the equations for dispersed phase are solved in a Lagrangian framework. The transportation of bubbles is tracked by integrating the bubble motion equation accounting for the forces acted by the liquid phase on the bubbles. A closure model in the liquid and gas momentum equations must be provided to account for liquid-bubble interactions. The bubble size distribution can be calculated as part of the solution at each time step and models are also required to account for break-up and coalescence of the bubbles. However, this kind of approach is quite computationally intensive, thus, it is still inhibitive for studying two-phase bubbly flows in large-scale bubble column reactors or at high void fraction system. As Eulerian-Eulerian two-fluid modelling describes two-phase mixture motion in a macro sense, the use of this approach may be preferable for industrial applications, especially for the case of a highly dispersed void fraction system such as bubble column reactors.
The use of Eulerian-Eulerian two-fluid LES modelling is more desirable because the adoption of low order turbulence models, such as the k-e and even for the Reynolds Stress Model (RSM), may not well capture the instantaneous eddy turbulence structures which will affect the bubble entrainment, breakage and coalescence. Deen et al. [1] did pioneering study on LES modelling for gas-liquid flow in a square cross-sectional bubble column. The simulation results were compared with employing the k-e model. It was concluded that the turbulent viscosity was overestimated by using the k-e model and only the low frequency unsteady flows could be validated. This is very likely attributable to the transient bubble-eddy interactions being ignored. It has been well accepted that LES can reproduce high frequency experimental data and predict the transient motion of the bubble plume, as experimentally observed. Lakehal et al. [2] pointed out that because Reynolds Averaging Navier-Stokes (RANS) models depend on time averaging, they appear to screen out the local fluctuations related to the turbulence as well as the fluctuations related to the interaction between the bubbles and surrounding eddies. However, these local fluctuations are at least partially remained in LES modelling, instead of time averaging, the spatial filtering is employed. It resolves turbulent eddies with the size larger than computational grid directly, while in sub-grid scale (SGS), the behaviour of eddies as well as the interaction between bubbles and the carrier phase are modelled [3]. Furthermore, RANS models are derived by assuming the turbulences are isotropic in different scale, though the bubbly flows are capable of high anisotropy in velocity fluctuations, especially in the direction of gravity [4]. Conversely, LES modelling only assume isotropy for unresolved scales instead of applying it to resolved scales at the same time.
Less studies have been done on the LES modelling of gas-liquid two-phase flow, which is more challenging, than RANS models. Bombardelli et al. conducted simulations of wandering motion in bubble plumes by using both k-e model and LES [5]. The results show that the wandering effect can be reasonably reproduced by employing Smagorinsky SGS model and bubble-slip model, however, the k-e model can only capture the wandering effect at the beginning of the simulation.
They also compared both modelling results and indicated that time-averaged axial velocities at different heights above the diffuser for the plane located at mid-thickness by LES still have some discrepancies between the modelling and experimental data. The LES Smagorinsky SGS model has been widely used in recent years, and it has been general practice that the turbulence dissipation introduced by the model is proportional to the Smagorinsky constant, ! . Deen et al. used constant value 0.1 for Smagorinsky constant in the SGS model to simulate the bubbly flow in a bubble column [1] [6], but the sensitivity of the simulation results regarding to ! is not investigated. It has been accepted that Smagorinsky's eddy viscosity model v " = ( # ∆) $ | ̅ |, where the model constant C % can be set to 0.1 is working in the most of single-phase shear flows. However, in subsequent fluid dynamics investigations, researchers have employed this formula directly into two-phase and three-phase flow by changing the model constant from 0.1 to 0.18 empirically, neglecting the multiphase fluid mechanisms behind them. The motivation behind this doing is very likely to consider the change in typical mixing length due to the hybrid effect of the presence of bubbles in the flow. Furthermore, the use of this type of Bossinesq's viscosity model to describe the turbulence energy dissipation in bubbly flow is inappropriate as the bubble wakes also feed the so-called bubble induced turbulence into the flow. Smith and Milelli [7] considered the contribution of bubble induced turbulence into the LES work on liquid phase in bubble plumes, and simulated the dispersed phase by using random dispersion model (RDM). The bubbly flow in a bubble column at low Reynolds number was also simulated by using LES accounting for the effect of bubble induced turbulence by Ma et al. [8]. A similar LES modelling for a cylindrical bubble column was conducted by Milelli et al. [9] using a relatively coarse cylindrical coordinate grid. They compared both constant and dynamic Smagorinsky SGS model in vertical shear layers laden with bubbles at very low overall gas hold-up, revealing that the dynamic Smagorinsky SGS model proposed by Germano et al. [10] did not have better performance than the constant Smagorinsky model. They also modified the SGS model by taking bubble induced dissipation into account and found that the new SGS model did not remarkably improve the simulation results. It has been revealed by the above researchers that the bubble-induced turbulence model introduced by Sato et al. [11] did not have much influence on the turbulent velocity fluctuations in a rectangular bubble column. However, these studies have indicated that the use of the Smagorinsky SGS model with and without considering BIT over-predicted the averaged axial velocity and gas hold up profile. The kinetic energy in sub-grid scale was obtained in a LES work of bubbly flows conducted by Ničeno et al. [12]. They derived the turbulent dispersion force as well as the contribution from bubble induced turbulence yielding the transported SGS kinetic energy. They also compared it with the simulation result employing dynamic Smagorinsky SGS model, an improvement was found, but there were still discrepancies with the experimental data. Liu and Li [13] used dynamic Smagorinsky SGS model to simulate the bubbly flow in a bubble column.
Different ratios of filter width to bubble diameter were tested to check the mesh independency and the results were compared with the published data. The authors analysed the obtained power spectral density from LES and reported that there exist two zones with slopes of -5/3 and -25/3, respectively. They claimed that the steep falling off of the slope in the power energy spectrum is due to the BIT as the result of injection of bubbles. Thus, it can be expected that a LES SGS model considering the contribution from BIT in bubbly flow may lead to the simulation results becoming better consistent with the experimental data but depending on how to appropriately consider the BIT. It should be pointed out that the turbulent viscosity model accounting for the BIT in twophase flow as proposed by Sato and Sekoguchi [14] was derived by assuming the flow about a fixed bubble as the flow about a cylinder. In reality, the bubbles will response to the turbulent eddies that entrain the bubbles. This response should exist at the sub-grid scale where the bubbles may not follow the turbulent eddy motion faithfully. When assessing this type of bubbles' dynamic response to eddies, one can consider the slip velocity between the bubbles and eddies to be influenced by the response of the bubbles to eddies. The bubbles appear to escape from the turbulent eddies where they are entrapped because of the buoyancy effect. Therefore, the instantaneous fluctuation of bubbles would always differ from it of the surrounding turbulent eddies in sub-grid scales, especially for the eddies having the similar size with the bubble diameter [15]. Secondly, in bubbly flow, the interfacial forces acting on the bubbles are strongly influenced by the interactions between the bubbles and the near eddies, and these forces have to be properly implemented in LES SGS modelling. In addition, the bubbles that rise in the bubble column will encounter those turbulent eddies generated by the shear turbulence and the preceding bubble wakes in their paths. The relative size difference between bubbles and eddies lead to different relaxation times. Therefore, the instantaneous filtered velocity of a bubble is strongly correlated to the turbulent eddy fluctuation velocity.
It should be mentioned here that bubble size distribution has a pivotal role in predicting gas-liquid interfacial area, which will further influence the prediction of the mass and heat transfer between phases. To obtain the bubble size distribution when using the E-E approach, additional equations accounting for bubble breakup and coalescence, together with bubble growth or shrinkage because of mass transfer are required. Studies on the bubble size distribution in various bubbly flows can be found in the literature. Different models derived from population balance equation were developed, such as the multiple size group model (MUSIG) [16], the interfacial area transport models [17][18][19][20][21]. When using the E-E modelling approach for practical applications [22], the dispersed bubbles are treated as a quasi-continuum with each computational grid in whole domain containing corresponding fractions of the carrier and dispersed phases. Separate momentum and continuity equations are solved together on the same grid for each phase. For LES, by applying the filtering at a filter scale (∆), the filter scale should fall into the inertial sub-range region in turbulent kinetic energy spectrum as a key criteria for a successful E-E LES as indicated by Niceno et al. [12], and thereby the scale of motion greater than ∆ can be resolved. However, it is noticeable that the bubble-induced large scale turbulent eddies with the size larger than ∆ cannot be resolved properly in LES, due to the missing information of interphase details above filter scale. In addition, bubble-induced turbulence not only affects the carrier phase liquid flow for the length scale smaller than bubbles but also has the impact on the large scale flow as reflected in the predicted turbulent kinetic energy spectrum obtained from the LES modelling. This may deteriorate the accuracy of predicting the resolved scale motion of the large eddies in LES. Thus the grid requirements may sometimes be in confliction when modelling with the E-E LES. Milelli [23] has proposed the requirement of gird size in the E-E LES and carried out a parametric investigation on different mesh resolutions and bubble diameters. It was suggested that 1.2 < ∆/d & < 1.5 (0.67 < d & /∆< 0.83) would be an optimum filter width [24]. The comparative study was conducted by Dhotre et al. [4] for ∆/dB = 1.2 and 2.5, and they found that good agreement with the experimental data can be still obtained for both grids. Niceno et al. [12] has trialed a grid refinement study on ∆/dB =1.25 and 2.5 but demonstrated similar quantitative vertical gas and liquid velocity to be shown by applying both grids. Liu and Li [13] also employed the E-E LES model with 5 different ∆/dB in their study on bubble column bubbly flows and have revealed that the adoption of the ratio of grid to bubble size ∆/dB =1.25 and 1.5 can yield the results agreed with the experimental results.
Generally speaking, the application of ∆/dB >1.0 for the E-E LES modelling is required based on the previous studies mentioned above. It should be noted that the E-E LES modelling of bubbly flows accounting for the bubble size distribution coupled with for bubble coalescence and breakage is still rarely reported, to the best knowledge of the authors.
This work aims to implement a modified Smagorinksy SGS model which accounts for the bubble response to the surrounding turbulent eddies through introducing a Stokes number into the LES simulation of a cylindrical bubble column. The mathematical modelling and numerical methods are presented in Section 2. The simulation results and related discussion will be then followed in Section 3, focusing on the effect of bubble response to the turbulent eddies in SGS scale and the BIT influence on the simulated turbulent kinetic energy power spectrum together with the correlations between the turbulent eddy structures and bubble distribution in the bubble column.
Finally, this paper will end with a conclusive summary of key findings.

Governing equation
The LES model is obtained by spatially filtering the equations of momentum. The current study employs a Eulerian-Eulerian two-fluid model with each phase being described by separate equations. For phase k, the filtered equations of mass and momentum can be expressed by ) ) + ∇ • ( ) ) ) ) ) = −∇ • ( ) ) ) − ) ∇ + ) ) + *,) where = and = represents for gas and liquid respectively. The terms on the right-hand side of Equation (2) stand for the contributions due to the stress, the pressure gradient, gravity and momentum exchange between each phase that arise from the actions from interfacial forces individually. The stress term can be defined by Equation (3), where ,--represents the effective viscosity for the continuous phase, which may be assumed to be consisted of the following terms: the molecular viscosity 1 , the turbulent viscosity / and an additional term to describe bubble induced turbulence 23 [4]. This is defined in Equation (4).
The bubble induced turbulence can be modelled based on Sato's model [11], which is given by Equation (5).
The momentum exchange term that introduces the interface forces is defined by Equation (6), where the terms including, from left to right, the interface drag force, lift force, turbulence dispersion force and virtual mass force, respectively.

Drag force
The interphase momentum transfer between continuous and dispersed phases because of the drag force contributed from both viscous shear (skin drag) and the pressure gradient (form drag) can be expressed by Equation (7), where the drag coefficient, 6 , can be expressed by using Equation (8) with regard to distorted bubbles [25], where the Eötvos number > = stands for the ratio of the buoyancy to the surface tension forces.

Lift force
Due to the radial velocity gradient in the bubble column, the lift force acts on the rising bubbles perpendicularly to the relative motion of gas and liquid phases. The correlation between slip velocity and the curl of liquid velocity is associated with the lateral lift force [26,27], which is given by where CL is the lift force coefficient and is based on the estimation according to Tomiyama et al.
where 2 is the bubble Reynolds number and = C =

Added mass force
The relative acceleration of the bubble and surrounding liquid is considered by the added mass force [29], which can be estimated by where CAM stands for the virtual mass coefficient and a constant 0.5 is used through this paper.

Turbulent dispersion force
Considering the effect of the random fluctuations of turbulent eddies, the turbulent dispersion force is considered in the current simulation. The formulation proposed by Burns et al. [30] to estimate the force is adopted, given by where CTD is the turbulent dispersion coefficient and is assumed to constant 1 in this work, νt is the turbulent kinematic viscosity and /6 represents the turbulent Schmidt number, /6 = 0.9 is adopted here.
In LES, the velocity in Equations (1) and (2) are described by Equation (13), where ) is the velocity that needs to be resolved in filtering, while b ) represents the instantaneous velocity and ) C stands for the unresolved part that requires the closure from the use of the SGS model in the LES simulations. It needs to be noted that the equations of mass and momentum are derived by time-averaging in RANS models. In LES modelling, these equations are solved by spatial filtering, hence, b ) and ) C are referred to grid scale and sub-grid scale velocity, respectively.
Following Garcia's work [31] but considering the contribution from the added mass on bubble translation, the relative velocity between the carrier fluid and the bubble can be obtained from Equation (14), where λ represents the different turbulent length scales in the range between the integral and Kolmogorov scales (L>λ>η) and CAM is the added mass coefficient. When the derivative of S equals zero at a certain , S will have a maximum, as defined by Equation (15). * ~+ (, * (-,-, ' ).
In this turbulent length scale range, ( ) K/0 can be regarded as the fluctuating velocity of the bubble. Thus, * can be expressed as, The size of the bubbles and their surrounding turbulence eddies are different, hence, bubbles will not response immediately to the flow motion of the eddies. Taking their slip velocity into account, the bubbles appear to get rid of the controlling from the eddies where they are entrapped [32]. As demonstrated in Equation (18) [15], the correlation between the fluctuating velocity of bubble and liquid in terms of the turbulent eddies with length scales in the inertia subrange can be expressed by Equation (18).
When Equation (18) is implemented into the sub-grid scale, the relationship can be defined as . In terms of liquid response time in SGS, 1,!5! = ∆/ ′ a,bcb , where ∆= k∆ d ∆ e ∆ ) l K/0 is the filter width and u a,bcb stands for the liquid velocity in local grid.
The turbulence dissipation due to the bubbles corresponds to the inertial subrange can be assumed that mainly occurs when approach to * , and the dissipation can be estimated by Equation (20).
Different values of the constants Cb have been trialled, but a value of 0.7 is employed which demonstrates good agreement with Camasara's results. 5 nnnn stands for the local gas hold-up after filtering. The total dissipation is given by Equation (21).
Employing the eddy viscosity model, the liquid-phase turbulence modified SGS viscosity can be modified as represented by Equation (22), where Cs is the Smagorinsky constant and represents the characteristic resolved strain rate tensor. l has been assumed as the filter length scale D in the range of the inertia subrange.

Numerical Modelling
In order to denmonstrate the reliability the proposed modified LES SGS model, the simulation of bubbly flow in the bubble column reactor, based on the work of Camarasa et al. [34] and Kulkarni et al. [35], have been carried out. The experimental settings are summarized in The solver of ANSYS CFX 18.0 was employed in the LES simulation. At inlet, a mass flow rate condition normal to inlet was used. In the current work, the inlet mass flow rate is given by ṁ= ? 6 ? , and the volume fraction for liquid and gas phases are specified as 1 = 0 and 5 = 1. The bubble diameter of 4 mm is adopted in this work, which is the typical value of gas-liquid bubble columns under the same pressure, superficial velocities and air-inlet distributor conditions.
At the top surface of the bubble column, a pressure-constant boundary (relative pressure of 0) is applied. No slip condition is employed for the wall. It should be noted that the bubble size will increase along the height of bubble columns and this change in bubble size is usually very small, thus leading to the negligence of such change in the numerical simulation, certainly in the currently available reported studies using Eulerian-Eulerian LES [1,3,44,12,13,58]. According to Zhao's et al. study on evolution of bubble size distribution from gas blowout in shallow water [59], the compressibility factor can be used for characterisation of bubble size changes with liquid depth.
By approximating of the cases of bubble column bubbly flows, the variation of the compressibility factor z for bubble size along a stationary liquid level height of 0.9m in the bubble column is very marginal. Also, the bubble rise-up in a stationary water tank is considered as an example.
Assuming that the water level height is 0.9 m and the ambient pressure to be 1 bar, one can roughly

Results and Discussion
The gas-liquid flow in the bubble columns have been simulated using a constant time step size of 0.001 s for resolving the temporal variations of the flow field. Since the bubbly flow in the bubble column is transient, the simulation was run for 100 s and the data was collected over the last 50 s until the turbulent flow field becomes statistically stationary. The results and findings based on the simulations are discussed as follows.

Grid independency study
In this section, the results obtained on various d & /∆ values from 0.57 to 0.9 were compared with the published experimental data by Camarasa et al. [33]. Figure 2 where n 2 is the average bubble velocity, N is the sampling number for collection of the instantaneous bubble velocity at the given radial position and Dt is the time step for the simulation.  Figure 2(b), the time-averaged gas hold-up is obtained by using the following expression: where T0 and T1 are the beginning and end time for sampling. As can be seen the figure that the 0.75 has been adopted in the present E-E LES modelling (around 1,200,000 mesh elements in total were adopted, especially in the center region of the bubble column. It should be noted here that this grid resolution is consistent with the condition used in Milelli's work [23] and it is expected that the turbulence with the scale being larger than the bubble diameter can be well resolved.

Predicted flow patterns, gas holdup and velocity distributions
Chen et al. [36] have In Figure 4, the averaged bubble and liquid velocity profiles at H/D = 6 obtained by applying both the modified SGS and standard Smagorinsky SGS models are presented. By applying the same method as evaluated using Equation (23), a quantitative comparison is made with the published experimental data obtained by Camarasa et al. [34]. It can be shown from Figure 4 that the averaged bubble axial velocities obtained from the LES with the modified SGS model are in good agreement with the experimental results [34], remarkably improved in comparison to the use of standard SGS Smagorinsky model. However, such consistency becomes poor close to the column wall. A likely reason is that the interaction between bubbles and turbulent eddies may be not well reflected in the modified SGS model as the bubble size is greater than the grid size. In terms of the grid set-up employed in the present simulation, Milelli's condition is only held in the central part of the bubble column, while 2 /∆ close to the wall region is much larger than 0.75 [27]. This causes a relative poor performance of the modified SGS model in the estimation of the eddy viscosity in the region near the wall. However, as LES imposes the requirement of 5 < V < 10, the use of Milelli's condition will violate this constraint. This remains for further investigation. Figure 5 shows the cross-sectional averaged gas hold-up variation along the bubble column height after time-averaging, which can be obtained by using Equation (24). As can be seen from Figure   5 that The time-averaged gas holdup distribution at H/D = 6 in the radial direction is presented in Figure   6. It can be obtained that the simulation is quantitatively consistent with the experimental data as reported in Camarasa's et al. work [34]. It should be noted, however, that the profile of the averaged gas hold-up by LES is over-predicted for the location between the column wall and core region. One explanation is that the inhomogeneity which causes bubble induced turbulence during

Turbulent kinetic energy contributed by shear turbulence and bubble induced turbulence characterized by LES with the modified SGS model
The liquid axial velocity-time series between 50 s and 100 s at the centre of H/D = 6, obtained by the modified SGS LES is shown in Figure 9. For comparison purpose, the experimental timedependent liquid axial velocity sampled corresponding to the same location but acquired by Kulkarni et al. [35] is also plotted in the figure. It needs to be noted here that the experimental data obtained by Kulkarni et al. [35] covering a period of 50 s has 8192 sample points while the results of the modified SGS LES have adopted 10,000 sample points. The LES simulation has reasonably recapture the transient fluctuations as exhibited by the experiment, even the amplitudes of the fluctuations predicted by the modified SGS LES are also consistent with the experimental ones.
The observed differences are likely caused by the different sampling rates and the unavoidable noise from the bubble detection in the experiment. As discussed in Kulkarni's work, since it is not guaranteed that all the bubbles pass entirely and centrally through the measurement volume, the chordal passage will cause refraction on bubbles that eventually lead to the relative higher amplitude liquid velocity [35]. Thus, the probability density function (PDF) and turbulent kinetic energy spectral analysis relating to the time-dependent liquid velocities would be able to provide the physical insight into bubble induced turbulence in the bubble column.  . Figure 10(a) also presents the one-dimensional spectrum of single phase flow in accordance with Pope's model [42]. The model is described as where L and h stand for integral and Kolmogorov length scale, respectively, and the model constants are C > = 0.49, > = 0 for a -5/3 spectral slope by default, C a = 6.78, v = 0.4, and β = 5.2 [43]. As shown in Figure 11(a), the turbulent kinetic energy spectrum predicted by the modified SGS model follows Pope's model well and the modified SGS LES model captures -5/3 scaling in the inertial subrange and recovers -3 scaling laws for the wave number greater than the typical wave number characterized by the bubble size, i.e. 2 = 2p/ 2 . It is interesting to note here that the representative bubble frequencies, estimated by 2 = | 5 − 1 |/2 2 , is 12 Hz when the bubble diameter is 4 mm [44]. In general wave equations, κ = x , where ν is the frequency of the wave, is the wavelength and is the mean liquid velocity [37]. Thus 2 can be converted to κ &K =302.80m EK . It can be observed from Figure 10 which is close to the representative bubble frequency. Lance and Bataille [45] have indicated that the eddies induced by bubble wake are very quickly dissipated by viscosity before turbulence spectral transfer can take place. Pope also stated that the directional information of the large scales is missing with the energy passing down the cascade. In the energy cascade (l < l yz ), the dominant process is composed of the energy transfer to successively smaller scales and viscous dissipation.
Here, l yz is the turbulence length scale between the anisotropic large eddies and the isotropic small eddies, which is hypothesised by Kolmogorov [46]. Thus, the input energy of bubbles will not take part in large length scales, which correspond to low wave numbers, but make a contribution towards higher wave numbers. The production of eddies with the size of the bubbles will contribute towards the dissipation in the higher frequency range as indicated by Lance and Bataille [45].
Since the bubble induced turbulence dissipation can be estimated by ( ) 0 , the drag force acting on the gas phase in the turbulent bubbly flow in the bubble column is roughly balanced by the bubble buoyancy on average and one can have the following estimation, given by This leads to E(k)~k -3 , which is also demonstrated by the existing experimental work as well as DNSs and consistent with our LES results as shown in Figure 10 [47][48][49][50][51][52][53][54][55]. Based on the modelling discussion, it can be assumed that the total turbulent dissipation rate e 1 is consisted of the turbulent dissipation rate due to shear turbulence and the dissipation due to the bubble-eddy interaction as described by Equation (27), viz.
The difference in the power spectrum for different SGS models is noted in Figure 11

Correlation between large eddy structures and local gas holdup
The present study has adopted the vorticity to characterise the large eddy behaviour in the bubble column, defined as the curl of the flow velocity field by = × .
It is expected that the large eddy structure development in the bubble column would be significantly affected by the entrained bubbles while this interaction between bubbles and turbulent (b) eddies has been accounted for in the modified SGS model. Thus, a correlation to reflect this coupling can be presented. Figure 12 shows the isosurface of 5 = 0.23 highlighted by the vorticity and | 1 | = 58 EK highlighted by air volume fraction at t = 90 s. It can be observed that the bubble volume fraction is strongly coupled with the vorticity. We propose the following spatial correlation between local gas hold-up and vorticity magnitude to characterize the variation of such correlation along the axial height of the bubble column, defined as

R(∆ )
Increment Δh (m) correlation coefficient distribution. As the combination shown in Figure 14, the occurrence of high gas hold-up is accompanied by the presence of a high spatial correlation coefficient.

Figure 14.
Combination of instantaneous gas hold-up and spatial correlation coefficient at t = 90s.

Interfacial mass transfer across bubbles using the modified SGS model
The volumetric mass transfer coefficient 1 is an important parameter to evaluate the efficiency of interfacial mass transfer. Besides, when employing the species transport equation to find the species concentration, the contribution from the interfacial mass transfer across the bubbles of the source terms can be calculated with the presence of 1 [56]. As the spherical bubble assumption was made in the present work, the interfacial area concentration a for the bubbles in the bubble column can be estimated by where n 2 is the bubble column volumetric averaged gas holdup. Thus, the coefficient of mass transfer 1 and the interfacial area concentration can be obtained individually from various models of mass transfer. Since the modified SGS model that taking the turbulence kinetic energy contribution from BIT and bubble interaction with the turbulence eddies into account, the relative velocity between gas and liquid and energy dissipation rate play significant roles in estimating the value of interfacial mass transfer coefficient 1 . The eddy cell model was reported by Lamont and Scott [57] and indicated that the very small scale of the turbulent eddies play significant roles in mass transfer and these motions lead to a sophisticate viscosity, once these small scale behaviors can be controlled and the surface renewal rate as well as the mass transfer mechanism can be then defined analytically and shown as, where 1 is liquid mass diffusivity of liquid phase, 1 is the turbulence dissipation rate. It can be figured out that 1 can be estimated based on the eddy cell model through a key parameter, turbulence dissipation rate 1 . The influences of the eddies induced by the bubble wakes and bubbles' dynamic responses to the surrounding liquid on the liquid turbulent kinetic energy spectrum were illustrated in Figure 10. Therefore, apart from the consideration of the simply shear turbulent dissipation, the effect of the interactions between bubbles and eddies is also needed to be addressed. After substituting Equation (27) into Equation (30) yields the estimation for 1 based on the eddy cell model can be expressed as.
Equation (31) shows that 1 is related to the kinetic energy integrated from the energy spectrum obtained in Section 3.3. In addition, the assumption can be reasonably made that the mass transfer between the bubbles and the eddies of the similar size or marginally larger is dominant in whole process.
As shown in Figure 15(b), a higher volumetric mass transfer coefficient seems to be more likely in the vicinity of the column wall when implementing the SGS model without considering bubble response to turbulent eddies. However, when the eddy size is slightly larger than the bubble and in the inertial subrange, the bubbles will be strongly entrained by eddies. This phenomenon is well demonstrated by Figure 15(a), where the volumetric mass transfer coefficient is more uniformly distributed compared to the SGS model without modification. Therefore, by employing the modified SGS model, the distribution behavior of mass transfer characterized by 1 inside the bubble column can be better analyzed.

Conclusions
LES Simulation of gas-liquid flow in a bubble column reactor has been carried out using the (1) It can be observed from the simulation resulting from the modified SGS model that the gas hold-up and velocity profiles demonstrate a better agreement overall with the experimental results [34,35] compared with the standard Smagorinsky SGS model, but both gas hold-up and the streamwise gas velocity are slightly over-predicted in the vicinity of the bubble column wall.
(2) The use of the modified SGS turbulence model is able to capture the detailed flow behaviour of bubbly flow in the bubble column. (3) The power turbulent kinetic energy spectrum of the axial liquid velocity indicates that the slope of classical -5/3 law can still be observed for the frequency range of f < 10.70 Hz, followed approximately by a -3 scaling law when the frequency f > 10.70 Hz, the representative bubble frequency calculated according to 2 = | 5 − 1 |/2 2 is 12Hz. This is consistent with the recent findings on the bubble induced turbulence as reported by Prakash et al. [44], Bouche et al. [47], Mendex-diaz et al. [48], Mercado et al. [49], Murai et al. [50], Riboux et al. [51,53], Bunner and Tryggvason [52], Roghair et al. [54] as well as Sugiyama et al. [55], indicating that the slope of -3 law has been also recovered by using the modified SGS LES for the bubble column.
(4) The spatial correlation between the cross-sectional averaged gas hold-up and local vorticity clearly indicates that the bubbles rising-up is strongly entrained by large spirally turbulent eddies with the trend of bubbles to cluster in the central region of the bubble column.
(5) Based on the eddy cell theory, the volumetric mass transfer coefficient estimated by using the modified SGS model can have better accuracy of estimation of the interfacial mass transfer between bubbles and liquid than that using the standard SGS model.

Acknowledgements
This work was financially supported by the National Natural