Abstract
The development of a turbulent mixing layer at mixing of two horizontal water streams with slightly different densities is studied by the means of numerical simulation. The mixing of such flows can be modelled as the flow of two components, where the concentration of one component in the mixing region is described as a passive scalar. The velocity field remains common over the entire computational domain, where the density and viscosity difference due to the concentration mainly affects the turbulent fluctuations in the mixing region. The numerical simulations are performed with the open source code OpenFOAM using two different approaches for turbulence modelling, Reynolds Averaged Navier Stokes equations (RANS) and Large Eddy Simulation (LES). The simulation results are discussed and compared with the benchmark experiment obtained within the frame of OECD/NEA benchmark test. A good agreement with experimental results is obtained in the case of the single liquid experiment. A high discrepancy between the simulated and the experimental velocity fluctuations in the case of mixing of the flows with the slightly different densities and viscosities triggered a systematic investigation of the modelling approaches that helped us to find out and interpret the main reasons for the disagreement.
1. Introduction
Interaction of liquid streams with different densities and temperatures is of particular interest in nuclear industry. Notable examples are turbulent mixing in junctions of the primary coolant piping or boron mixing in the reactor core of the pressurized water reactor. The turbulent wake mixing zone, in which two streams interact, depends on the physical properties of the two liquids involved. Computational Fluid Dynamics (CFD) simulations are commonly used to predict the turbulent phenomena and may successfully complement the experimental data. The thoroughly validated CFD models are particularly useful for the investigation of key turbulent mechanisms. Similar to the experiments, also the numerical simulations can be affected by uncertainties, in large part due to the increasing complexity of the physical models that are very sensitive to the uncertain boundary conditions [1]. The CFD models need to be validated on the detailed experimental data, but it can be also the other way around; the detailed numerical simulations can be of great support to verify the regularity of experimental methods in order to avoid systematic experimental errors [2].
The focus of this paper is the simulation of the turbulent mixing wake at mixing of two water streams with a small difference in densities. The experimental data of the GEMIX (Generic Mixing Experiment) experimental test facility located at the Paul Scherrer Institute (PSI) are used to validate the accuracy of the numerical simulations. The measured data are obtained in the frame of participation within the OECD/NEA benchmark exercise [3]. The numerical simulations are performed with the open source code OpenFOAM [4] using two different approaches for turbulence modelling, the Reynolds Averaged Navier Stokes equations (RANS) and the Large Eddy Simulation (LES).
2. Mixing Experiment
The confined wake flow water mixing experiments have been carried out at the GEMIX facility [3]. The Y-shaped squared flow channel (on Figure 1(a)) is used for the experiments. The inlet section consists of the upper and the lower leg. The upper and the lower streams have the same volumetric flow rate and the same shape of the inlet velocity profile. The two inlet legs, separated by 3° splitter plate, lead to the mixing section of the channel. The mixing section with the square cross-sectional area of 50 mm by 50 mm is 600 mm long (see Figure 1).
(a)
(b)
Two experimental cases are considered that differ in the properties of the flow entering into the upper and the lower leg. Hereinafter they are named as the “single liquid” and the “two-liquid case.” In the “single liquid” case, a tap water at 20°C is used for both inlet legs. In the “two-liquid” case, the tap water at 20°C is flowing through the upper inlet leg and the deionized water at 22.5°C with added sucrose is flowing through the lower leg. The mass fraction of the sucrose in conjunction with a temperature adjustment is used to increase the density of the water flowing through the lower leg by 1% while the kinematic viscosity of both streams differs by approximately 2%. The physical properties of the liquids used in the experimental cases are listed in Table 1.
In the mixing section a particle image velocimetry (PIV) and a laser induced fluorescence (LIF) are applied to measure the mixing process. One wire mesh sensor at a time is installed to measure the cross-sectional concentration profile. The five measuring planes located 50, 150, 250, 350, and 450 mm downstream the splitting plate are shown in Figure 1(b). The measurement details are provided in [3].
During the experiment the velocity values in and direction were measured and the RMS velocity fluctuations were calculated. The experimentalists assumed that the velocity fluctuations are the same in the and the direction and consequently they did not measure the velocity fluctuations in the direction. The turbulence kinetic energy is calculated from the velocity fluctuations measured in the experiment aswhere and are the velocity fluctuations in the and the direction, respectively.
3. Simulation Model
Numerical simulations are performed with the OpenFOAM CFD code, version 1606+ [4]. The TwoLiquidMixingFoam solver is designed for modelling two miscible incompressible fluids mixing and is used to simulate both experimental cases. This solver solves one momentum equation, equation for turbulence transport, and one additional equation for the transport of the passive scalar (concentration ):where is the velocity vector, is the molecular diffusion coefficient, is the turbulent (eddy) viscosity, and is the turbulent Schmidt number. From the solved concentration field, a corrected density field is obtained as
The new density is then used in momentum and turbulence transport equations. The results include single velocity, pressure, and concentration fields.
The two experimental cases described in Section 2 are being simulated. In further text they are described as the “Single liquid” and the “Two-liquid” case. The diffusion coefficient and the turbulent Schmidt number were set to 2.0 × 10−5 m2/s and 1.25, respectively, for both cases.
Two different approaches are used to model the turbulence. In the first approach, the Reynolds Averaged Navier Stokes (RANS) equations are used to resolve the flow [5]. Within the RANS approach the entire spectrum of turbulent fluctuations is modelled using the - SST (Shear Stress Turbulence) model with the default coefficients and the scalable wall functions [6]. The second approach, the Large Eddy Simulation (LES), applies space filtering on transient momentum equations and uses the subgrid scale Wall Adapting Local Eddy (SGS WALE) viscosity model for proper scaling in the near-wall region [7]. The mesh for the considered LES calculation is sufficiently dense that the later contribution is small enough and can be neglected. The turbulence kinetic energy for the LES calculation is thus calculated from the resolved velocity calculations in the same way as in the experiment (1).
3.1. Computational Mesh and Boundary Conditions
Computational domain for the RANS calculations includes 50 mm long inlet legs before the junction and 450 mm long mixing section. The origin of the coordinate system in the domain is located at the tip of the splitter plate (see Figure 1(a)). The coordinates , , and represent the streamwise, the crosswise, and the spanwise direction, respectively. The mesh used for the RANS calculations is shown in Figure 2. The mesh cells are uniformly distributed in the streamwise direction, while they are refined in the wall-normal direction near the channel walls and in the middle of the mixing section (crosswise and spanwise) using a linear Stretching Ratio (SR). The convergence of the results during mesh refinement is studied on four different meshes, listed in Table 2. The finest mesh A and the coarsest mesh D consisted of 1,500,000 and 60,350 elements, respectively. The wall-resolved RANS using a second order method usually requires for the mesh near the wall. This requirement is fulfilled for the meshes A and B (see Table 2) and shows up when comparing the turbulence kinetic energy profiles just after the splitter plate (at = 0.05 m) in Figure 3. The mesh B shows no major discrepancies comparing to the results with the finest mesh A; hence it is used for all further RANS simulations in this study.
The computational domain used for the LES calculations includes the 100 mm long inlet legs before the junction and the 450 mm long mixing section. The subgrid scale model in the LES simulation requires denser mesh with the mesh cells having much lower aspect ratio (AR) than in the case of RANS. The mesh parameters for the LES are listed in Table 3. The mesh cells are uniformly distributed in the streamwise direction and refined in the wall-normal direction near the channel walls and in the mixing section (crosswise and spanwise) to satisfy the condition at the walls. The mesh parameters are comparable to the ones used in the wall-resolved LES simulations of the channel flow [2].
Velocity and turbulence kinetic energy boundary conditions at the inlet for the RANS calculations were extracted from the experimental profiles provided within the benchmark with the average inlet velocity value of 0.6 m/s. As already mentioned, the LES domain uses longer inlet legs in order to achieve a fully developed flow at the inlet. The mapped boundary condition implemented in the OpenFOAM uses the recycling method to generate the inflow turbulence from the initial disturbances in the velocity field. For both types of calculations the no-slip boundary condition is applied on the walls and the Neumann boundary condition (gradient of the variable is equal to zero) is used at the outlet.
Using a Flow Though Time (FTT) measure, which is the ratio of the streamwise length of the domain to the mean velocity magnitude, the time averaging in the case of LES calculations started after 6 FTT, at statistically steady state flow conditions, in order to nullify the effect of initial transient events at the inflow. To assure the convergence of the results approximately 11 FTT (10 s of the flow time) are simulated. Residuals are set to 10−7 for all quantities.
4. Results and Discussion
4.1. On the Inlet Flow Conditions
In the case of RANS the experimental profiles for the inlet velocity and the turbulence kinetic energy are directly imposed at the inlet boundaries of both inlet legs. This is however not possible in the case of LES approach. Instead, the fully developed flow at both inlet legs can be reproduced by using the recycling method for generation of the inflow turbulence. After some appropriate length () from the inlet, the velocity field was recycled back to the inlet. The recycling distance must be long enough, so that the large flow structures are not affected. In our case the used recycling length was 5 cm. The adequate length in the streamwise direction was verified by calculating the time averaged autocorrelation function of the instantaneous velocity components along the central line of the inlet leg (Figure 4) using
One of the key parameters for the wall-bounded flows is the mean friction velocity (), which can be calculated from the simulation results using the local wall-shear stress (near-wall velocity gradient), and the friction Reynolds number ():where is the hydraulic diameter of the inlet leg that amounts to m. The mean friction velocity calculated from the results is m/s and the friction Reynolds number is . The bulk Reynolds number based on the mean streamwise velocity and is = 20,000.
The velocity profiles and the fluctuations in the channel flows are usually compared with the results from the similar studies in the nondimensional form using the normalization with the friction velocity : where and represent the nondimensional wall distance and the nondimensional mean velocity, respectively. Mikuž and Tiselj [2] used LES approach to study similar channel flow with friction Reynolds number . The obtained friction velocity was m/s. Zhang et al. [8] used DNS approach, to study the flow in a duct of square cross section with friction Reynolds number . Figure 5 shows the comparison of our LES streamwise velocity profiles (located at the inlet leg 5 cm before the splitter plate), with the LES velocity profiles obtained by Mikuž and Tiselj [2] and DNS results of Zhang et al. [8]. The profiles are similar near the wall; some smaller discrepancies can be observed at above 100, where the velocity profiles in duct are slightly steeper.
Figure 6 compares our LES simulations of the three normalized components of the inlet velocity fluctuations with the LES velocity fluctuations obtained in [2] and with DNS fluctuations by Zhang et al. [8]. As discussed in [9], the resolved velocity fluctuations obtained by the LES should be smaller than those calculated by the DNS or the measured ones. This is shown also in Figure 6, at least for the fluctuations in the crosswise and the spanwise direction. The streamwise fluctuations for the both LES calculations are however somewhat higher than for the DNS results. In general, the discrepancies between the 3 cases in Figure 6 are relatively small. These results show that the LES inlet flow conditions are reasonably well predicted and approximately correspond to the fully developed inlet flow.
(a)
(b)
(c)
The comparison of the inlet flow conditions between the experiment and those obtained by the RANS and the LES simulations is shown in Figure 7. It can be seen that the streamwise velocity profile in the experiment is more flat (due to the installed flow straightener grids) as in the case of the LES simulation. The turbulence kinetic energy near the wall is much higher in the case of LES. It should be noted that the experimental inlet boundary conditions cannot be simply reproduced in the case of LES simulations, where the fully developed flow conditions were used instead. As shown in Figure 7, the experimental profiles of inlet velocity and can be directly imposed as inlet boundary in the case of RANS.
(a)
(b)
4.2. Single Liquid Case: Mixing Region
Figure 8 shows the streamwise velocity profiles in the mixing region of the channel at different measuring positions. It can be seen that the streamwise velocity profiles are relatively well predicted by the RANS and the LES approach.
The turbulence kinetic energy profiles in the mixing section are shown on Figure 9. The profiles seem to be well predicted by both simulations at all measuring positions except for the last one, located 0.45 m away from the splitter plate. At locations = 0.05 m and = 0.015 m the simulated profile predicts a slight local drop near , which is not observed in the experiment. This seems to be the consequence of the eddy viscosity modelling approach that predicts lower at locally lower velocity gradients. If we compare the streamwise velocity profiles on Figure 8 and profiles on Figure 9, it can be seen that, in the narrow region around = 0, the velocity profile is relatively flat with a small local velocity gradient that leads to the local minimum in the profile. Similar results, with the deficit in the profiles, were observed also in the measured data of the previous benchmark experiment performed at GEMIX facility [10]. The RANS results show somewhat better matching with the experimental results, whereas the comparison of the LES results shows some discrepancy, especially in the region near the walls. This seems to be the consequence of different inlet profiles of the velocity and in the case of LES calculation. The higher from the inlet propagates into the mixing zone and with the increasing distance from the splitter plate slowly approaches to the experimental values.
For the LES case the most credible is the comparison of the directly calculated velocity fluctuations between the simulations and the experimental results. Figure 10 shows the comparison of the streamwise and the crosswise velocity fluctuations in the mixing zone at several locations along the channel. Crosswise components of the velocity fluctuations are reasonably well predicted by the simulations, except in the last measuring position. On the other hand, streamwise components are much better predicted away from the splitter plate, while significant overprediction just after the flows merging can be observed, especially in the proximity of the walls. This further shows that the inlet flow condition can be the main reason for discrepancy.
4.3. Two-Liquid Case: Mixing Region
Since the Reynolds number, the density and the kinematic viscosity are very similar for the single liquid and the two-liquid cases (1% difference in the density and 2% difference in the viscosity), the experimentalists did not measure the inlet conditions for the mixture of water and sucrose assuming that the velocity profiles and fluctuations should be also similar.
The comparison of the measured and the calculated results (RANS and LES) between the single liquid and the two-liquid case at the location just behind the splitter plate are shown in Figure 11. The calculated profiles of the velocity and the turbulence kinetic energy are indeed rather similar for both cases. The experimental velocity profiles also match for both cases, but a huge difference can be observed for the experimental turbulence kinetic energy profiles in the mixing region, where the turbulence kinetic energy for the two-liquid case is more than 10 times higher than the one from the single liquid case.
(a)
(b)
In Figures 12, 13, 14, and 15, the calculated profiles of the two-liquid case are compared with the experiments at different locations along the mixing region. The calculated velocity profiles in the mixing region for the two-liquid case match the profiles from the experiment as can be seen on Figure 12. Figure 13 shows the comparison of concentration profiles, where the simulated profiles are somewhat wider than the experimental ones. However, the calculated turbulence kinetic energy profiles strongly underpredict the experimental data, as shown on Figure 14. The underpredicted turbulence kinetic energy is the main reason why we continued to use the LES approach also in the two-liquid case. Namely, the LES approach resolves the velocity fluctuations and offers better insight into the behaviour.
The measured and the calculated velocity fluctuations for the two-liquid case are shown in Figure 15. The measured velocity fluctuations in the centre of the channel are approximately the same in both directions, as are the velocity fluctuations calculated with the LES. The fluctuations are well predicted for the single liquid case (Figure 10) but are greatly underpredicted for the two-liquid case (Figure 15). Despite the small difference in the density and in the kinematic viscosity in the two-liquid case, only 1% and 2%, respectively, the measured turbulence kinetic energy is by the order of magnitude higher than for the single liquid case. Several possible reasons for such turbulence increase are addressed and discussed in the sensitivity analysis section.
4.3.1. Sensitivity Analyses
Trying to understand the discrepancies, the sensitivity calculations with different diffusion coefficient and turbulent Schmidt number values were also performed. These two parameters influence only the concentration profile, while the other flow quantities remain unaffected. Figure 16 shows the influence of the turbulent Schmidt number on the shape of the mixing layer. It can be observed that the smaller the turbulent Schmidt number, the wider the mixing zone.
Despite the different inlet velocity profile shapes (uniform, fully developed) used in the calculations, the velocity fluctuations do not increase. To investigate the possible influence of the turbulence closure model, alternative turbulent models (standard - and Launder-Reece-Rodi - LRR) were also tested within the RANS approach, but no improvements were obtained, as can be seen in Figure 17. The same conclusion can be drawn when different solver methods within the OpenFOAM code were tested. The turbulence kinetic energy was still by the order of magnitude lower than in the experiment.
4.3.2. Buoyancy Effect on the Production of the Turbulence Kinetic Energy
The additional turbulence kinetic energy production term due to the buoyancy [11] can be written aswhere , and are eddy viscosity, numerical constant, density, and gravitational acceleration, respectively. The implementation of the buoyancy production term in the turbulent equation within the OpenFOAM is not possible without making significant changes to the source code. Therefore, the buoyancy effect was estimated by superposition of the buoyancy source term on the values. Based on the current simulation results, the increase of due to the buoyancy was calculated from (7). The increase of in the mixing zone (around m) amounts to 18% at most, which is much too small to explain the underprediction of the measured values. It should be noted that here the buoyancy effect is simply superimposed on the values calculated without this term. In the real simulation the feedback effect on the momentum equation appears. Namely, the increased turbulence kinetic energy due to the buoyancy would increase the turbulent viscosity, which would dampen the velocity gradients in the momentum equation. The lower velocity gradients would in the next step enter into the turbulent equations resulting in the lower turbulent kinetic energy values. Hence, the estimated increase of due to the buoyancy can be assumed as the highest possible.
4.3.3. Refractive Index Effect in the Experiment
The above sensitivity analyses offer no plausible explanation that the high disagreement between the simulated and measured profiles in the mixing region may arise from the modelling effects. Another possible source of error may arise also from the measured data, which is briefly discussed here.
When the two miscible liquids are mixed high concentration gradients may occur. According to Heidcamp [11] the refractive index of the two liquids involved in this benchmark differs by 4 × 10−3. When the two liquids with the different refractive indexes are mixed, the refractive index gradient appears, which affects the light propagation through the liquid. Concentration and the refractive index fluctuations may cause fluctuations in the dissipated light propagation and this could have the effect on the overestimated values of the velocity fluctuations.
Particle image velocimetry camera sees the particle at one position, which is corrupted due to the refractive index gradient. On the next measurement, the camera sees the same particle at different position, which is also corrupted but in a different direction. In the reality, the particle in the time interval between the measurements travelled a very short distance, but because of the light diffraction, the camera could record greater or smaller distance and consequently measure the higher velocity fluctuation. Based on the comparison of simulation results with experimental data, the problem of refractive index effect in the two-liquid mixing case has been already exposed in [12]. This was later confirmed by experimentalists [13] who observed that the distortions in the measurements grow with the distance between the measuring position and the camera lens. At the measurement of the velocity components this effect is not significant due to the averaging process. However, in the case of the measured velocity fluctuations this effect may lead to a significant experimental error, which cannot be easily estimated at the current state of the knowledge [13].
5. Conclusions
A numerical simulation of the evolution of the turbulent mixing layer in the horizontal confined two-component flow was performed with the OpenFOAM code, version v1606+.
The simulated velocity and the turbulence kinetic energy in the single liquid case obtained with the RANS and the LES simulations match the experimental results relatively well. The RANS simulations provided somewhat better agreement with the single liquid experiment due to the possibility of imposing the experimental inflow condition. Nevertheless, the LES approach has proved to be a very useful method when a better insight into the development of primary velocity fluctuations is required.
In the two-liquid case only the calculated velocity and the concentration profiles satisfactorily agree with the experimental data. On the contrary, the simulated velocity fluctuations and consequently the turbulence kinetic energy are by order of magnitude lower than experimental values, for both, RANS and LES simulations. The comparison of the single liquid and the two-liquid simulation results shows very similar velocity and velocity fluctuation profiles, since the density and the viscosity differences between the two mixing flows in the two-liquid case are very small, only 1% and 2%, respectively. Despite the small differences in the density and in the viscosity, the velocity fluctuations measured in the two-liquid experiment are by order of magnitude higher than in the single liquid experiment. A sensitivity analysis of modelling approaches was carried out in order to explain the deviation in the fluctuations. The simulation results indicated that the possible reason for such high disagreement may be attributed to the systematic experimental error. It was found that the effect of refractive index difference may lead to the experimental error in the measured fluctuations for the two-liquid case.
This study shows the tight relationship between the experimental and simulation results. The simulation models need to be validated by experimental data, but it can be also the other way around; the simulations can be used to verify the accuracy of the experiments.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge the financial support from the Slovenian Research Agency (research core funding no. 26 P2-0026 “Reactor engineering”).