Abstract
Large eddy simulations (LES) of dispersed gas-liquid flows for the prediction of flow patterns and its applications have been reviewed. The published literature in the last ten years has been analysed on a coherent basis, and the present status has been brought out for the LES Euler-Euler and Euler-Lagrange approaches. Finally, recommendations for the use of LES in dispersed gas liquid flows have been made.
1. Introduction
Gas-liquid flows are often encountered in the chemical process industry, but also numerous examples can be found in petroleum, pharmaceutical, agricultural, biochemical, food, electronic, and power-generation industries. The modelling of gas-liquid flows and their dynamics has become increasingly important in these areas, in order to predict flow behaviour with greater accuracy and reliability. There are two main flow regimes in gas-liquid flows: separated (e.g., annular flow in vertical pipes, stratified flow in horizontal pipes) and dispersed flow (e.g., droplets or bubbles in liquid). In this work, we consider only dispersed bubbly flows.
Dispersed Bubbly Flow. The description of bubbly flows involves modelling of a deformable (gas-liquid) interface separating the phases; discontinuities of properties across the phase interface; the exchange between the phase; and turbulence modelling. Most of the dispersed flow models are based on the concept of a domain in the static (Eulerian) reference frame for description of the continuous phase, with addition of a reference frame for the description of the dispersed phase. The dispersed phase may be described in the same static reference frame as the continuous, leading to the Eulerian-Eulerian (E-E) approach or in a dynamic (Lagrangian) reference frame, leading to the Eulerian-Lagrangian (E-L) approach.
In the E-L approach, the continuous liquid phase is modelled using an Eulerian approach and the dispersed gas phase is treated in a Lagrangian way; that is, the individual bubbles in the system are tracked by solving Newton’s second law, while accounting for the forces acting on the bubbles. An advantage here is the possibility to model each individual bubble, also incorporating bubble coalescence and breakup directly. Since each bubble path can be calculated accurately within the control volume, no numerical diffusion is introduced into the dispersed phase computation. However, a disadvantage is, the larger the system gets the more equations need to be solved, that is, one for every bubble.
The E-E approach describes both phases as two continuous fluids, each occupying the entire domain, and interpenetrating each other. The conservation equations are solved for each phase together with interphase exchange terms. The E-E approach can suffer from numerical diffusion. However, with the aid of higher order discretization schemes, the numerical diffusion can be reduced sufficiently and can offer the same order of accuracy as with E-L approach (Sokolichin et al. [1]). The advantage here is that the computational demands are far lower compared to the E-L approach, particularly for systems with higher dispersed void fractions. We review these approaches here with respect to the turbulence descriptions.
Turbulence Modelling. The major difficulty in modelling multiphase turbulence is the wide range of length and time scales on which turbulent mixing occurs. The largest eddies are typically comparable in size to the characteristic length of the mean flow. The smallest scales are responsible for the dissipation of turbulence kinetic energy. The Direct Numerical Simulation (DNS) approach, with no modelling, resolves all the scales present in turbulence. However, it is not feasible for practical engineering problems involving high Reynolds number flows. The Reynolds-Averaged Navier–Strokes (RANS) approach is more feasible; it models the time-averaged velocity field either by using turbulent viscosity or by modelling the Reynolds stresses directly.
The large eddy simulation (LES) falls between DNS and RANS in terms of the fraction of the resolved scales. In LES, large eddies are resolved directly, that is, on a numerical grid, while small, unresolved eddies are modelled. The principle behind LES is justified by the fact that the larger eddies, because of their size and strength, carry most of the flow energy (typically 90%) while being responsible for most of the transport, and therefore they should be simulated precisely (i.e., resolved). On the other hand, the small eddies have relatively little influence on the mean flow and thus can be approximated (i.e., modelled). This approach to turbulence modelling also allows a significant decrease in the computational cost over direct simulation and captures more dynamics than a simple RANS model.
In RANS models often the assumption of isotropic turbulence is made for the core of the flow, which is not valid in dispersed bubbly flows; that is, the velocity fluctuations in the gravity direction are typically twice those in the other directions. This assumption is not made in LES for large structures of the flow, giving LES an advantage over RANS for the core regions of the flow. However, the situation is different close to the walls, where LES’ assumption of isotropic turbulence is heavily violated, due to the absence of large eddies close to the walls.
2. LES for Dispersed Bubbly Flows
In dispersed bubbly flows, the large-scale turbulent structures interact with bubbles and are responsible for the macroscopic bubble motion, whereas small-scale turbulent structures only affect small-scale bubble oscillations. Since, large scales (carrying most of the energy) are explicitly captured in LES and the less energetic small scales are modelled using a subgrid-scale (SGS) model, LES can reasonably reproduce the statistics of the bubble-induced velocity fluctuations in the liquid.
There are three important considerations for modelling of dispersed bubbly flows.(1)Separation of length scales of the interface, that is, micro-, meso-, and macroscales. The separation of these scales forms the basis for “filtering” the Navier–Stokes equations and applying proper model equations for multiphase situation. Important for dispersed flow is to identify the scales at which the governing equations are to be applied; microscales, that is, scales which are small enough to describe individual bubble shapes; mesoscales, which are comparable to bubble sizes; and macroscales, which entail enough bubbles for statistical representation. (2)The grid-scale equations. Depending on the ratio of the length scales introduced above, with the grid resolution we can afford, on a given computer hardware, a proper form of the governing equations must be chosen. For instance, if the mesh size is in the micro-scale order, one can use single-fluid, interface tracking techniques to solve the problem. If, on the other hand, the grid size is large enough for statistical description of bubbles, the E-E approach can be used. Should the grid size be comparable to the meso-scales, we are in a limiting area for both approaches, and special care must be taken in order to solve equations which describe the underlying physics consistently. (3)The physical models. Depending on the selected grid-scale equations, physical models of various complexities must be employed. The options here are numerous, whether they concern turbulence modelling or interphase modelling, but these models are generally simpler in case more of the microscales are resolved.
In the following sections, we describe each of these three elements to model turbulent dispersed bubbly flow.
2.1. Filtering Operation
The aim of filtering the Navier-Stokes equations is to separate the resolved scales from the SGS (nonresolved). The interface between the phases, and the level of detail required in its resolution/modelling, defines the filter in a multiphase flow.
When LES is applied at a micro-scale, filtering of turbulent fluctuations needs to be combined with interface tracking methods. These methods have been developed and used in both dispersed flow and free surface flow by Bois et al. [2], Toutant et al. [3, 4], Magdeleine et al. [5], Lakehal [6], and Lakehal et al. [7]. These methods require that all phenomena having an influence on space and time position of the interface are also simulated. For the amount of details required and the large size of practical problems of interest, these types of models should merely be seen as a support for the modelling and validation of more macroscopic approaches and cannot address a real industrial-scale problem (Bestion [8]).
When LES is applied at a macro-scale, the interface resolution is not considered. However, in practical simulations, these would require too coarse grids, leading to poor resolution of turbulence quantities. Much more often we are in the meso-scale region, in which the mesh size is comparable to bubble sizes. This pushes the main assumptions of the E-E approach to its limit of validity, and the grid is not fine enough for full interface tracking. In other words, the mesh requirement for E-E multiphase modelling conflicts with the requirements by LES approaches [9].
The issue of the requirement of the mesh size was first addressed by Milelli et al. [11] who carried out a systematic analysis and performed a parametric study with different mesh sizes and bubble diameters. They showed that for case of a shear layer laden with bubbles it was possible to provide an optimum filter width , where is the filter width and is the bubble diameter (shown in Figure 1). This means that the grid space should be at least 50% larger than the bubble diameter. The constraint imposed on the ratio implies that the interaction of bubbles with the smallest resolved scales is captured without additional approximation.
2.2. Grid-Scale Equations
The principle of the LES formulation is to decompose the instantaneous flow field into large-scale and small-scale components via a filtering operation. If denotes the filtered or grid-scale component of the variable that represents the large-scale motion then where is the variable of interest, subscript refers either to the liquid or the gas phase. In the remainder of this paper, we omit the bars of all resolved variables for the sake of simplicity. The following filtered equations are obtained:
The right hand side terms of (3) are, respectively, the stress, the pressure gradient, gravity, and the momentum exchange between the phases due to interface forces.
The SGS stress tensor which reflects the effect of the unresolved scales on the resolved scales is modelled as where is the effective viscosity.
In the E-E approach, separate equations are required for each phase (see (3), , ), together with interphase exchange terms (for details, Drew [12]). In most of the investigations, turbulence is taken into consideration for the continuous phase by SGS models. The dispersed gas phase is modelled as laminar, but influence of the turbulence in the continuous phase is considered by a bubble-induced turbulence (BIT) model.
In the E-L approach, there are two coupled parts: a part dealing with the liquid phase motion and a part describing the bubbles motion. The dynamics of the liquid are described in a similar way as in the E-E approach, whereas the bubble motion is modelled through the second law of Newton.
Since, the governing equations for the liquid and gas phase are expressed in the Eulerian and Lagrangian reference frames, respectively; a mapping technique is used to exchange interphase coupling quantities. Depending upon the volume fraction of the dispersed phase, one-way (e.g., ) or two-way coupling between gas phase to liquid phase ( prevails. In both cases, bubble-bubble interactions (i.e., collisions) can be neglected, but the effect of the bubbles on the turbulence structure in the continuous phase has to be considered for higher volume fraction and does not play any role in lower volume fraction of gas phase Elgobashi [13]. The work reviewed here considers the two-way coupling which consists of the following.
2.2.1. Forward Coupling (Liquid to Bubble)
In the forward coupling, calculated liquid velocities, velocity gradients, and pressure gradients on an Eulerian grid are interpolated to discrete bubble locations for solving the Lagrangian bubble equation motion.
2.2.2. Backward or Reversed Coupling (Bubble to Liquid)
The forces available at each bubble’s centroid need to be mapped back to the Eulerian grid nodes in order to evaluate the reaction force . The two-way interaction (forward and backward) is accomplished with a mapping method, for example, PSI-cell method [14], modified PSI-wall-method [15], or mapping functions discussed by Deen et al. [16].
2.3. Interfacial Forces
The motion of a single bubble with constant mass can be written according to Newton’s second law:
The bubble dynamics are described by incorporating all relevant forces acting on a bubble rising in a liquid. It is assumed that the total force, , is composed of separate and uncoupled contributions originating from pressure, gravity, drag, lift, virtual mass, wall lubrication and wall deformation turbulent dispersion:
For each force the analytical expression or a semiempirical model is used, based on bubble behaviour observed in experiment or in DNS.
To summarize, the influence/contribution of these forces are as follows.(1)The modeling of the lift force for capturing bubble plume meandering and bubble dispersion is important. However there is an uncertainty regarding appropriate value or correlation representing lift coefficient. There is also recommendation that bubble size-dependent lift coefficient should be chosen [17].(2)The value of the lift coefficient can be different than the one used in RANS approach. It is because of different handling of factors responsible for bubble dispersion, that is, the interaction between the bubbles and influence of turbulent eddies in the liquid phase. In RANS approach, they are considered by means of the lift and turbulent dispersion force, with uncertainty of exact contribution of the individual forces. Most of the investigators use a constant value of the lift coefficient (), while the value of the turbulent dispersion coefficient is varied (0.1 to 1.0) to get good agreement with the experimental data. However, in LES, bubble dispersion caused by liquid phase turbulent eddies is implicitly calculated, and a more realistic contribution of the lift force can be used. The coefficient for the effective lift force thus may vary between the two approaches [18].(3)The virtual mass force is proportional to the relative acceleration between the phases and is negligible once a pseudosteady state is reached. It has little influence on the simulation results for bubble plumes [19], Milelli [20]. It is mainly because of the acceleration and deceleration effects are restricted to small end regions of the column. A constant coefficient is used in almost all investigations. (4)In LES, through filtering, velocities are decomposed into a resolved and a SGS part. The resolved part of the turbulent dispersion is implicitly computed. However, in case of a bubble size smaller than the filter size, turbulent transport can be present at SGS level and should be considered [9]. This can be done using a one-equation model, wherein it can be modelled by replacing the total kinetic energy by SGS contribution (). By the same argument, other forces also need modelling at SGS level.
The values or expressions for the coefficient of drag, lift and virtual mass force used by different investigators are given in Tables 1 and 3.
2.4. SGS Models
It is well known that in turbulent flow energy generally cascades from large to small scales. The primary task of the SGS model therefore is to ensure that the energy drain in the LES is same as obtained with the cascade fully resolved as one would have in a DNS. The cascading, however, is an average process. Locally and instantaneously the transfer of energy can be much larger or much smaller than the average and can also occur in the opposite direction (“backscatter”).
2.4.1. Smagorinsky [21] Model
The simplest, well-known, and mostly used Smagorinsky [21] model is based on the Boussinesq hypothesis. It requires the definition of time and length scales and a model constant. Smagorinsky used the following expression to calculate the turbulent viscosity, that is, the SGS viscosity: where is the (laminar) dynamic viscosity, is the Smagorinsky constant, is the characteristic strain tensor of filtered velocity, and Δ is the filter width, usually taken as the cubic root of the cell volume.
In the single-phase flow literature, the value of the constant used is in the range from = 0.065 (Moin and Kim [22]) to = 0.25 (Jones and Wille [23]). The value of used in gas-liquid flows varies from that of single phase flow and is in the range of 0.08 to 0.12 [11, 20, 24, 25]. The lower range of value, compared to single phase, could be attributed to the interphase coupling term, which acts as a form of SGS model and can make contribution to the turbulent kinetic energy dissipation. The sensitivity analysis carried out for value shows that larger values can produce excessive damping effect to the liquid velocity field and eventually leads to a steady-state solution [26, 27].
The main reason for the frequent use of the Smagorinsky model is its simplicity. Its drawbacks are that the constant has to be calibrated and its optimal value may vary with the type of flow or the discretization scheme. Moreover, the model is purely dissipative and hence does not account either for the small-scale effect on the large scales adequately (by neglecting the “backscatter” of turbulent energy), while it acts purely as a drain for the turbulent kinetic energy.
The dynamic model, originally proposed by Germano et al. [28], eliminates some of these disadvantages by calculating the Smagorinsky constant as a function of space and time from the smallest scales of the resolved motion.
2.4.2. Dynamic SGS Model
The dynamic SGS model assumes SGS turbulent energy to be in local equilibrium (i.e., production = dissipation). The eddy viscosity is estimated from (7) but with a as a local, time-dependent variable.
The basic idea is to apply a second test filter to the equations. The new filter width, twice the size of the grid filter, produces a resolved flow field. The difference between the two resolved fields is the contribution of the small scales whose size is in between the grid filter and the test filter. The information related to these scales is used to compute the model constant. The advantage here is that no empirical constant is needed and that the procedure allows the negative turbulent viscosity implying energy transfer from smaller to larger scales (energy back-scatter). This effect, in principle, allows both an enhancement and attenuation of the turbulent intensity introduced by the bubbles.
The model has a few drawbacks; wide fluctuations in dynamically computed constants can cause stability issues, along with additional computational expense.
2.4.3. One-Equation Model
In spite of the fact that dynamic SGS model calculates model constant , thus making a constant-free model, it lacks the information on the amount of SGS turbulent kinetic energy, a datum which may prove useful in modelling some aspects of dispersed flows (e.g., SGS bubble-induced turbulence).
The essence of the one-equation model is to solve additional transport equation for SGS turbulent kinetic energy: Here, is production of SGS turbulent kinetic energy and is defined as and SGS viscosity is obtained from
The availability of the SGS turbulent kinetic energy allows for modelling of SGS interphase sorces such as bubble-induced turbulence and turbulent dispersion at SGS. The application of one-equation SGS model for bubbly flows is illustrated in more detail in sections below.
2.5. Effect of Bubble-Induced Turbulence (BIT)
In the E-E approach, the turbulent stress in the liquid phase is considered to have two contributions, one due to the inherent, that is, shear-induced turbulence that is assumed to be independent of the relative motion of bubbles and liquid and the other due to the additional bubble-induced turbulence (Sato and Sekoguchi [37]). For BIT there are two modelling approaches. The first approach is proposed by Sato and Sekoguchi [37] and Sato et al. [42]: with as a model constant which is equal to 0.6 and as the bubble diameter. Milelli et al. [11, 24] found that the modelling of the bubble-induced turbulence did not improve the results. They tried two different formulations: the Tran model and the Sato model and found that they have negligible effect. This was attributed to fact that the bubble-induced viscosity (and turbulence) is not crucial, the turbulence being mainly driven by the liquid shear, and a low void fraction (≈2% leading to ) did not significantly modify the situation. It was thought that in a case in which the bubbles actually drive the turbulence (via buoyancy and/or added mass forces), the situation would be different. However, in subsequent studies, similar observations were made in bubble plumes simulated by Deen et al. [19], Dhotre et al. [20], Ničeno et al. [9].
The second approach for the modelling of BIT allows for the advective and diffusive transport of turbulent kinetic energy. This model incorporates the influence of the gas bubbles in the turbulence by means of additional source terms in the equation and is taken to be proportional to the product of the drag force and the slip velocity between the two phases. This approach was used in work of Niceno et al. [10] through the use of a one-equation model. They found significant influence of the additional source terms as used by Pfleger et al. [43], as shown in Figure 2.
(a)
(b)
Figure 2 shows the comparison of the liquid kinetic energy obtained for the case of a bubble plume rising in a square column. It can be seen that the simulation without BIT underpredicts the turbulent kinetic energy. The use of the Sato model reproduced the double-peaked profile for kinetic energy. The Pfleger model also reproduced the experimental data very well. Figure 2(b) shows the ratio of the modelled SGS energy to the resolved energy. With no BIT, this ratio has the lowest value, whereas the Sato model yields more SGS energy, while the Pfleger model gives a ratio that is roughly twice as high, which is particularly pronounced in the middle of the column. Table 2 gives a summary of BIT models proposed by various investigators.
3. Numerical Details
Crucial parameters for obtaining reliable LES results are the time step selection, the total time for gathering good statistics of the averaged variables, and discretization schemes for the variables. The time step choice is determined by the criterion that the maximum Courant-Fredrichs-Levy (CFL) number must be less than one ().
For flow variables, central difference should be used for discretization of advection terms and avoid using diffusive upwind schemes. However, for scalars variables, high-order schemes (MUSCL, QUICK, or Second-Order) may be tolerable to avoid nonphysical solutions (e.g., negative volume fractions). An alternative to high-order schemes are the bounded central differences. The risk with use of all but central scheme is their diffusivity. Their influence on LES may exceed the modelled SGS transport.
It is necessary to follow the initial phase of the simulation, wherein the turbulent strutures develop starting from initial condition and to reach a statistiacally steady state. The duration of this phase depends on the flow characteristics. The simulation must be run for a total time long enough to allow all turbulent instabilities that develop during this phase to be convected across the region of interest. However, the convecting velocities of the turbulent structures and the regions of interest are not always known as a priori. This is why it is recommended to run the simulation a multitude (typically 5 times) of the slowest integral time scales, which often is the flow through time defined as the ratio of the system height over the bulk (superficial) velocity.
4. LES Prediction of the Flow Pattern for Dispersed Bubbly Flows
Here, we review different LES studies that were performed using the E-E and E-L approaches for simulating flow patterns in gas-liquid bubbly flows. Table 1 gives a summary of key numerical parameters (filter size, number of grids, SGS model, bubble diameter, coefficient for interfacial forces) and experimental details (geometrical dimension, sparger design, range of superficial gas velocity) used by investigators.
4.1. Euler-Eulerian (E-E) Approach
4.1.1. Milelli et al. [11, 24, 49]
Milelli et al. reported for the first time two-phase LES with E-E approach. They first investigated statistically 2D flow configuration and then free bubble plume.
They addressed important concerns related to the two-phase LES simulation. For instance, they found that the optimum ratio of the cutoff filter width (i.e., the grid) to the bubble diameter () should be around 1.5. That means mesh size should be at least 50% larger than the bubble diameter (Figure 1) so that (a) bubble size determines the largest scale modelled (b) and its interaction with the smallest calculated scale above the cut-off is captured. This is also supported by the scale-similarity principle of Bardina et al. [50].
Milelli [49] investigated LES for a free bubble plume and compared their predictions with the experiment of Anagbo and Brimacombe [51]. Here, they found that the mean quantities were not strongly affected by the different SGS models. Moreover they found little impact of the dispersed phase on the liquid turbulence, from the turbulent energy spectrum taken in the bubbly flow region which revealed a power-law distribution oscillating between −5/3 and −8/3 in the inertial subrange. The results conform to previous studies, which attributed the more dissipative spectrum to the presence of the dispersed phase. Hence, they found no influence of modifying the SGS model to account for bubble-induced dissipation.
Further, they observed in simulation that the lift coefficient value plays a major role in capturing the plume spreading and the used lift coefficient may differ for an LES compared to the one that is justified in an RANS approach. The plausible explanation here is from different handling of two factors responsible for bubble dispersion, that is, interaction between the bubbles and influence of turbulent eddies in the liquid phase.
4.1.2. Deen et al. [19]
Deen et al. [19] reported LES for gas-liquid flow in a square cross-sectional bubble column for the first time. They investigated the performance of RANS and LES approaches, influence of the interphase forces, and bubble-induced turbulence.
They found that RANS approach (k-ε model) overestimated the turbulent viscosity and could only predict low frequency unsteady flow. On other hand, LES as shown in Figure 3 reproduced high frequency experimental data and predicted the strong transient bubble plume movements as in an experiment.
Furthermore, they also identified that the lift force is responsible for transient spreading of the bubble plume and in absence of it, only with drag force, the bubble plume showed no transverse spreading.
They considered the effective viscosity of the liquid phase with three contributions: the molecular, shear-induced turbulent (modelled using Smagorinsky model), and bubble-induced turbulent viscosities [37]. Like in the work of Milelli, they confirmed the marginal effect of the BIT on the predictions. The effect of virtual mass force on the simulated results was also found to be negligible.
4.1.3. Bove et al. [29]
Bove et al. [29] reported LES with E-E approach for the same square cross-sectional bubble column as used by Deen et al. [19]. They studied the influence of numerical modelling of the advection terms and the inlet conditions on LES performance. The upwind first-order and higher-order Flux Corrected Transport (FCT) schemes for both the phase fraction equations and the momentum equations were employed. The simulations using a second-order FCT scheme showed relatively good agreement with the measurement data of Deen et al. [19]. The authors showed that the proper discretization of the momentum and volume fraction equations is essential for correct prediction of the flow field.
Further, the LES results were found to be very sensitive to inlet boundary conditions (Figure 4). Three different inlet configurations simulated showed that the inlet modelling influences the predicted fluid flow velocity (as in Figure 4(a)) and an important fluid flow parameter, the turbulent viscosity (Figure 4(b)). In this work, the sparger (a perforated plate) was not modelled due to the difficulty in adapting the mesh grid to the geometry. They also suggested that near wall region description in the SGS models is important, and the lack of the near wall modelling can lead to erroneous prediction of frictional stresses at the wall.
(a)
(b)
They used drag model for the contaminated water which gave a better prediction of the slip velocity; however, the velocity profile was underestimated for both gas and liquid phase. Reason for the underprediction was not clear, whether it was due to drag model or an improper value of the lift coefficient used or an error in the near wall modelling. Need for further work in this direction was suggested.
4.1.4. Zhang et al. [26]
Zhang et al. [26] reported LES in a square cross-sectional bubble column. They investigated the Smagorinsky model constant and carried out a sensitivity analysis. It was found that higher values led to higher effective viscosity which dampens the bubble plume dynamics leading to a steep mean velocity profile (as shown in Figure 5). They obtained a good agreement with the measurements with in range of 0.08–0.10. They also confirmed that the lift force plays a critical role for capturing the dynamic behaviour of the bubble plume.
(a) Axial liquid velocity
(b) Axial gas velocity
They extended the work of Deen et al. [19] and predicted the dynamic behaviour in the square bubble column using a k-ε turbulence model extended with BIT.
4.1.5. Tabib et al. [17]
Tabib et al. [17] reported LES using E-E approach in a cylindrical column for a wide range of superficial gas velocity. In accordance with the earlier work, they confirmed the importance of a suitable lift coefficient and drag law. Moreover, they studied the influence of different spargers (perforated plate, sintered plate, and single hole) and turbulence models (k-ε, RSM, and LES) using the experimental data of Bhole et al. [52]. The main findings from the study were that the RSM performs better than the k-ε model; the LES was successful in predicting the averaged flow behaviour and was able to simulate the instantaneous vortical-spiral flow regime in the case of a sieve plate column, as well as the bubble plume dynamics in case of single-hole sparger. Finally, they concluded that LES can be effectively used for the study of the flow structures and instantaneous flow profiles.
4.1.6. Dhotre et al. [20]
Dhotre et al. [20] reported LES with an E-E approach for a gas-liquid flow in a square cross-sectional bubble column. They studied the influence of SGS models: Smagorinky and Dynamic models of Germano et al. [28]. It was found that both the Smagorinsky model ( = 0.12) and the Germano model predictions compared well with the measurements.
They further investigated the value of obtained from the Germano model. Reason for similar performance of both models was clear from the probability density function of (from Germano model) over the entire column. As shown in Figure 6, the value of has the highest probability in the range of 0.12–0.13. Like Zhang et al. [26], the authors confirmed that with a proper BIT model, RANS also performed well for mean quantities of flow variables. Figure 7 shows the comparison of the predicted instantaneous vector flow field for axial liquid velocity from all the three models (Smagorinky, Germano and RANS).
(a) Germano
(b) Smagorinsky
(c) RANS
It was further concluded that the Germano model can give correct estimates for the configuration under consideration and, in general, can be used for other systems where is not known as “a priori” from previous analysis.
4.1.7. Niceno et al. [10]
Niceno et al. [10] investigated LES with E-E approach for a gas-liquid flow in a square cross-sectional bubble column. They demonstrated the applicability of a one-equation model for the SGS kinetic energy (). The predictions showed that the one-equation SGS model gives superior results to the Germano model with the additional benefit of having information on the modelled SGS kinetic energy: with = 0.07 a model constant. They studied the influence of two approaches for bubble-induced turbulence: approach of an algebraic model (Sato et al. 1975) and extra source terms (as used in Pflger et al. 1999) in the transport equation for SGS kinetic energy approach. It was found that the latter approach improved the quantitative prediction of the turbulent kinetic energy (as shown in Figure 2(a)). The modelled SGS kinetic energy for the Pfleger model found to be much higher than for the Sato model (Figure 2(b)), indicating the Pfleger model needs a more appropriate constant for LES.
They suggested that the modelled SGS information can be used to access the SGS interfacial forces, in particular the turbulent dispersion force. In their work, the effect of SGS turbulent dispersion force could not be determined as the bubble size was almost equivalent to the mesh size.
4.1.8. Dhotre et al. [18]
Dhotre et al. [18] extended LES with E-E approach for a gas-liquid flow in a large-scale bubble plume. The predictions at three elevations were compared with the measurement data of Simiano [55] and an RANS prediction. The LES approach was shown superior in capturing the transient behaviour of the plume (Figure 8) and predicts second-order statistics of the liquid phase accurately.
(a)
(b)
They emphasized the crucial role of the lift force in the prediction of the lateral behaviour of the bubble plumes. In the RANS approach the turbulent dispersion force is required to reproduce the bubble dispersion; however, in LES, bubble dispersion is implicitly calculated by resolving the large-scale turbulent motion responsible for bubble dispersion. The dependence of the bubble dispersion with the value of lift coefficient was also observed in Milelli et al. [11, 24], Deen et al. [19], Lain and Sommerfeld [56], Van den Hengel et al. [31], Tabib et al. (2008), and Dhotre et al. [20]).
Dhotre et al. [18] found good agreement with the measurement data at higher elevation, while discrepancies were observed at lower elevation, near the injector. The reason for the discrepancies was attributed to the absence of modelling bubble coalescence and breakup. This was also found in the work of Van den Hengel et al. [31], wherein the authors showed that most of the coalescence occurs in the lower part of the column and recommended to consider bubble size distribution and coalescence and breakup models for reproducing the bubble behaviour near the sparger.
4.1.9. Niceno et al. [10]
Niceno et al. [10] reported LES with E-E approach for a gas-liquid flow in a square cross-sectional bubble column. They compared two different codes (CFX-4 and Neptune) and two subgrid-scale models (as in Figure 9). The prediction from the Smagorinsky model in the Neptune CFD code and the one-equation model of CFX-4 was compared with the measurement data of Deen et al. [19]. Agreement between the predictions from the two SGS models was found to be good, and it was concluded that the influence of the SGS model was small. This is in contradiction with earlier work of Van den Hengel et al. [31], where they showed significant contribution of the SGS model (Figure 10), which is discussed in more detail in section (4.2). It remains to be seen if this was due to the fine mesh used by the authors (). Niceno et al. [10] argued that with the known flow pattern in a bubble column, that is, a dominant bubble plume meandering between the confining walls, the biggest eddy having most energy is of the size of the domain cross section. Thus, the grid used in their work was a compromise between sufficiently fine to capture the most energetic eddies, and sufficiently coarse to stay close to the Milelli criterion [11, 24]. Furthermore, they pointed out the limitations of LES with E-L or E-E approach without resolving interface; they indicated that the most influential interfacial forces (drag and lift) are modelled for the large-scale field and their effect from the small scale remains a question. On the other hand, they recommend large-scale simulation, as in the works of Lakehal et al. [25], which explicitly resolves the large-scale part of the interfacial forces and models the part at the SGS level, where the effects are smaller and hence less influential on the accuracy of the results.
(a)
(b)
(c)
4.1.10. Tabib and Schwarz [30]
Tabib and Schwarz [30] extended the work of Niceno et al. [9] and attempted to quantify the effect of SGS turbulent dispersion force for different particle systems, where the particle sizes would be smaller than the filter size. They used LES with E-E approach.
They used the formulation of Lopez de Bertodano [57] to approximate the turbulent diffusion of the bubbles by the SGS liquid eddies for a gas-liquid bubble column system [17]. The bubble size was in range of 3–5 mm. The mesh used in simulations was coarser than the bubble diameter. They found a high contribution from the SGS turbulent dispersion force, when compared with the magnitude of the other interfacial forces (like drag force, lift force, resolved turbulent dispersion force, and force due to momentum advection and pressure). Finally, Tabib and Schwarz concluded that for LES with E-E approach, when the mesh size is bigger than bubble size, the SGS turbulent dispersion force should be used, and a one-equation SGS-TKE model overcomes a conceptual drawback of E-E LES model.
4.2. Euler-Lagrangian (E-L) Approach
4.2.1. Van den Hengel et al. [31]
Van den Hengel et al. [31] reported LES with E-L approach for a gas-liquid flow in a square cross-sectional bubble column. The liquid phase was computed using LES, and a Lagrangian approach was used for the dispersed phase. They used a discrete bubble model (DBM) originally developed by Delnoij et al. [58, 59] and extended it to incorporate models describing bubble breakup and coalescence. The mean and fluctuating velocities predicted in the simulations showed a good agreement with the experimental data of Deen et al. [19].
Authors studied the influence of the SGS model on the predictions and found that without SGS model, the average liquid velocity and liquid velocity fluctuations are much lower compared to the case with a SGS model. This was due to the lower effective viscosity in this case, which led to less dampening of the bubble plume dynamics and subsequently to flatter mean liquid velocity profiles (as shown in Figure 10).
In this work also, the authors confirmed the important role of the lift coefficient in capturing the plume dynamics. They considered two lift coefficients ( = 0.5 and 0.3) and found that a smaller value of the lift coefficient led to higher average velocity and velocity fluctuations and less spreading of the plume, which resulted in overprediction of the average velocity in the centre of the column.
4.2.2. Hu and Celik [32]
Hu and Celik [32] studied LES with an E-L approach for the gas-liquid flow in a flat bubble column. The liquid phase was computed using LES, and a Lagrangian approach was used for the dispersed phase. The authors developed a mapping technique called particle-source-in-ball (PSI-ball) for coupling the Eulerian and Lagrangian reference frames. The concept is a generalization of the conventional particle-source-in-cell (PSI-cell) method as well as a template-function-based treatment [14].
They reported second-order statistics of the pseudo-turbulent fluctuations and demonstrated that a single-phase LES along with a point-volume treatment of the dispersed phase could serve as a viable closure model.
Hu and Celik reported that the predicted mean quantities (such as mean liquid velocity field) were in good agreement with the experimental data of Sokolichin and Eigenberger [54], as shown in Figure 11, and further gave an accurate prediction of the instantaneous flow features, including liquid velocity fluctuations and unsteady bubble dispersion pattern. Hu and Celik also studied the influence of the Smagorinsky constant and found that the constant for multiphase systems falls in a relatively smaller range than for single-phase flows. Higher values of the showed an excessive damping effect to the liquid field, which led to a steady-state solution. This observation is in accordance with other investigators [26, 31]. Furthermore, authors proposed to use as a modeling parameter rather than a phyiscal constant, as the interphase coupling terms used as well as the high frequency turbulent fluctuations contribute to the turbulent kinetic energy dissipation.
(a)
(b)
(c)
(d)
4.2.3. Lain [27]
Lain [27] reported an LES with E-L approach for a gas-liquid flow in a cylindrical bubble column. He used LES for the liquid phase, and a Lagrangian approach for the dispersed gas phase. The interaction terms between liquid and gas phases was calculated using the particle-source-in-cell (PSI-cell) approximation of Crowe et al. [14]. The bubbles were considered as a local source of momentum, and source term was added.
A simple model for the subgrid liquid fluctuating velocity to account for the BIT considered in this work was found to have no influence on the predictions. As in previous works, authors confirmed a strong dependency of the bubble dispersion in the column on the value of transverse lift force coefficient used. He concluded that the lift coefficient depends on the bubble-liquid relative velocity and was the main mechanism responsible for the spreading of bubbles across the column crosssection. He further compared the simulation results with particle image velocimetry (PIV) measurements (Border and Sommerfeld [60]) and k-ε calculations.
4.2.4. Darmana et al. [33]
Darmana et al. [33] used the LES with E-L approach for simulating the gas-liquid flow in a flat bubble column and validated the model with experimental data of Harteveld et al. [61]. They investigated seven sparger designs and their influence on the flow structure. It was found that the model captures the influence of different gas sparging very well (e.g., Figure 12 shows one such case simulated). However, in all cases simulated, authors found systematic overprediction of dispersed phase distribution (25%), which was attributed to an inaccuracy of the drag force and the turbulence model at high gas void fractions.
(a)
(b)
4.2.5. Sungkorn et al. [34]
Sungkorn et al. [34] reported LES with the E-L approach for a gas-liquid flow in a square cross-sectional bubble column. They modelled the continuous liquid phase using a lattice-Boltzmann (LB) scheme, and a Lagrangian approach was used for the dispersed phase. For the bubble phase, the Langevin equation model [41] was used for estimating the effect of turbulence. The bubble collisions were described by a stochastic interparticle collision model based on the kinetic theory developed by Sommerfeld [40]. The predictions showed a very good agreement with the experimental data for the mean and fluctuating velocity components. Figure 13 shows the sanpshots of predicted the bubble dispersion patterns.
(a)
(b)
(c)
(d)
It was also found that their collision model leads to two benefits: the computing time is dramatically reduced compared to the direct collision method and secondly it also provides an excellent computational efficiency on parallel platforms. Sungkorn et al. [34] claim that the methodology can be applied to a wide range of problems. The investigations are valid for lower global void fraction, and further work is required to consider it for higher void fraction systems.
5. Application of LES
5.1. Preamble
The investigations discussed in earlier sections dealt with the use of LES for predicting the flow patterns. In the published literature, the knowledge of flow pattern has been employed for the estimation of equipment performance such as mixing (Joshi and Sharma [62], Joshi [63], Ranade and Joshi [64], Ranade et al. [65], and Kumaresan and Joshi [66]), heat transfer (Joshi et al. [67], Dhotre and Joshi [68]), Sparger design (Dhotre et al. [69], Kulkarni et al. [70]), gas induction (Joshi and Sharma [71], Murthy et al. [72]), and solid suspension (Raghava Rao et al. [73], Rewatkar et al. [74], and Murthy et al. [75]). Joshi and Ranade [76] have discussed the perspective of computational fluid dynamics (CFD) in designing process equipment with their views on expectations, current status, and path forward. The LES simulations provide substantially improved understanding of the flow pattern. Therefore, in this section, the application of LES for design objectives like mixing, heat transfer, and chemical reactions by some investigators will be reviewed. The LES simulations have also been used in the identification of turbulent structures, their dynamics, and the role of structure dynamics in the estimation of design parameters. The LES simulations have also been used in the estimations of terms in k-ε and RSM models such as generation, dissipation and transport of turbulent kinetic energy (k), the turbulent energy dissipation rate (ε), and Reynolds stresses. These estimations have improved the understanding of RANS (k-ε and RSM) models. These two applications of LES are also described briefly.
5.2. Mass Transfer and Chemical Reaction
5.2.1. Darmana et al. [77, 78]
Darmana et al. [77, 78] used LES with E-L approach to simulate flow, mass transfer, and chemical reaction in flat bubble column. They considered mass transfer, rate in liquid-phase momentum equation and reaction interfacial forces in the bubble motion equation.
Also, the presence of various chemical species was accounted through a transport equation for each species. Darmana et al. estimated the mass transfer rate from the information of the individual bubbles directly. They used the model to simulate the reversible two-step reactions found in the chemisorption process of CO2 in an aqueous NaOH solution in a lab-scale pseudo-2D bubble column reactor (e.g., Figure 14). They found good agreement between simulation and measurement for the case without mass transfer. In absence of an accurate mass transfer closure, the authors found that the overall mass transfer rate was lower compared to the measurement. However, the influence of the mass transfer on the flow agreed well with experimental data.
5.2.2. Zhang et al. [79]
Zhang et al. [79] followed a procedure similar to that used by Darmana et al. [78], although in this case an E-E approach was used to simulate flow, mass transfer, and chemical reactions in square cross-sectional bubble column [19]. Zhang et al. studied physical and chemical absorption of CO2 bubbles in water and in an aqueous sodium hydroxide (NaOH) solution. They used a bubble number density equation for coupling of flow, mass transfer, and chemical reaction. The authors demonstrated the influence of the mass transfer and chemical reaction on the hydrodynamics, bubble size distribution, and gas holdup.
5.3. Mixing and Dispersion
5.3.1. Bai et al. [36]
Bai et al. [36] used LES with E-L approach to investigate the effect of the gas sparger and gas phase mixing in a square cross-sectional bubble column. The liquid phase was computed using LES, and a Lagrangian approach was used for the dispersed phase. They used the DBM and investigated the effect of two SGS models: Smagorinsky [21] and Vreman [80]. They compared the vertical liquid velocity and turbulent kinetic energy of the liquid phase at three different heights with PIV data and found that the model proposed by Vreman performed better than Smagorinsky model.
They further investigated the effect of the gas sparger properties (sparged area and its location) on the hydrodynamics in a bubble column and characterized the macromixing of the gas phase in the column in terms of an axial dispersion coefficient. They compared the predicted liquid phase dispersion coefficient with the literature correlations as shown in Figure 15. The range of superficial gas velocity investigated in work is low compared to what is common in industrial application. For large-scale reactors at high superficial velocities, Bai et al. recommended to extend the discrete bubble modelling with bubble coalescence and breakup.
5.4. Estimation of the Turbulent Dispersion Force
In the RANS approach, the drag and lift forces depend on the actual relative velocity between the phases, but the ensemble equations of motion for the liquid only provide information regarding the mean flow field. The random influence of the turbulent eddies is considered by modelling a turbulent dispersion force. By analogy with molecular movement, the force is set proportional to the local bubble concentration gradient (or void fraction), with a diffusion coefficient derived from the turbulent kinetic energy. The value of the turbulent dispersion coefficient is chosen to get an agreement with the measurement data and is not known as a priori.
In LES, the resolved part of the turbulent dispersion is implicitly computed, and hence one can use information from LES for calculating the magnitude of this force. The methodology depends on scales at which LES is to be applied. For instance, at the mesoscale, in the E-L approach, bubbles dispersed by drag and lift through turbulent eddies can be computed. At micro-scale LES, one might need to consider bubble coalescence and breakup phenomena along with a reasonable number of bubbles. It can be computationally expensive, but in view of increasing available computer power, this should become feasible soon.
5.5. Dynamics of Turbulent Structures and the Estimation of Design Parameters
The turbulent flows contain flow structures with a wide range of length and time scales which control the transport processes. The length scales of these structures can range from column dimensions (highest) to Kolmogorov scales (lowest). However, not all the scales of turbulence contribute equally to different transport rates and mixing. If only mixing is the important design criterion, then the knowledge about the mean flow pattern (large-scale structures) would generally suffice the purpose (Ekambara and Joshi [81]). However, for the prediction of the gas holdup, bubble size distribution, true mass transfer coefficient, and heat transfer coefficient, the knowledge about all the scales is important [82, 83]. Hence, it is imperative to identify the scales and dynamics of turbulent flow structures and their relationship with the rates of different transport process. The present empirical design practices do not consider these basic mechanisms and conceales the detailed local information about the relationship between the turbulence and the equipment performance.
The subject of quantification of local turbulent flow structures and reliable estimation of transport properties has been reviewed by Joshi et al. [84] and [82, 83]. The velocity and pressure data from LES were analyzed using the mathematical techniques such as multiresolution analysis [85], wavelet transforms (discrete and continuous), proper orthogonal decomposition (POD), and hybrid POD-wavelet techniques (Tabib and Joshi [86], Tabib et al. [87], Sathe et al. [88], and Mathpati et al. [89]). These techniques give the size, shape penetration depth, and energy content of all the flow structure in the system. This flow structure information can also be used for the construction of energy spectrum and for examining the scaling laws for turbulence in bubble columns. Such understanding of turbulence is expected to provide better insights into the transport phenomena. One such attempt has been reported by Deshpande et al. [90, 91].
5.6. Comparison of Turbulence Models
CFD provides detailed flow information within single- and multiphase reactors. Most popular and computationally inexpensive models such as k-ε model and Reynolds stress model (RSM) are widely used to predict the mean flow pattern. These models can give reliable estimation about the liquid phase mixing. However, they do not accurately predict the turbulence parameters such as turbulent kinetic energy and the dissipation rate due to inbuilt modelling assumptions as well as complexity of flow [11, 24]. These models are time averaged, and hence the information related to different turbulent structures is lost.
It is known that a large number of simplifying assumptions are made while deriving the k-ε and RSM models. Therefore, it is important to understand the gravity of these assumptions on the quantitative values of transport rates of k and ε due to convection, diffusion, and turbulent dispersion. It is also important to know the quantitative estimation of production and dissipation rates of k and ε. Therefore, it is important to estimate these five terms using k-ε, RSM, and LES models. From the LES simulations, the time series of velocity and pressure can be stored. These are subsequently used for the detailed comparison of k-ε, RSM, and LES models in terms of the rates of transport (convection, molecular and turbulent diffusion) and the rates of production, and the dissipation of k and ε for the case of dispersed bubbly flows [92].
6. Summary and Suggestions for Future Work
(1) E-E and E-L LES are promising approaches for predicting unsteady, buoyancy-driven flow inducing large-scale coherent structures for gas-liquid dispersed flow. Care should be taken to clearly identify the scales (micro, macro, or meso) at which LES should be applied, in order to decide the level of interface resolution and modelling required. The approach of LES at mesoscales (i.e., without explicitly tracking interface) using E-E and E-L description has been reviewed for gas-liquid dispersed flow.(2) Pioneering work of Milelli et al. [11, 24] has initiated the LES approach for gas-liquid dispersed flows. The main contribution comes from insights in the cutoff filter requirement and SGS modelling.(3) The simulation and the experimental measurement of Deen et al. [19] in a square cross-sectional bubble column have triggered a systematic development of the two-phase LES for both E-E and E-L approaches.(4) The concept behind the LES is very simple but characterized by a large number of choices (regarding numerical and physical modelling) that all have significant influence on the results. However, it offers great potential in terms of determination of statistical quantities and instantaneous information about flow structures. This information can be extremely useful for the prediction of other physical processes behaviour (e.g., transport of scalar (temperature, concentration), chemical reactions).(5) From LES simulation with E-E/E-L approaches that were reviewed in this work, it is recommended that:(a) The grid or filter size selection based on filter size to bubble diameter ratio Δ/ of 1.2 gives reasonable results. (b) The Smagorinsky constant, , is a modelling parameter rather than a physical constant. Although the constant value of the parameter gives satisfactory results, for unknown configuration, it should be estimated with Germano dynamic procedure (using the overall distribution of the constant through probability density).(c) The lift force is the main mechanism for the dispersion, and the lift coefficient should be estimated though sensitivity of interfacial forces on values of slip velocity and gas holdup. The lift coefficient in LES can be different from that in RANS.(d) The central difference scheme should be used for the discretization of advection terms for flow variables and high-order schemes (MUSCL, QUICK, or Second-Order) can be used for scalar variables.(e) The minimum time for gathering statistics should be at least one flow through time (as defined as ratio of the system height over the bulk (superficial) velocity). (6) In advent of computer hardware, the E-L approach appears very promising for the near future. Further work in mapping functions for two-way coupling can expedite the development of this approach that can be used as a means of both predicting the properties of specific turbulent flows and providing flow details that can be used like data to test and refine other turbulence-closure models.(7) The approach for BIT with extra production terms into the SGS-turbulent kinetic energy equation (following the procedure described by Pfleger and Becker [38]) has shown to be more effective than the approach involving a bubble-induced viscosity [37]. It can be that the enhanced eddy viscosity in LES does not represented as realistic physical model, as the SGS turbulent kinetic energy. Nonetheless, it is an interesting issue, and more work in investigating the BIT should be undertaken.(8) Treatment of the interphase forces needs more attention.(a) The drag and nondrag forces (lift, virtual mass force) can be modelled using resolved field approaches. The modelling of these forces for the SGS and their effect on the overall simulation results need to be evaluated. (b) One finds strong dependency of the bubble dispersion on the value of transverse lift force coefficient. The transverse lift, which depends on the bubble-liquid relative velocity, seems to be the main mechanism responsible for the spreading of the bubbles. It will help if one can estimate the separate contributions of each of these forces. (c) The virtual mass force has little influence on simulation results. So far, a constant coefficient has been used in all the investigations; however, dependence on void fraction has been shown in experiments. It would be good to have a correct description in order to improve results near the inlet where bubble acceleration effects are important.(9) The strong coupling between subgrid-scale (SGS) modelling and the truncation error of the numerical discretization can be exploited by developing discretization methods where the truncation error itself functions as an implicit SGS model. Such attempt can be useful and go in the direction of finding a universal SGS model.(10) In order to use LES for reliable predictions at minimum computational costs, understanding of the influence of discretization methods, boundary conditions, wall models, and numerical parameters (e.g., convergence criterion, time steps, etc.) is essential. The contribution focusing on these aspects should be undertaken for both E-E/E-L approaches.(11) Substantial development has been achieved in LES in the last decade for understanding bubbly gas-liquid dispersed flow. However, it is mainly restricted to low superficial gas velocities and gas fractions. Future work should focus on industrially relevant large-scale reactors at high superficial gas velocity. The modelling of bubble coalescence and breakup might be necessary, along with further clarity in filtering operations.(12) Joshi and coworkers have used LES for the identification of flow structures and their dynamics. They have proposed a procedure to use this information for the estimation of design parameters. Substantial additional work is needed for finding 3D information on the structure characteristics such as size, shape, velocity, and energy distributions
Nomenclature
| : | Drag force coefficient |
| : | Lift force coefficient |
| : | Virtual mass coefficient |
| : | Turbulent dispersion coefficient |
| : | Model parameter in turbulent dissipation energy equation (=1.44) |
| : | Model parameter in turbulent dissipation energy equation (=1.92) |
| : | Smagorinsky model constant |
| : | Model constant (Sato and Seguchi [37]), (=0.6) |
| : | Mean bubble diameter |
| : | Diameter of the column (m) |
| : | Force originating due to pressure (N/) |
| : | Gravitational force per unit volume of dispersion (N/) |
| : | Lift forceper unit volume of dispersion (N/) |
| : | Virtual mass force per unit volume of dispersion (N/) |
| : | Turbulent dispersion force per unit volume of dispersion (N/) |
| : | Wall lubrication force per unit volume of dispersion (N/) |
| : | Wall deformation force per unit volume of dispersion (N/) |
| : | Height of column (m) |
| : | Turbulent kinetic energy per unit mass () |
| : | Subgrid-scale turbulent kinetic energy () |
| Characteristic of strain tensor of filtered velocity | |
| Instantaneous axial velocity (m/s) | |
| : | Width of column (m). |
| Δ: | Filter width |
| : | Simulation time step |
| : | Grid-scale component of scalar |
| : | Resolved component of scalar |
| : | Filtered component of scalar |
| : | Turbulent energy dissipation rate per unit mass () |
| : | Fractional phase holdup |
| : | Fractional gas phase holdup |
| : | Density (kg/) |
| : | Density of liquid (kg/) |
| : | Effective viscosity of phase (Pa s) |
| : | Effective viscosity of liquid phase (Pa s) |
| : | Molecular viscosity of liquid phase (Pa s) |
| : | Bubble-induced viscosity (Pa s). |
| : | Phase, = gas phase, = liquid phase |
| BI: | Bubble-induced. |
| BIT: | Bubble-induced turbulence |
| E: | Eulerian |
| L: | Lagrangian |
| SGS: | Subgrid-scale. |
Acknowledgment
N. G. Deen would like to thank the European Research Council for its financial support, under its Starting Investigator Grant scheme, contract number 259521 (cutting bubbles).