Abstract
In this work, we report on Euler-Euler large eddy simulation (EELES) of dispersed bubbly flow in a square cross-sectional bubble column. Simulations are performed using the Neptune_CFD package, and results processed using the SALOME platform. The motivation to undertake this study is to check our implementation of the Smagorinsky subgrid-scale (SGS) model into NEPTUNE_CFD. We outline all the physical models used, and we present instantaneous realizations of velocity and void fraction fields in order to illustrate the structure of the turbulence field, and long-time averaged results, to compare with analogous simulations performed using the CFX-4 code and experimental data. The same physical models and constants have been used in both the CFX-4 and NEPTUNE_CFD codes, except the SGS model, which is Smagorinsky in case of NEPTUNE_CFD and a one-equation model in CFX-4. The results obtained with EELES compare reasonably well with experiment, meaning in particular that the implementations have been successful. Some perspectives on the further use of EELES are also given.
1. Introduction
The aim of the NURESIM [1] integrated project is to provide the initial step towards a common European software platform for analysing, modelling, and exchanging data for nuclear reactor simulations. One of the main goals of the project is the integration of state-of-the-art physical models in a common, open-software platform, including the latest advances in reactor core physics, thermal-hydraulics, and coupled multiphysics modelling capabilities. The objective of the thermal-hydraulic subproject within NURESIM is the improvement of the predictive capabilities of the simulation tools for key two-phase-flow, thermal-hydraulic processes that can occur in nuclear reactors, focusing on two high-priority issues: critical heat flux (CHF) [2] and pressurized thermal shock (PTS) [3]. As a part of this activity, the Paul Scherrer Institute is investigating the potential of the EELES approach for simulating bubbly flows, which are important for understanding CHF and PTS, using the Neptune_CFD code [4], the CFD component of the NURESIM platform. In view of this, it was thought desirable to carry out EELES in a three-dimensional, square-cross-sectional bubble column, based on the experiments of Deen et al. [5], to validate the implementation of the SGS model in the code. In the present work, for reasons of simplicity, the performance of Neptune_CFD with the EELES paradigm is assessed based on the Smagorinsky SGS model. These results, and those obtained with the CFX-4 code, are compared against experimental findings [5]. We point to the potential weaknesses of the current approach, and propose directions for future research and development in the area.
Previous work on EELES (often referred to as two-fluid LES) for bubbly flows include Milelli et al. [6], who used Smagorinsky and dynamic SGS models to study the motion of a bubble plume in a cylindrical test section 50 cm in diameter, 40 cm high, for a maximum superficial gas velocity of 6.11 mm/s. A similar approach was taken by Deen et al. [7] to simulate a bubble plume in a test section with a square cross-section of 15 cm 15 cm, a height of 45 cm, and superficial gas velocity of 4.9 mm/s. The authors compare the EELES prediction with those obtained using a model and conclude that EELES results are closer to experimental measurements. Turbulent bubbly shear flows have also been analyzed [8] for a square channel 30 cm wide, 4 cm deep, and 60 cm high. A splitter plane was used to divide the inlet channel in two equal parts. A mixture of water and air, with different velocities and void fractions, was introduced at the inlet sections. More recently, an LES calculation incorporating the Smagorinsky SGS model has been reported for a rectangular column [9]. LES predictions were compared with transient results obtained using a mixing length model and the model; all three models produced similar results.
1.1. The Deen Bubble Column Experiment
In the Deen series of experiments [5], the hydrodynamics of the bubble plume is studied in the same apparatus, but for different inlet conditions. In the test to be investigated here, bubbles are created by injecting air with a superficial gas velocity of 4.9 mm/s through a central sparger at the bottom of the vessel of dimension 5 cm 5 cm. The experiment operates at atmospheric pressure. A sketch of the Deen experiment is given in Figure 1.
2. Physical Models
2.1. Euler-Euler Approach
The governing equation for conservation of mass, with no mass exchange between the phases, takes the form: in which is time, is the volume fraction, is the phase density, the phase velocity, and is the phase indicator. Conservation of momentum is governed by the relation: where is shear stress, is pressure, and is the gravity vector. The last term on the right-hand side represents the momentum exchange between the phases due to interface forces and requires modelling. The models we have used are given in Section 2.3.
2.2. LES of Turbulence
Turbulence is modelled in the liquid phase of the flow. The first step in the derivation of the governing equations for LES is the decomposition of the velocity field into resolved (grid scale) and unresolved (SGS) parts. For each component, we write where the resolved velocity on the left side is equal to the difference between the true instantaneous velocity (first term on the right) and the unresolved velocity (second term on the right). The filtering of the nonlinear advection term on the left-hand side of (2) leads to an additional stress-like term, which has to be modelled. Hence, we write The first term represents the laminar stress, with the strain rate tensor defined as and the second term derives from the turbulence effects and needs closure. In this work, we have used two closure models: Smagorinsky model (in Neptune_CFD), and a one-equation model (in CFX-4).
2.3. Smagorinsky Model
To complete (4), we use a formulation [10] for the turbulent stress:The second term on the left-hand side, the trace of the SGS stress tensor, is implicitly added to the pressure, while the deviatoric part is modelled by the expression on the right-hand side, that is, by the Boussinesq eddy viscosity () concept. The Smagorinsky model is an algebraic one, with the eddy viscosity calculated from: In (7), is a filter width, associated with the cell size. In the present work, it was estimated as the cube root of the cell volume: although other definitions are possible. is the Smagorinsky constant, and was set here to 0.12 since the flow is dominated by the large-scale structures generated by the buoyancy force, as pointed out in [11]. The eddy viscosity, defined by (7), is used in the stress tensor appearing in the momentum (4), giving
2.4. One-Equation SGS Model
A disadvantage of the Smagorinsky model is the loss of information resulting from use of the deviatoric part of the SGS stress tensor only. In CFX-4, we have implemented a transport equation for SGS kinetic energy , defined by of the following transport form [12]: where is the turbulence production term defined as: Using this model, the eddy viscosity is then calculated from: with model constants and . A detailed description of the one-equation SGS model for dispersed bubbly flows can be found in [12].
2.5. Interfacial Forces
To close the momentum equation set for the two phases, the various interphase exchange terms have to be modelled. The interphase exchange terms relevant for the configuration considered in this work are drag, buoyancy, virtual mass, and lift forces. In the framework of the EELES presented in this work, interfacial forces are applied only to the resolved field, while their SGS counterparts are neglected.
The drag force is modelled from the resolved velocity field as follows [13]: where is the drag coefficient, defined here in terms of the Eötvös number (): In the present application, a bubble size of 4.0 mm is specified (the value reported in [5]), and the coefficient of surface tension N/m, giving and .
For the virtual mass force, we use the model of Drew and Lahey Jr. [14]: where is the virtual mass force coefficient, here set to 0.5.
The lift force accounts for the transverse migration of bubbles under the influence of the liquid shear field. Following earlier work [14, 15], it is modelled here according to in which the lift coefficient, , is set to 0.5.
3. Requirements for the Numerical Grid
In the present work, we couple the Euler-Euler approach for multiphase flow with LES, and therefore have to consider the resolution requirements of both techniques simultaneously in order to choose a satisfactory grid. A basic requirement is that the control volume size should be large enough to encompass all the interface details. This is the intrinsic assumption in the derivation of the Euler-Euler model equations, and strictly has to be satisfied at the discrete level as well.
In LES, the SGS model is often very simple, and only drains energy from the resolved field without feed-back. Therefore, our goal in LES is to resolve as much of the flow field as possible, and to have as fine a grid as feasible for the available computer hardware. Since the Euler-Euler approach specifies the minimum control volume size, whereas for LES we are invariantly seeking as fine a grid as possible, the requirements for the numerical grid may sometimes be in conflict. The point is illustrated in Figures 2 and 3, which show the turbulence spectra under different modelling conditions. For a successful LES, we must have a filter width (Δ) in the inertial subrange region, and all scales of motion larger than that (left of Δ in Figures 2 and 3), must be accurately resolved on the numerical grid. If, however, we have a bubble diameter larger than Δ (see Figure 2), it is obvious that they would induce some large-scale motions which are not properly accounted for by the EELES, since there is no information on the interface details or their influence on the resolved large-scale motions. If in addition we use a model for bubble-induced turbulence, it would drain the energy from the resolved field, further deteriorating the accuracy of the resolved field. This is illustrated by the saw-like shading in Figure 2. This influence on the resolved part of the spectra is not acceptable for LES. This is not just a conceptual consideration, as discussed in [6], and it is explained in more detail in Section 3.1 . The situation which is safe for EELES is shown in Figure 3, where is smaller than Δ, and all bubble-induced scales which cannot be calculated using EELES approach fall into the SGS part of the spectra.
3.1. The Milelli Condition
The grid size considerations discussed above, are not new. In the present work, we have indicated why the bubbles must be smaller than the cell size from the point of view of the energy spectra and modelling closures for the interfacial forces. A systematic a posteriori analysis of the minimum ratio of the bubble and cell sizes for LES modelling of free bubble plumes is reported by Milelli et al. [6], which resulted in the following criterion: Equation (18) states that the cell size must be at least 50% larger than the bubble diameter for accurate LES. The situation is illustrated in Figure 4. In the Deen experiment, the mean bubble size is 4 mm. Since the flow is dominated by the energetic, large-scale structures in the core of the flow, with wall effects having a smaller impact on the overall flow field, we have created a uniform grid with 30 30 100 cubic cells. This results in a cell size of 5 mm, hence, which is below the optimum value. A coarser mesh which satisfies the Milelli condition has also been used, but with no significant change in the computed results [12].
It might be argued at this point that the grid used in the present simulations is too coarse for LES, at least from the point of view of capturing all the relevant (large) scales. Judging from the flow patterns that can be expected in a square column such as this, with the bubble plume meandering from one direction to the other, it is quite conceivable that the largest, and therefore the most energetic, eddies will be of the size of the domain cross-section. Therefore, we are confident that the grid resolution we employ (30 30 cells in the cross-section) is sufficient to resolve these large structures. Should the topology of the flow be different, and should the flow be dominated by small, near-wall scales of motion, this resolution would not be adequate. The adequacy of the grid resolution used is proved in the following section.
4. Simulation Details
In this section, we briefly summarize all the relevant simulation details.
Inlet Superficial Velocities:
for liquid: = 0.0 m/s;
for gas: = 0.0049 m/s.
Spatial Discretization:
number of cells in computational domain: 90 000;
advection scheme: bounded central scheme for Neptune_CFD;
quick scheme for
CFX-4.
Temporal
Discretization:
the linear backward time differencing for Neptune_CFD;
quadratic backward differencing for CFX-4;
CFL = 1.0, variable time step,
approximately equal to 0.05 second for Neptune_CFD constant time step of 0.05 second for CFX-4, leading to CFL 1.0;
integral time of simulation: 400 seconds;
simple time-average is used in both, not
weighted by volume fraction.
5. Results
In this section, the results for our EELES with Neptune_CFD are presented. We start with the instantaneous velocity and void fraction fields to show the flow patterns in the bubble column flow, and continue with time-averaged velocity and fluctuating liquid velocity profiles. Time-averaged profiles obtained with Neptune_CFD and CFX-4 are compared with the experimental measurements of Deen et al. [5].
5.1. Instantaneous Results
Instantaneous velocity fields are shown in Figure 5. Velocity vectors are coloured with the velocity magnitude. Two things are immediately apparent. First, the velocity field does not show any signs of the growing instabilities to which LES computations are very often prone. This is not surprising, since the advection scheme is bounded. Another thing worth noting is that the velocity field pattern is changing dramatically in time. Bubbles are moving up through the large turbulent structures, driven by their buoyancy, resulting in different bubble paths at each time instant. This behaviour was also reported for bubble plume experiments [16]. Figure 6 shows two isosurfaces of constant void fraction. The blue one corresponds to the low-void fraction of and the red to . The threshold for the blue isosurface gives a good illustration of the bubble plume shape, whereas the red one indicates the position of the free surface.
(a)
(b)
(c)
(a)
(b)
(c)
5.2. Time-Averaged Results
In this section, we compare the simulated profiles of liquid and gas vertical velocities and liquid turbulence intensities against the reported experimental data [5]. The measurements were taken along the centreline of the horizontal plane at height 25 cm.
Figure 7 shows comparisons of vertical liquid and gas vertical velocity profiles. The liquid velocity profile is somewhat under-predicted by both codes. In contrast, the gas velocity prediction compares very well with experiment in both cases. From these time-averaged results, it might be concluded that the drag force is underestimated in both codes, probably due to a wrong assumption for the drag coefficient , or for the relative velocity between the two phases [16].
(a)
(b)
The vertical component of the resolved vertical fluctuating liquid velocity is plotted in Figure 8(a). The qualitative comparison is encouraging for both codes, since the twin-peaked shape seen in the experiments is reproduced. Quantitatively, however, both codes under-predict the magnitude to some extent. Figure 8(b)shows the resolved lateral fluctuating liquid velocity, which is predicted very well. The total turbulent kinetic energy, plotted in Figure 9, in spite of the under-prediction of vertical fluctuating liquid component, is well predicted by both codes in the middle of the column and somewhat over-predicted close to the walls. This over-estimation of the near-wall kinetic energy might be attributed to insufficient resolution for the near-wall structures.
(a)
(b)
6. Conclusions and Perspectives
Implementation of the Smagorinsky SGS model in the Neptune_CFD code has been validated using data from the Deen bubble column experiment. Results have also been compared with predictions obtained from simulations performed using the code CFX-4. Generally, liquid velocity profiles are under-predicted compared with measured data, whereas the gas velocity profiles are predicted very well. This can only be attributed to a wrong assumption for the drag coefficient. The liquid velocity fluctuations predicted by both codes compare well with experiments, except for some under-prediction of the vertical component. Overall, the results obtained with the two codes are consistent, regardless of the fact that different SGS models have been used. This is a clear indication that the SGS model influences results only to a very small extent.
When envisaging the potential of LES for modelling multiphase flows in general, one should clearly distinguish between the scales at which LES might be used, since this implies the level of detail of accurate interface resolution to the degree of modelling that can be tolerated. Simulations at mesoscales imply an Euler-Euler description of the interface between the phases, such as the one described in this work. However, if LES for multiphase flows is applied at microscales, explicit interface tracking procedures would be needed. This has already been proposed in [17], and is generally referred to as large-scale simulation (LSS).
The principal advantage of EELES over LSS is that since the interface details are not calculated explicitly, the simulations may be carried out at lower cost. The principal disadvantage of EELES is that the most influential interfacial forces (lift and drag) are modelled for the large-scale field, meaning that the question of how to model these forces remains as open for EELES as it is with RANS. LSS, on the other hand, explicitly resolves the large-scale part of the interfacial forces, leaving the modelling at the SGS level, where the effects are smaller and hence less influential on the accuracy of the results. LSS modelling comes at a heavy price, due to the need to resolve all relevant interface details, imposing huge demands on computing power. Applying LSS to industrial-scale problems is beyond the current state of computing resources.
Another disadvantage of EELES stems from the different resolution requirements imposed by the Euler-Euler description of the two fluids and that of the LES approach itself, expressed in [6]. Since for an Euler-Euler description, the minimum cell size must be larger than the interface detail, accurate LES requires as fine a grid as possible. This requirement is particularly stringent in the near-wall regions, where the turbulent structures (e.g., streaks) are very small, and a very fine grid is needed to capture all relevant details. The disadvantages of EELES just outlined are not so pronounced for flows in which the bubbles are small compared to the grid size, such as for the bubble column studied in this work. Generally, one should use EELES if the flow is likely to be unsteady, with large coherent structures generated either by buoyancy or obstacles in the flow. EELES, as LES itself, should be avoided for wall-dominated flows. In other words, large eddy simulation should be used in the cases in which large eddies are present. Flows with density stratification in large channels (relevant to the PTS issue), in addition to bubble column devices, are examples where EELES has the potential to give accurate insight into the phenomena taking place and the ability to quantify them.