Euler/Euler large eddy simulation of bubbly flow in bubble columns under CO 2 chemisorption conditions

It has been generally recognised that the turbulent dispersion force plays an important role in interphase momentum transfer. However, the effect of the added mass stress on the momentum and mass transfer in bubble column bubbly flow has not been addressed appropriately. As the turbulent eddies in the surroundings of bubbles interact strongly with the rising bubbles in bubble column bubbly flow, such interactions will bring out the change of interfacial areas between the bubbles and carrier fluid, consequently leading to changes in the interfacial mass transfer. When employing large eddy simulation for modelling bubbly flow coupled with the chemisorption process


Introduction
Various carbon-neutralisation technologies have been intensively discussed in recent years.A typical example is CO 2 absorption or recycling in aqueous NaOH using bubble column.It has been widely recognised that the adoption of bubble column reactors is effective for CO 2 chemisorption due to their simple structures and high interfacial contact areas that can give high mass and heat transfer rates.However, the chemisorption process in the bubble column involves complex phenomena with the interfacial mass transfer, chemical reactions and hydrodynamics being strongly coupled and numerical simulation of such reactive gas-liquid two-phase flow in the bubble column remaining a challenge (Fig. 2).The key issue lies in a suitable modelling method for bubble volume fraction, volumetric mass transfer coefficient, interfacial area concentration, equivalent bubble diameter, enhancement factor and reaction rate constant.At present, two approaches to numerical computations of gas-liquid two-phase flows are employed: the two-fluid approach (Euler/Euler method) where both phases are treated as a continuum [43,57], and the Euler/Lagrange approach in which the fluid phase is treated as a continuum while the disperse phase involving a large number of bubbles is traced through the previously calculated flow field by taking into account relevant fluid forces.The two-fluid approach (Euler/Euler method) allows for efficient numerical computations of large and industrial scale processes such as bubble columns, but strong approximations and sophisticated closures are usually required for accurately describing the phase interaction.The adoption of Euler/Euler large eddy simulation for modelling of bubble column bubbly flows has received more attentions in recent years.Deen et al. [1] conducted the Euler/Euler LES simulation of bubbly flow in a rectangular bubble column and they have indicated that the use of LES for modelling of bubbly flow can well reproduce the results that are consistent with the experimental data in the high wave number (refer to small length scale) for transient bubble plume.In contrast, the use of RANS modelling approach for such bubbly flow has to employ an enhanced turbulent eddy viscosity and only the flow patterns corresponding to small wave numbers (large scale turbulent eddies) can be captured.By performing the Euler/Euler LES of bubbly flow, Lakehal et al. [33] pointed out that the Euler/Euler RANS modelling relies on the time averaging, which results in the local fluctuations caused by both local turbulence and dynamic interactions between phases being not properly described.However, these local fluctuations can be partly resolved to some extent in the Euler/Euler LES modelling as the spatial filtering is used instead of the time averaging.The use of LES directly resolves flow details with the scales greater than the filter width, typically equivalent to the imposed computing grid, while those flow details of the small eddies and the interactions between the bubbles and carrier phase at the sub-grid scale (SGS) are modelled [42].Ničeno et al. [45] employed the one-equation SGS turbulent kinetic energy Euler/Euler LES and hypothesized that the model's sub-grid scale kinetic energy can be used to estimate the SGS dispersion turbulent force.Ma et al. [2] presented detailed Euler-Euler Large Eddy Simulations (LES) of dispersed bubbly flow in a rectangular bubble column.The physical models describing the momentum exchange between the phases including drag, lift and wall force were chosen but based on previous experiences of the authors.Experimental data, Euler-Lagrange LES and unsteady Euler-Euler Reynolds-Averaged Navier-Stokes results were compared.They have indicated that the use of the proposed model combination provides good agreement with experimental data for the mean flow and liquid velocity fluctuations.The energy spectrum obtained from the resolved velocity of the Euler-Euler LES was also presented in their work.Ma et al. [3] further conducted Euler-Euler Large Eddy Simulations of dispersed bubbly flow in a rectangular bubble column at a low Reynolds number.The emphasis of the study was the analysis of bubbly flows concerning the influence of the bubble-induced turbulence model.Liu and Li [34] have presented the work on Euler-Euler large eddy simulation (EELES) of transient dispersed turbulent bubbly flows in a laboratory scale square crosssectioned bubble column.The sub-grid modelling based on the Smagorinsky kernel with dynamic C S constant was implemented.The bubble induced turbulence and various interfacial forces including drag, lift, and virtual mass forces are incorporated in the model.Good quantitative agreement with previous experimental measurements by a two-camera particle image velocimetry (PIV) was obtained both for the fluctuating velocities and the mean velocities.The classical − 5/3 law of power spectrum densities (PSD) of liquid vertical velocity was found to be held properly in the low frequency region but a scaling around − 25/3 was found for high frequency region.The classical − 5/3 law and − 3 scaling were found in both low and high frequency regions, respectively, in energy spectra, which was observed in the Euler/Euler LES study for bubbly flow in a cylindrical bubble column [36].Fard et al. [17] studied the performance of both one-equation SGS kinetic energy model and standard Smagorinsky SGS model in predicting the vertical-spiral and churn-turbulent flow in bubble column.Mohammadi et al. [4] presented an efficient mathematical model for conducting large eddy simulation (LES) of turbulent bubbly flows at low gas hold-ups.The mathematical model is based on the filtered multi-fluid equations derived by application of the component-weighted volume-averaging process of each phase.The filtered equations were formulated in generalized curvilinear coordinates to simulate turbulent gas-liquid bubbly flows in complex geometries.The Germano model [19] was used to account for sub grid turbulent stresses.The accuracy of the proposed method is assessed via comparison with the experimental data and other numerical results in the literature.Tabib and Schwarz [61] have examined the incorporation of SGS turbulent dispersion (TD) force in Euler/Euler LES work and determined that by employing a smaller magnitude of SGS-TD force, the results of a coarser mesh can be improved.However, in all of the aforementioned studies of Euler/Euler LES modelling, the effects on the filtering of the momentum interfacial exchange terms on the prediction of hydrodynamics, turbulence characteristics and mass transfer in bubble column bubbly flows have been rarely investigated.
With respect to the mass transfer in bubble column bubbly flows, both the Euler/Euler (E-E) and Euler/Lagrange (E-L) approaches have been used to describe the chemisorption process of carbon dioxide in bubble columns [1,11,12,21,23,26,27,29,30,31,39,62,63,70].For the Euler/Lagrange approach, Darmana et al. [12] used a discrete bubble model to simulate the CO 2 absorption in aqueous NaOH in the frame of Euler-Lagrange LES.Their simulation results indicated that the overall mass transfer rate was underestimated compared with the experimental data.The pioneering experimental and modelling work carried out by Darmana et al. was recovered by many other researchers [6,24,31,62,70].Taborda and Sommerfeld [62] studied the effect of bubble dynamic model initially proposed by Sommerfeld et al. [58] on the mass transfer during the CO 2 chemisorption, indicating that the mass transfer in CO 2 chemisorption can be better predicted by accurate estimation of the surface area of bubbles through considering bubble oscillation behaviour in the bubble column.After comparing with three different cases, they revealed that the case with consideration of the bubble oscillations and full dynamic model can deliver better results being consistent with the experimental data.While using the Euler/ Euler modelling approach, Chen [9] conducted a simulation on the carbon dioxide recycling process by considering four absorption systems in the bubble column.They indicated that the factors that have the impact on the mass transfer coefficient can be sequentially identified by the pH value, the concentration of CO 2 , the temperature and the introduced gas-flow rate.The mass transfer enhancement factor in their model had to be adjusted based on the experiments reported by Fleischer's group in order to achieve the better prediction results.It should be noted here that the enhancement factor plays an important role in evaluating the efficiency of the mass transfer for the reactive flows.The effect on the mass transfer using different enhancement factor models has been assessed in detail by Krauβ and Rzehak [31].The authors have proposed a simplified expression of the enhancement factor model which has shown to be more reasonable after comparing the simulation results.Furthermore, an engineering calculation applicable model which well expresses the chemisorption reaction between CO 2 and the water apart from OH -was implemented into their Euler/Euler modelling [30].However, the time-dependent change of the pH value was not well reflected in their work, which may be attributed to the fact that the coupling between the hydrodynamics and chemisorption reaction dynamics is still not well addressed.In addition, the mass transfer occurs in the region where the species concentration differences (Fig. 2), and is highly related to the estimation of the mass transfer coefficient and local interfacial area density.The bubble size variation will be directly influenced by the local interfacial area density.The effect of bubble size distribution at the inlet of the reactive bubble column on bubble dynamics was discussed by Hlawitschka et al. [24] who conducted Euler/Euler LES simulation of the reactive bubble column and compared the modelling results with their experimental data.They found that the change in the bubble size was not significantly affected by the reaction, and the predicted concentration distribution was overestimated for the region of top 1/3 of the bubble column.The adoption of the population balance model for accounting the bubble size changes taking the likely bubble coalescence and break-up was conducted by Buffo et al. [6] in their Euler/Euler two-fluid model simulations.They have reported that a generally similar trend and consistency with the experimental results on bubble volume fraction profiles but they had to employ an inlet bubble size distribution that has a standard deviation of 15% to the measured mean experimental value.Zhang et al. [70] conducted the Euler/Euler LES by incorporating the bubble number density model, revealing the coupling between the shear turbulence and mass transfer in the bubble column.
The correct prediction of the mass transfer in bubble columns is highly dependent upon both the turbulence modelling of shear turbulence and the description of the bubble-induced turbulence.The RANS (Reynolds averaged Navier-Stokes) modelling approaches have been widely used in conjunction with a turbulence model such as k-ε or k-ω.The turbulence models for describing the bubble-induced turbulence (BIT) have also been analysed by Zhang et al. [69], showing that no major differences existed when including the BIT.However, the impact on the mass transfer by the BIT was not assessed in these previous studies.It should be mentioned that a more recent work of Magolan et al. [5], considered the BIT-correlations but has indicated the existence of remarkable differences in the computed turbulent kinetic energy for the formulations considered.As large eddy simulations (LES) are increasingly being used for the calculation of bubbly flows, it has been accepted that the fine turbulent structures can be described by comparatively simple models, e.g. the use of Smagorinsky model (Deen et al., 2001); [11,20,25,48,60], and the effects on the transfer process can be better evaluated.For Euler/Euler two-fluid LES modelling, the turbulence vortex structures are resolved up to the filter width, corresponding approximately to the grid size, and thus the bubble transport by these vortices is resolved.Most of the cases employing such LES modelling have neglected the influence of the fine structure turbulence (SGS: subgrid scale) on the bubble motion and the effect of the entrained bubbles on the SGS [48] since their contribution to the energy of the fluid fluctuation is comparatively low.However, as indicated in the review article by Dhotre et al. [14], the use of LES for disperse bubbly flows can correctly reproduce the turbulent kinetic energy only when modelling bubble-induced turbulence (BIT) is conducted.This has clearly indicated consideration of the BIT has a decisive influence on the mass transfer and a subsequent chemical reaction in the bubble column.
From the available research, the effects of consideration of turbulent dispersion and added mass stress on the interfacial mass transfer in the Euler/Euler LES modelling of bubble column bubbly flow are rarely discussed, i.e., the term ∇•(α k τ k ).The appearance of the additional added mass stress is the consequence when filtering one of the interphase forces, the added mass force.The present work will particularly assess the impact of including the added mass stress in the LES modelling on the mass transfer of the CO 2 chemisorption in the bubble column.In addition, since the range of the bubble size in the chemisorption process in the present study falls into the integral turbulent length scale (Taylor length scale ~4 mm) and inertia sub-range, the effects of such interaction between the bubbles and the surrounding eddies on the eddy effective diffusivity will be taken into account.This paper will be organised as follows.Section 2 presents the mathematical modelling used in present work and related formula derived for those modified terms that will be particularly addressed.The numerical set-ups for the validation of the Euler/Euler LES modelling of CO 2 chemisorption reactive bubbly flow are described in Section 3 while Section 4 presents the predicted results obtained from Euler/Euler LES modelling by employing the three different mass transfer models together with the comparisons with the available experimental work plus the detailed discussion.The conclusions reached from the present study are given in Section 5. Two appendices are also provided with the correlations for reaction kinetics and physico-chemical properties that were used in the study and the empirical expression based on the fitting for the enhancement factor of an instantaneous irreversible second order reaction.

Governing equations
The two-fluid model based on the Euler/Euler LES turbulence modelling was adopted in this work.It was assumed that each fluid (or phase) is treated as a continuum in any size of domain under consideration while both phases share the domain and can interpenetrate as they are transported within it.The Eulerian modelling frame work employs the phasic function of presence to account for the space occupation by either the liquid or bubbles.This allows the mass and momentum transport equations for each phase to be directly derived by applying LES filtering.In this LES model, each parameter φ is characterised as the combination of the Favre phase-weight filtered part φ that needs to be resolved during the filtering process and the unresolved part φ ' that needs to be modelled by using the SGS model for closure.For convenience, the double overbar symbol "=" will be dropped for convenience hereafter.However, when analysing the added mass stress that will be discussed in later of this section, the double overbars will be adopted again for clarification of the derivation process.The phase-weight filtered conservation equations for continuity and momentum can then be written as.

∂ ∂t
where index k = G, L stands for the gas (bubble) and liquid, respectively, u k refers to the local instantaneous velocity vector and α k stands for the volume fraction satisfying α G + α L = 1.On the left hand side of Equation ( 1), the terms represent the change rate of the phasic mass and transport due to advection.The term on the right hand side indicates the mass transfer occurred between phases, which is given by.
where the carbon dioxide chemisorption process is concerned.The terms on the right-hand side of Equation ( 2) represent the stresses due to viscous and turbulent shear, the pressure gradient, gravity and the filtered momentum exchange between the phases, caused by the actions of the interface forces, respectively.The stress term of phase k can be expressed as: where μ eff L is the effective viscosity for the carrier phase, which is assumed to be composed of the contributions from the molecular viscosity μ L , the shear induced turbulence viscosity μ T and the extra viscosity due to bubble-induced turbulence μ BI , given by.
The viscosity due to the bubble induced turbulence can be modelled following the work of Sato and Sekoguchi [55], which is given by Equation (6), It should be pointed out that the turbulent viscosity model accounting for the BIT in two-phase flow as proposed by Sato and Sekoguchi [55] was derived by assuming the flow about a fixed bubble as the flow about a cylinder.In reality, the bubbles will interact with the turbulent eddies and response to the turbulent eddies that have entrained the bubbles, which will generate the slip fluctuations.This response should also 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.Therefore, the instantaneous fluctuation of bubbles would always differ from it of the surrounding turbulent eddies in subgrid scales, especially for the eddies having the similar size with the bubble diameter [32].Based on our previous work [36] that has followed Garcia's work [18] but also considering the contribution from the added mass on bubble translation, the relative velocity between the two S. Long et al. phases can be evaluated by, u λ (ελ) where λ represents the different turbulent length scales in the range between the integral and Kolmogorov scales (L > λ > η).When the derivative of u λ equals zero at a certain λ, u λ will have a maximum, as defined by.
Substituting λ * into Equation ( 7) yields. ( In this turbulent length scale range, (dε) 1/3 can be regarded as the fluctuating velocity of the bubble.Thus, u * λ 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.Considering the bubble response to the eddies and the interaction between bubbles and eddies with a Stokes number, defined by St SGS = τbubble τL,SGS , the Smagorinksy model of sub-grid eddy viscosity can be modified based on the correlation between the fluctuating velocities of the bubbles and liquid in terms of the turbulent eddies with the length scales falling into the inertia subrange [32], expressed by Equation (11).
As demonstrated in Equation ( 11), the instantaneous fluctuation of bubbles would always smaller from the surrounding turbulent eddies' fluctuation in sub-grid scales, especially for the eddies having the similar size with the bubble diameter.When Equation ( 8) is implemented into the sub-grid scale, the relationship can be replaced by where the bubble response time scale can be estimated by as proposed by Sommerfeld et al. [58].The SGS turbulent eddy turn over time can be estimated by τ L,SGS = Δ/u ′ L,SGS .Using u ′ 2 L (λε) 2/3 and substituting Equation ( 11) into (10), the modified slip velocity can be expressed by Equation (12): The turbulence dissipation due to the bubbles corresponds to the inertial subrange can be assumed that mainly occurs when the eddy integral scale λ approximately equals Δ and the dissipation can be represented by Equation (13): Thus introducing Equation ( 12) into (13) yields.
where the value of C b can be determined using the trail simulation.In the present study, C b with the value of 0.9, 0.7, and 0.5 has been tested.It was found that the use of a value of 0.7 may give rise to the predictions that have better consistency with the experimental data [12].We thus have adopted C b = 0.7 throughout all the Euler/Euler LES modelling in the present work.The total turbulence dissipation due to liquid phase shear induced turbulence and the turbulence due to bubble response to SGS eddies can be estimated by.
By employing the eddy viscosity model, the liquid-phase turbulence coupled with the modified Smagorinsky SGS viscosity [56] due to bubble response to the eddies can be written as.
where Cs is the model constant.The SGS filtering length scale Δ takes the value based on the grid control volume as Δ = (V cell ) 1 3 .The modified SGS eddy viscosity model considering bubble dynamic responses to eddies has been implemented into our Euler/Euler LES modelling.
For estimating the effective viscosity of gas phase, Equation ( 17) is adopted.
In the present LES model, large eddies with the length scale being greater than the filter size are directly resolved while those small eddies with the size smaller than the filter size are modelled using the SGS model.It thus can be imagined that the eddies smaller than Δ = (V cell ) will give rise to the local fluctuations, implying the potential impact on the interfacial mass transfer when the bubbles are entrained in the bubble column.If no coalescence and breakup is assumed with the condition of low bubble volume fraction, α G < 5%, the shape of bubbles are more likely to be spherical or ellipsoidal with little occurrence of breakup and coalescence [10,64].In such case, the bubble number density can be defined by Equation (18): Therefore, the interfacial mass transfer occurred between phases would subsequently give rise to the bubble size change as the bubble equivalent size variation can be associate with the zero-th moment of bubble size distribution (BSD) through

Interfacial momentum exchange modelling
In Euler/Euler LES modelling of bubbly flows, the total interfacial force arising as the action of the interfacial momentum exchange between the two phases is usually characterised by the contributions from several independent physical effects so that the interfacial force M F can be written as.
The forces indicated in Equation ( 19) represent the contributions from the interphase drag force, lift force, added mass force and turbulent dispersion force.These forces have been obtained by using the phase weighted filtering.The drag force is due to the resistance experienced by the bubble transport in the liquid.Both viscous and turbulent stresses generate the skin drag and pressure distribution around the moving bubbles gives rise to the form drag, especially when the boundary layer on the bubble surface separates to generate the bubble wakes.The lift force arises from the net effect of pressure and stress acting on the surfaces of the bubbles, which are strongly associated with the averaged shear gradient acting on the bubbles and the orientation of slip velocity.The turbulent dispersion force can be regarded as the turbulent diffusion of the dispersed phase by those turbulent eddies that are strongly interact with the bubbles.It should be noted here that the turbulent dispersion force is the result of phase-weight filtering of the drag force acting on the bubbles.Thus, for LES modelling of the gas-liquid system, the contribution from the added mass force needs to be evaluated by considering the filtering process itself.Most of recent work on Euler/ Euler LES modelling of bubbly flows have only included the filtered added mass force [8] [13,17,69].However, in the case of large eddy simulation, there would be an additional force at the SGS level in addition to the turbulent dispersion force, i.e., added mass stress (SGS-AMS), when phase-weight filtering the added mass force term.An attempt is made in the present study to take account of sub-grid-scale added mass stress (see Fig. 2).The interphase momentum transfer due to drag force after taking filtering is given by.
The drag coefficient for closure is the correlation proposed by Tomiyama [65], where the Eötvos number σ is the ratio of the bubble buoyancy force to the surface tension.The lift force acting on bubbles can be estimated by.
where Re B is the bubble Reynolds number measured based on the equivalent bubble diameter in case of deformed bubbles being concerned and The filtered added mass force can be estimated using Equation (23), where the added mass coefficient C AM takes the typical value of 0.5 which is suitable for a spherical bubble.
The common used turbulent dispersion model is proposed by Burns et al. [8] employing Farve averaged drag force model in RANS modelling.After employing the spatial filtering to the drag force in this Euler/ Euler LES work, the turbulent dispersion force can be obtained using the eddy viscosity hypothesis.With the consideration of the bubble-eddy dynamic interactions in the eddy viscosity model in Equation ( 16), the turbulent dispersion force can be expressed as.
where the SGS turbulent eddy viscosity used is presented in Appendix A (10).As it is postulated that the sub-grid fluctuations correspond to those SGS turbulent eddies, the SGS turbulent eddies would significantly affect the bubble transport and dispersion when they act on the surface of the bubbles, leading to the bubble continuous deformation or oscillation when observing to follow their trajectories accordingly (Fig. 1).As the added mass force acting on the rising bubbles corresponds to the liquid around the front part of the bubble, it can be imaged that such SGS turbulent eddies will have an important effect on the interfacial momentum exchange and mass transfer.We are now employing the phaseweighted filtering to the added mass force in Euler/Euler LES modelling.
The instantaneous added mass force can be expressed as.
By introducing the characteristic function of the phase presence χ K (K = G for bubbles, K = L for liquid) and taking the filtering, one has.
where the 'overbar' has been used to express the filtering.With the phase-weight filtered and fluctuating components of u K and u ′ K at given point being respectively expressed by.
the phase-weight filtering of the added mass force can be expressed as the contributions from the filtered added mass force and the added mass stress force, i.e.
) ) The last term in the right-hand side of equation ( 28) can be neglected when comparing its contribution with those of the first two terms as the term χ G u ′ L can be modelled by assuming to be proportional to the bubble volume fraction gradient [67] and this assumption can be also applicable to LES filtering, i.e.
where ν T,SGS,L is the liquid turbulent eddy viscosity and σ is the bubble diffusivity.Substitution of equation ( 29) into the last term of the righthand side of Equation ( 28) yields.
It can be seen from the above expression that the terms in Equation ( 3) are small in the magnitudes in the fully developed bubble column bubbly flow.Thus Equation ( 28) can be simplified as, ) ) where the double-overbar '=' is dropped again for consistency.It should be noted that the consequence of applying spatial filtering to the added mass force would deal with the correlation u ' ki •∇u ' kj as indicated in the second part of the right-side of Equation (31), which functions like the Reynolds stress but also correlates with the local bubble volume fraction fluctuation, thus being referred to as the SGS added mass stress (SGS-AMS).By employing the eddy diffusivity hypothesis, the SGS added mass stress (SGS-AMS) can be formulated, which can be given by.
where τ L and τ G are defined by u ' L,i u '  L,j and u ' G,i u ' G,j , respectively and are modelled based on Equation (15) during the simulation.As the bubbles are subjected to the eddy fluctuations, it can be expected that the SGS-AMS force will have a significant impact on the interfacial mass transfer between the bubbles and liquid phase and this effect will be discussed in Section 4.

Chemisorption process description and the involved interfacial mass transfer
Generally two reactions can be assumed to take place in chemisorption of carbon dioxide bubbles in the NaOH solution using the bubble column reactor.Firstly, there is a physical absorption of CO 2 from gas phase to liquid, CO 2 (g)→CO 2 (aq) (33) in which two reversible reactions simultaneously take place, described by.
The reaction rates can be evaluated by.
where k i± represents the reaction rate constant for the above reversible reactions, i = 1, 2 denote the first and second reaction, respectively, while "+" and "-" stand for the forward and backward reaction separately.The detailed estimation process of the reaction rate constant can be found in the Appendix.The interfacial mass transfer occurring in the above mentioned chemisorption process can be described by using the mass fraction Y j of each species j in the liquid mixture and the species transport equation for the liquid phase is given by.

∂ ∂t
where D j L,eff is the liquid phase effective diffusivity.Considering that the original expression for the turbulent mass effective diffusivity is usually expressed by the sum of molecular diffusivity and eddy mass diffusivity, i.e.
ρSct , it is crucial to estimate the turbulent eddy viscosity and Schmidt number.In order to mimic the real condition, the eddy viscosity model used in this study has implemented the modified SGS Smagorinsky model into the simulation, which has taken the response of bubbles to the SGS eddies into account.In other words, the eddy viscosity used in Euler/Euler LES should include the contributions from the SGS shear induced viscosity and from the interactions between bubbles and SGS eddies.As the eddy turbulent fluctuations will have a significant impact on the interfacial mass transfer for bubble size which have an equivalent size of the SGS grid and falls into the inertia subrange, the mass effective diffusivity that takes the eddy turbulent diffusivity into account can be expressed as,.
S. Long et al.where Sc t is the turbulent Schmidt number.The related parameters are listed in Table A1.The relationship between the molecular diffusivity of carbon dioxide in water adopts the one proposed by Ratcliff and Holdcroft [49], given by.
The mass transfer rate of species CO2 from bubbles to the liquid can be defined as.
where k L is the mass transfer coefficient for liquid side CO 2 , dB .The interfacial area concentration can be estimated by a = 6α G /d B if the bubbles are assumed to be spherical.The correlation which relates the mass fraction of CO 2 (aq) for both sides of the bubble surface can be written as Equation ( 44) by using the henry constant and the equivalent mass fraction of aqueous carbon dioxide in the liquid side can be thus defined as.
where the solubility of CO 2 in water, characterised by Henry constant H CO2 applies to the condition that the solute concentration is low.It should be pointed out that the use of the ratio of the CO 2 concentration in the liquid to the CO 2 concentration in the gas at equilibrium is appropriate from perspective point of view of the numerical simulation.
For the estimation of Henry constant H CO2 , Versteeg and van Swaaij (1988) have proposed a correlation for the temperature dependency of pure water based on their own experimental data, given by.
The interfacial mass transfer occurred can be characterised by the overall mass transfer coefficient, k L , which can be obtained from the following Sherwood number correlations for bubbles that are dependent on Reynolds and Schmidt numbers [2,35]; Brauer and Soudack [4].For non-spherical bubbles, the correlation is derived by accounting for the stochastic deformations of the interface caused by turbulent eddy fluctuations acting on the bubble surface.The overall mass transfer coefficient estimation in the present work is based on Equation (46): The interfacial mass transfer is also influenced by the pH value of the solute.Such an effect is usually accounted for by introducing the enhancement factor.The enhancement factor dependency on pH value was investigated by Fleisher et al. [9] who have replicated the same experiment on CO 2 chemisorption in NaOH solution as done by Darmana et al. [12].They revealed that the chemical reactions in CO 2 chemisorption has an influence on the absorption process.They proposed the enhancement factor which can be estimated by the following expression given by Equation ( 47): Thus, using Equation ( 40) coupled with Equations ( 41)-( 47), the mass fraction distributions of each species,   Y CO2l can be obtained using the modified LES SGS model.

Numerical simulation
The modified LES models proposed in the present study were validated by comparing with the experimental data carried out by Darmana et al. [11,12] and our collaborate research team, herein referred to as Case 1, Case 2 and Case 3. The schematics of the configurations of three bubble columns studied in this work are illustrated in Fig. 3.In the LES modelling, the reference position of the coordinates was chosen to be the centre of the cross-section at the bottom of the bubble columns.The columns in three cases were assumed to be filled with a NaOH solution with the given liquid level.The gas aerator in case 1 was located at the centre of the bottom plane of the reactor via 21 needles within a 5 mm square pitch.Air or carbon dioxide was injected to the system via a perforated plate with 40 holes evenly distributed at a 6.25 mm square plate for Case 2. While for Case 3, the gas was introduced into the column via a capillary gas sparger equipped with 13 needles with an inner diameter of 0.57 mm.The detailed information is listed in Table 1.
The initial bubble sizes adopted in the LES modelling were 5.5 mm (Case 1, 3) and 4 mm (Case 2), combined with the use of bubble number density equation to compute the bubble size, which is essential factor in calculating the related coefficients in evaluating the interfacial mass transfer.The simulation domains for the three cases were partitioned with the computational grids of d B /Δ ≈ 0.733 for Case 1, 0.72 7 for Case 2 and 0.675 for Case 3 in the central region of the bubble columns with a growth rate of 1.2 from the walls, which satisfy the constraints suggested by Milelli's criterion for bubbly flows using Euler/Euler two-fluid LES modelling [42].The mesh set-ups with 48,000, 49,500 and 52,920 mesh cells were adopted for three cases, respectively.In order to clearly illustrate the impact from the chemisorption process on the hydrodynamics in the bubble column bubbly flows, nitrogen was assumed to be injected into the column at the beginning and CO 2 was then supplied.Similarly, the LES simulation also adopted the procedure with N 2 being supplied without chemical reaction for the first 20 s of the simulation and the solutions of the species transport equations being turned on afterwards.The simulation time for each case was set according to the published work on the numerical simulation work and the experimental work.For example, for Case 1 [12], the corresponding numerical work by Darmana et al. [12] has adopted 250 s; while for Case 2 [11], the numerical validation work conducted by Zhang et al. [70] used 250 s.For Case 3 carried out by the collaborator [63] is also experimentally repeated in this work (details can be found in Appendix-B), the corresponding simulation work done by Taborda et al. [63] was last for 250 s.The simulation results obtained by the LES during the given time interval were collected for time averaging.To better capture the transient dynamics of the turbulent eddy development in the bubble columns, the time step adopted in the simulation was chosen in terms of the CFL criteria, i.e., min( |uL|δtE Δ , |uG− uL|δtE Δ ) < 1.0, but was gradually increased from 0.0005 to 0.001 s.In the LES simulations, the contact between the liquid phase and the bubble column wall was set as a no-slip condition, whereas a free-slip condition was used for the dispersed phase with the assumption of minimal direct contact between the bubbles and the walls.The turbulent wall function with assuming a smooth wall of the bubble column was utilized to eliminate the requirement to resolve the viscous sublayer for a very small y + .Constant relative static pressure (P rel = 0) was used for the outlet at the top of the bubble columns.At the inlet, the normal gas superficial velocity and mass fraction were specified according to the experimental gas superficial velocities and the gas fluxes.The mass fractions were computed by analysing five species transport equations for CO 2 (aq), CO

Results and discussion
In order to validate the reliability of the proposed modified SGS-AMS model coupled with the consideration the effect of the modification by bubbles dynamic response to the surrounding eddies on the bubble induced turbulence eddy diffusivity, the LES simulations were conducted for both rectangular (Cases 1 and 2) and cylindrical (Case 3) bubble columns based on the experimental data as reported by Darmana et al. [11,12] and the experiment conducted by Taborda et al. [63] and our experimental data which repeated Taborda et al.'s experiment by using the same experimental set-up and the experimental conditions, detailed in Appendix-B).Corresponding to each experimental case, the simulations were conducted by including (a) conventional turbulent dispersion and mean added mass forces (hereafter referred to as Model A); (b) the modified SGS turbulent dispersion (SGS-TDF) and mean added mass forces with coefficient C AM = 0.5 (hereafter referred to as Model B) and (c) the modified SGS turbulent dispersion (SGS-TDF) and the SGS added mass stress force (SGS-AMS) (hereafter referred to as Model C).In the following subsections, the overall mass transfer in the CO 2 chemisorption process will be comprehensively discussed, focusing on the influences of accounting for the added mass stress force term arisen as the consequence of spatial filtering on the interfacial mass transfer, time-averaged velocities, bubble volume fraction profiles, liquid phase turbulent kinetic energy spectrum and the species concentration spectrum.

Effect of SGS-AMS on CO 2 chemisorption process
The time history of the species concentration and pH variation obtained by including the SGS-AMS in the LES modelling are compared with the two experimental cases reported from Darmana et al. [12].The time-dependent predicted pH and species concentration profiles acquired at middle point of z = 980 mm using three different models A, B and C are shown in Fig. 4. In the figure, green dot lines represent the predicted concentration by using the standard force model (Model A), blue dash-dot lines stand for the profile using the standard added mass force model and modified SGS-TDF (Model B), while the red dashed lines are the one with modified SGS-TDF and SGS-AMS (Model C).Generally, the predicted species concentration evolution and the pH value variation obtained in three LES simulations with and without the modifications are in good agreement with the experimental data.During the bubble rising-up, CO 2 may be physically dissolved into the liquid across the interface, leading to an increase in CO 2 (aq) concentration.As a result, this process triggers the chemical reactions in the liquid phase.It can be   seen from Fig. 4 that for the first 75 s of simulation time, carbonate is accumulated during the fast consumption of hydroxyl ions, where the dissolved CO 2 are fully interacting.This phenomenon is well reflected by the apparent change in slope of pH curve.The products are carried by the large circulation among the bubble column to the top region as well as the downwards recirculation near the wall region.When the initial OH − is totally consumed, the carbonate concentration reaches the maximum and begins to drop, the bicarbonate concentration starts to increase until the initial concentration of hydroxide ions reached.At the same time, the aqueous carbon dioxide starts to store since all of the hydroxide has been used.By comparison, the simulation results predicted by the model without including the SGS-AMS modification have an obvious delay in the species concentration variation when the same reaction rate constant, mass transfer coefficient were used.A faster reaction was found when using the SGS-TDF and SGS-AMS models (Model C) while the use of the SGS-TDF and SGS-AMS without modifying the eddy viscosity (Model B) predicts a reaction with the rate which lies in the predicted reaction rates by using Model A and Model C.This phenomenon indicates that the influence of the modification is well reflected in the mass diffusivity estimation when solving the transport equation.By considering the SGS-TDF and SGS-AMS, the dispersion of the bubbles in transverse or radial direction may be better estimated and bubble cluster wobbling which can lead to a longer residence time.Also, the inclusion of the SGS-AMS term in the LES modelling seemed to yield a better estimation for momentum exchange due to those turbulent eddies in the zones characterised by added mass surrounding the bubbles, giving rise to a higher interfacial mass transfer rate, consequently leading to the predictions for an earlier dissolution and stimulating the following reactions.Thus, the predicted concentration profiles by employing the Model C clearly have a closer trend to the experimental data.This reveals an important fact that the contribution of the filtered SGS-AMS to the interfacial mass transfer when using the LES modelling approach cannot be neglected.
For Case 2, the evolution of domain-averaged pH and CO 2− 3 concentration predicted by employing three models A, B and C were also compared with the LES simulation results reported by Zhang et al. [70] as shown in Fig. 5.The productions of CO 2− 3 predicted by our three approaches are similar to their results but a faster reaction progress can be observed when employing both filtered TDF and AMS models.It was also found that when neglecting the contributions from the SGS-TDF and SGS-AMS (Model A), the predicted reaction progress has a significant delay and reaches the equilibrium much later.
Apart from the results for Cases 1 and 2 which are based on the bubble column with the rectangular cross-section, similar finding was also observed in the case of cylindrical bubble column (Case 3).The pHvalue time history was predicted by using the three models as well, which is shown in Fig. 6.It can be seen from the figure that by using the modified models (Case 3-Model B and Model C), the predictions are noticeably improved with a better performances by using the modified models comparing to the results obtained without considering the contribution of the SGS-TDF and SGS-AMS.It should be noted that the slope of the pH time history curve changes with the time evolution, which is caused by the shift in equilibrium of the first reaction mechanism in terms of the bicarbonate expressed in Equations ( 34) and (35).It seems that applying Model B and Model C, an improvement in predicting the reaction progress is strongly evidenced with better prediction being achieved using Model C.

Effect of inclusion of SGS-TDF and SGS-AMS on hydrodynamics
The better performance of adopting the SGS-TDF and SGS-AMS models (Model C) in the predicted species concentration evolution may be attributed to a better description of the bubble dispersion and liquid shear stress so that the bubble residence time in the bubble column can be better predicted.To assess this speculation, the timeaveraged bubble axial velocity at z = 750 mm predicted by using the three models are compared with the experimental results, as shown in Fig. 7.The time-averaged bubble axial velocity is calculated based on the following relationship: For the non-reactive case, one can observe that there is not significant differences in the magnitude of the maximum averaged velocity using the three models but the large difference can be observed in the positions away from the central in the radial distribution for Model A. Nevertheless, for the reactive case with CO 2 injected, it can be seen from the figure that a good agreement with the experimental data is attained for the Model B and Model C. For the simulation conducted without considering the SGS-TDF and SGS-AMS, the bubble axial velocity seems to be much over-predicted for Model A, implying that the lateral bubble distribution being not well described due to the underestimated lateral bubble dispersion.Consequently, the higher bubble volume fraction in the core region will induce higher bubble rising velocity so that an overestimation of the bubble axial velocity takes place.However, when considering the effect of SGS-TDF and SGS-AMS, a decrease in the bubble axial velocity is found, giving rise to an adequate estimation of the bubble axial velocity distribution in the transverse direction.The comparisons using the three models clearly indicates that the SGS-AMS plays an important role in the bubble dispersion and bubble dynamics.
Figure 8(a) depicts the predicted time-averaged laterally bubble axial velocity distribution at H/D = 7 for the non-reactive case with injection of N 2 , comparing with the experimental data in this work.To better consider the effect of the cross-section of the bubble columns (Model A, B and C) on the hydrodynamics, the bubble lateral velocity profiles are obtained by the following average method, ) dxdy where T 0 and T N are the start and end times of sampling from the LES modelling simulation, which take T 0 -T N = 250 s in the time averaging.It can be found from the figure that the higher bubble axial velocity in the center region is predicted (Case 3-Model A), which has a maximum value of 0.42 m/s at the axis.The maximum averaged bubble axial velocity predicted for Case 3-Model B and Case 3-Model C has respective 4.707% and 7.693% reductions compared with the case of using Model A, which is more consistent with the experimental data.However, the prediction performance becomes poorer in the region 0.95 < r/R < 1.0, which may be attributed to the conflict between of mesh refinement in the LES modelling and Milelli's limitation [42], i.e., 1.2< Δ /d B < 1.5, for two-fluid bubbly flow simulation.Although the relative smaller bubbles are more likely to be pushed towards the wall due to actions from the lateral dispersion and lift force, a finer mesh near wall region required by LES (5 < Δr + 〈3 0) still lead to a<1 mm grid imposed at the wall boundary, which will affect the accuracy of the LES modelling results.This remains to be resolved in further studies using the two-fluid Euler/Euler LES modelling for bubbly flow problems.As a key indicator in evaluating the correct estimation of the mass transfer in chemisorption of CO 2 , the LES predicted bubble volume fraction distribution profiles are compared with the experiment data.As shown in Fig. 8(b), the numerical results obtained by using both model B and C show good agreements with the experimental data.It can be seen from the figure that a small reduction in the central region and a consistent trend in the slope of bubble volume fraction profiles at 0.7 < r/R < 0.9 with the experimental data can be found, this may indicate the local eddy fluctuation induced shear has been well captured by the use of SGS-AMS model, resulting in a larger bubble surface area.In the present study, .To further investigate the potential reasons of the variations, the focus may shift to the predicted local bubble size distribution profile.
The local bubble size distribution has a significant impact on the hydrodynamics and overall interfacial mass transfer in the bubble column bubbly flows.To assess the effect of LES modelling implementing the SGS-TDF and SGS-AMS terms on the interfacial mass transfer in the bubble column, the cross-sectional averaged and time-averaged bubble mean Sauter diameter variation along the height for Case 2 is shown in Fig. 9.It can be seen from the figure that similar trends for bubble mean Sauter diameter variation are observed for all the models A, B and C, i.e., the bubble diameter decreases with the increase in the axial height.However, comparisons of all the three models with the experimental data reported by Darmana et al. [6] show that the predicted bubble diameter is not well matched with the experimentally observed bubble diameter in the beginning part of the bubble column where bubbles are injected into the bubble column.Darmana et al. [6] have observed that the bubbles within the injection region z < 0.4 m frequently cluster in the central zone, which may result in the difficulty in detection of accurate bubble diameter in the experiments.As the bubbles rise up, the influence from the large induced eddies by the four corners of the rectangular reactor becomes weaker and the bubble size measurement may become more reliable, especially when z > 0.4 m.Without any modification to the conventional SGS-TDF and SGS-AMS models, a small reduction in the predicted bubble diameter is found as comparing with the simulation result reported by Darmana et al. [6] After taking the SGS-TDF bubble dispersion into account, it can be seen from the figure that a steady decline in the predicted equivalent diameter was achieved in our LES simulation.Consistent with the bubble axial velocity prediction results, the implementation of the SGS-TDF gives a better bubble lateral dispersion estimation, resulting in a higher bubble volume fraction gradient in core region, which corresponds to the higher species concentration gradient and thus higher interfacial mass transfer.It should be noted that a consistency in the LES predicted bubble diameter with the experimental data using model C is achieved especially for the region in higher part of the bubble column.This may be attributed to the , where the association between the bubble diameter and the local bubble volume fraction can be identified.Thus, the higher the axial height from the gas sparger, the larger of a difference between the simulation results and the experimental observation without the modification on eddy viscosity and the turbulent dispersion.Furthermore, when using Model C, the correlation of turbulent eddy fluctuation induced shear and local bubble volume fraction for the bubbles ∇•(α k τk) α k ρ k would give rise to the local bubble volume fraction change, implying that the bubble diameter may be affected and this would be reflected in the estimation of the bubble Reynolds number, Re Based on the correlation proposed by Brauer [5], the overall interfacial mass transfer rate will be enhanced.It should be noted that the difference in the estimation of the bubble size between the use of model A and model C can reach about 11% of the original bubble size.Although the bubble size is still over-predicted by using model C, the comparison among the three numerically generated curves of d B /d initial reveals that the inclusion of SGS-TDF and SGS-AMS in the LES delivers the result that is consistent with the experimentally measured bubble size in the trend of bubble size change along the bubble column height except for the lower part of the bubble column.The reason for the overestimated profile indicates that the modified model still under-predicts the mass transfer rate in the Euler/Euler LES of chemisorption of CO2 in the bubble column.As can be seen from the equation: ṁj , the mass transfer rate is strongly related to the bubble interfacial area and the local species concentration while the estimation of Y j* L is also related to local bubble-eddy interaction, i.e. local turbulence generation.As pointed out by Yeoh and Tu [68], and Zhang et al. [70], the transport equation used to describe the interfacial area concentration change may be simply replaced by solving the bubble number density transport equation if no bubble breakage or coalescence takes place.By doing so, it becomes easier to track the bubble size using the bubble number density since the bubble number does not change with the given interfacial mass transfer.This allows to treat the bubble number density as a passive scalar.In the present study, the local .As the cases for study are in low bubble volume fraction, α G < 5%, the shape of bubbles is more likely to be spherical or ellipsoidal with little occurrence of breakup and coalescence [10,64], no coalescence and breakup is assumed in the present Euler/Euler LES.Thus, under-prediction of Y j* L would require a larger interfacial surface area tomaintain the mass transfer rate, which corresponds to an increase in the prediction of the bubble size.The negligence of bubble coalescence and breakup may affect the estimation of interfacial area concentration, but the implementation of bubble breakup and coalescence model into Euler/Euler LES remains challengeable.

Effect of inclusion of SGS-TDF and SGS-AMS on bubble dynamics
Figure 10 shows the instantaneous velocity vector highlighted by the local air volume fraction at different time steps with the results predicted using Model A, B and C from left to right.The transient meandering behaviour of the rising bubble group wobbling is captured using the proposed models.It can be seen from the figure that relative higher bubble volume fraction takes place in the central region of the bubble column, indicating that the bubble transport is significantly affected by the large turbulent eddies and the bubbles are entrapped by these large eddies to form the meandering phenomena.It can be seen from the snapshots that the colours is weaken in the region close to the bubble column wall, which implies fewer bubbles being existed or lower bubble volume fraction, when using Model A. However, the spread of bubbles is apparently evidenced when employing Model C in the modelling.Even when using Model B without considering the SGS-AMS, the lateral bubble dispersion can be observed as a result of considering the SGS-TDF.
To highlight the species concentration evolution in the chemisorption reaction process, the contours of CO 2 concentration and pH value in the X-Z cutting plane through the bubble column are displayed in Fig. 11 by using Model C for several time instants.It can be seen from the figure that at the beginning of the chemisorption process, only traceable amounts of dissolved CO 2 were detected over the whole region of the liquid solution in the bubble column, which gives out a slow reaction.At this point, the pH value of the liquid is obviously to indicate alkali-like, and very minor variations were noticed over time intervals <50 s (see also Fig. 6).As a result of the reactions, the distribution of the pH value appears to be more homogeneous than the CO 2 distribution does.As the mass transfer progresses, the hydroxide is quickly consumed and the pHvalue is further reduced (see Fig. 6).Due to the restraint of the bubble column cross-sectional shape, bubbly flow at z < 0.4H in the bubble column exists for the large recirculation vortices around the bottom corner of the bubble column, which results in poor interfacial mass transfer while a better interfacial mass transfer occurs in the upper region of the bubble column.With the chemisorption reaction process continuing, the bubbly flow patterns in the bubble column become dynamically stable for the period of 50 to 100 s.At 150 s, the dissolved CO 2 is found to be well distributed in the liquid NaOH solution and the predicted pH value throughout the bubble column distributes uniformly.
The instantaneous bubble locations observed from Darmana's et al. work [12] are shown in Fig. 12 as a qualitative contrast to the contours by the present LES simulation for Case 2 using Model C. In Fig. 12, the illustrations from the left to right are the instantaneous bubble locations, the predicted bubble volume fraction, the iso-surfaces of bubble volume fraction of 0.03, the bubble velocity, the liquid velocity and the pHvalues, respectively.The bubble volume fraction contours clearly show the bubble group being wobbling and exhibiting the trace of "S" shape, which matches the experimental observation.The gradient of gas hold-up near the wall, on the other hand, is higher in the simulation than in the experiment.As the chemisorption process is going on, the bubble volume fraction profiles steadily reduce the magnitude with increase of the height.It can be also seen from Fig. 12 that the pH variance is coupled with the dissolved CO 2 distribution and large turbulent eddy evolution.

Turbulent liquid kinetic energy spectrum and chemical species concentration spectrum
The turbulent kinetic energy power energy spectrum based on the axial liquid velocity fluctuations obtained at the centreline of the bubble column at Z = 0.75H are shown in Fig. 13(a).The predicted axial liquid velocity fluctuations were converted to a format of two-point correlation based on time development.The power spectrum density per frequency was then obtained by taking the Fast Fourier Transform of the timecorrelation.To have a deeper analysis of the power energy spectrum, the bubble representative frequency was also estimated for Case 2-Model C as reported by Prakash et al. [47], which is given by.Thus, it can be thought that the bubble induced turbulence energy was fed into in the system with the above induced eddy frequency.Following Risso and Ellinggsen [52], the typical wave number can be estimated by.
As revealed by the accessible experimental studies or DNS simulation data [37,47,54] the turbulence generated in bubble column bubbly flow can be characterized by the composition of the shear turbulence induced by the liquid velocity gradient and the bubble-induced turbulence (BIT).The appearance of − 5/3 scaling on the turbulent kinetic energy spectrum is very likely the effort of liquid phase shear turbulence but homogenised by the bubbles while the occurrence of − 3 scaling is the result of bubble induced turbulence (BIT).It should be noted that in a recent work of LES modelling bubble column bubbly flow [34], a − 25/3 scaling law was identified and the authors attributed it to the contribution from the BIT.As can be seen from Fig. 13(a), there exists an obvious transition in the slope of the predicted energy spectrum E 11 (κ) at around κ B1 ≈ 125m − 1 .For the wavenumber smaller than κ B1 , the − 5/3 scaling which is located in the inertial sub-range is well recovered while for the wavenumber greater than the representative wavenumber corresponding to the bubble size, the slope was found to approach − 3, being consistent with the experimental observations reported previously by other researchers [3,7,40,41,44,50,51,53,59].We believe that this effect can be partially described by consideration of the modified SGS-TDF and SGS-AMS terms in the LES modelling.It is cautiously mentioned that the bubble induced turbulence due to the rising bubble wakes may decay quite rapidly prior to the onset of turbulence spectral transfer.Pope (2000) also noted that the major mechanism in the energy cascade (l < l EI ) is the kinetic energy transfer to successively smaller scales (shear and bubble-induced viscosity dominants) and viscous dissipation (molecular viscosity dominants) where l EI denotes the turbulence length scale between anisotropic large eddies and isotropic tiny eddies.Thus, the bubble induced turbulence energy will not participate in the kinetic energy transferring in larger length scales associated with low wave numbers instead it will contribute to the turbulent kinetic energy transfer from the wave numbers associated with higher eddy frequencies.These eddies that are strongly related to the bubble size contributes to the turbulent dissipation in the higher frequency range as Lance and Bataille (1993) pointed out.
As can be seen from Fig. 13(b), the cut-off length scale in our LES modelling of the species concentration spectrum is slightly larger than the wavenumber of order    found in the predicted species mass fraction spectrum, following a − 5/3 Kolmogorov scaling law in the inertia subrange.In the present work, the species concentration fluctuations of OH − and aqueous CO 2 at z = 0.75H of the centreline of the bubble column are traced in the LES modelling.After applying two-point correlation and Fourier fast transformation, the concentration density spectrum was obtained as shown in Fig. 13(b).
It can be seen from the figure that a clear − 5/3 scaling law is found in the inertia subrange, corresponding to κ < 1/Δ cutoff and an approximate − 1 scaling law is recovered after this transition position.This further supports the argument that the mass transfer occurring in CO 2 chemisorption process in the bubble column is strongly affected by the turbulent eddies acting on the bubble surfaces as schematically indicated in Fig. 1.The wavenumber corresponding to the transition point in the species concentration spectrum right seemed to correspond to the thickness of the liquid film enclosing the bubbles.Another interesting point needs to be noted that a transition in the scaling for E Co2(aq) (κ) takes place at around κ 2 ≈ 100m − 1 , where the wavelength is slightly larger than κ B .This may indicate that a slower chemisorption occurs when the bubbles are entrapped by a similar or relative larger eddies following a fast process when the CO 2 gas penetrates across the bubble surface to enter the surrounding liquid film hit by far too small turbulent eddies.The transition for E OH − (κ) to take place at around finally reaches an equilibrium in the entire domain of the bubble column with the size of the largest eddies being the order of bubble column diameter.

Conclusions
Euler/Euler LES simulations of the carbon dioxide chemisorption in the sodium hydroxide solution in three different size bubble columns with low volume fractions were conducted using three different momentum exchange modelling models, i.e., A-conventional force models; B-modified SGS-TDF model; C-modified SGS-TDF and SGS-AMS models.The influences of considering the terms of SGS-TDF and SGS-AMS on the interfacial mass transfer of chemisorption process were assessed in the LES modelling.The bubble dynamic response to the turbulent eddies was also considered through the modification in the turbulent eddy viscosity and the mass eddy diffusivity.In the meantime, the bubble number density model was also used in the LES modelling of the mass transfer.The LES simulations using the three models have been conducted based on the two rectangular and one cylindrical bubble column experimental cases.The main conclusions reached as the results of the present study can be summarised as follows: (1) The LES simulation for prediction of the evolution of the species concentration and pH-value using the joint models of SGS-TDF and SGS-AMS (Model C) yields a similar degree of agreement with the experimental data.An apparent delay in the chemisorption reaction progress was found when using the models without considering the SGS-AMS modification.A possible explanation for this noticeable difference can be attributed to the fact that the model excluding the SGS-TDF and SGS-AMS (Model A) gives the poor estimation of the bubble lateral dispersion, resulting in an over-prediction of bubble volume fraction in the central region of the bubble column and interfacial mass transfer.The corresponding bubble dynamics in the core region would lead to a shorter bubble residence time, which may reduce the interfacial mass transfer across the phases, causing a delayed consumption on the hydroxide.
(2) The predicted time-averaged axial bubble velocity profiles in the lateral and radial directions for three models were compared with the experimental and the simulations conducted by other researchers.Quantitatively, a good agreement overall with the experimental data with N 2 supplied and with CO 2 supplied was obtained.However, the predicted distribution of the bubble axial velocities by considering the SGS-TDF and SGS-AMS (Model C) in the core region of the bubble columns was found to be lower than the ones using the models (Model A and Model B), which was more consistent with the experimental data, indicating the necessity of considering the SGS-AMS for an accurate description of bubble dispersion and bubble dynamics in the bubble columns.It should be noted that the consistency for the predicted bubble velocity profiles compared with the experimental data become poorer at the near-wall region, which needs further investigation on how the LES modelling with SGS-TDF and SGS-AMS models can implement the Milelli's criterion for mesh requirement.
(3) The cross-sectional averaged and time-averaged equivalent bubble diameter along the axial height of the bubble columns exhibits a stepped reduction trend as measured from the gas sparger.It was found that such predicted a bubble diameter change along the height has a smaller deviation from the experimental data when considering both SGS-TDF and SGS-AMS (Model C) in the LES model However, the bubble size was still overestimated and this may require the effect of the anisotropic SGS-Reynolds stress in the LES model to be considered.
(4) The turbulent kinetic energy spectrum and the concentration spectrum of hydroxide and aqueous CO 2 obtained in the LES modelling still present a typical − 5/3 scaling and − 3 scaling laws for the former while the transition position in the slope was found to be close to the estimated representative bubble wavenumber.This indicates that the bubble induced turbulence only contributes the interfacial mass transfer with those eddy length scale smaller than the equivalent bubble size.For the species concentration spectrum, the typical − 5/3 scaling law was also identified with the eddies falling into the inertial sub-range, following by a scaling which approaches − 1 while the transition cutoff length scale was found to be slightly smaller than the Kolmogorov scale η, which is consistent with the work reported by Lundgren [38].
The mechanism of the CO 2 chemisorption related to the interfacial mass transfer as reflected from the concentration spectrum indicates that the turbulent eddies with the size smaller than the equivalent bubble diameter have a major impact on the interfacial mass transfer, again affirming the importance of SGS-TDF and SGS-AMS.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.height of 2.00 m (Fig. B1).The bubble column is filled with NaOH solution (c NaOH L = 32mmoll − 1 ) with a static liquid height of 1.40 m.Nitrogen was injected with superficial velocity of 0.5 cm/s into the system to provide a steady flow profile, the gas hold-up was controlled at 2.1%.Then the gas changed to carbon dioxide to start the chemical absorption reaction.The gas was introduced in to the column via a capillary gas sparger equipped with 13 needles with an inner diameter of 0.57 mm.In order to obtain the bubble velocity and bubble volume fraction distribution, both PIV (Dantec Dynamics) and high-speed camera (Phantom VEO 1310) were applied.The images are acquired by high-speed camera at a rate of 2000 frames per second and a total scanning period of 10 s.Fig. B2(a) shows a typical image of bubbles captured at the central region of the bubble column.Fig. B2(b) displays the method of calculating the bubble volume fraction based on the images.For example, the second bin in Fig. B2(b) contains 9 bubbles, each individual bubble is assumed to be approximately an ellipsoidal sphere with its long axis and short axis recorded in order (Table B1).By averaging the ratio of the totally bubble occupied area to the area of each observation window for 50 samples, the local bubble volume fraction is estimated.The bubble velocity profile was obtained by evaluating the double frames of the images using the PIV, as shown in Fig. B3.The velocity vectors obtained from the instantaneous images obtained can be used to calculate the mean and fluctuation, which can give out the turbulence characteristics of the bubble flow in the bubble column.The time averaging has taken 1000 frames.

Fig. 1 .
Fig. 1.Schematic of turbulent eddy fluctuations around the bubbles using the LES spatial filtering in bubble column bubbly flows.

Fig. 2 .
Fig. 2. Schematic diagram of the mass transfer between the rising-up bubble and the surrounding liquid phase (NaOH solution) in the bubble column.

Fig. 3 .
Fig. 3. Schematic diagrams of the experimental set-up of the three reactive bubble column reactors.

2 − 3 , HCO − 3 ,
OH − and H 2 O.To initialise the LES simulation, the mass fraction of OH-in liquid phase in the bubble columns was calculated from the given pH while the rest related species mass fractions Y CO2, L , Y HCO − 3 , L and Y CO 2− 3 were set to 1 × 10 − 50 .H 2 O was set as the constrained species obeying ∑ N i=1 Y i,L = 0.

Fig. 4 .
Fig. 4. Time history of the predicted (a) local species concentrations and (b) pH evolution at middle point at z = 980 mm compared with the reported data of Darmana et al. [12] (solid lines).

Fig. 6 .
Fig. 6.Time history of predicted pH-value at middle point at H/D = 1.

Fig. 7 .
Fig. 7. Time-averaged radial distribution of axial bubble velocity at z = 750 mm (a) with N 2 supplied (b) with CO 2 supplied.
the local d 32 is obtained by solving d 32 = (

Fig. 9 .Fig. 10 .
Fig. 9. Variation of cross-sectional averaged and time-averaged equivalent dimensionless bubble along the axial height of the bubble column.

Fig. 13 .
Fig. 13.Concentration spectrums of hydroxide ions and aqueous CO 2 and liquid turbulent axial velocity spectrum obtained at the centreline of the bubble column at z = 0.75H.

Fig. B2 .
Fig. B2.Typical image of bubbles in the (a) central region and (b) right-half region of the bubble column.The focal length used in the experiment was 94 mm.

Fig. B3 .
Fig. B3.Instantaneous velocity vectors of liquid flow field in bubble column, obtained by PIV.

Table 1
Geometry parameters and the experimental conditions of the three bubble columns.

Table B1
Details of bubble measurement based on the images of bubbles captured.