1. Introduction

This work concerns the evaluation of the capabilities of the CFX-11 software for the numerical predictions of gas entrainment in the case of a plunging jet configuration. The configuration of an impinging jet occurs in different scenarios of reactor safety analysis.

In the scenario of an emergency core cooling (ECC), water is injected into the cold leg. The pipe may only be partially filled with hot water, if a loss of coolant accident occurs. In this case, the injected cold water impinges as a jet on the surface of the hot water. Depending on the velocity of the jet, steam bubbles may be entrained below the surface by the impinging jet. These bubbles contribute to heat exchanged and mixing of the fluids. Heat transfer between cold and hot water and mixing in the cold leg play an important role since the mixed water enters the reactor pressure vessel and may cause high temperature gradients at the wall of the vessel. These gradients cause mechanical stress in the wall due to thermal shock, which can have a negative effect on the durability of the reactor vessel.

An impinging jet may also occur, when an emergency coolant tank is filled up with water and the initial water level is below the inlet. Here, the mixing of the injected water and the water in the tank is a point of interest if the temperatures or the boron concentrations are different.

Another scenario for the occurrence of plunging jet phenomena can be found in the case of a break, when insulation material of components is released by the break. The fibrous material is transported into the reactor sump and might there perturb the core cooling system. During this situation, the reactor sump is partially filled with water. The jet from the break impinges at the sump water surface and causes a fluid flow in the sump, which influences the transport of the fibrous insulation material towards the sump strainers. The gas entrainment and its influence on the fluid flow field and the transport of the fibrous insulation are of particular interest.

Generally for the CFD modeling of large hydrodynamic configurations with multiphase flow, the Euler-Euler approach is used. The physical process of bubble generation near a plunging jet occurs on a very small scale, which cannot be resolved in a mesoscale simulation. Therefore, the gas entrainment has to be physically modeled in simulations of plunging jets. The aim of this study is to find an approach for the simulation of plunging jet, where the gas entrainment can be deliberately tuned to some extent (e.g., in terms adjusting free parameters), in order for a physical model or correlation of the entrainment process to be implemented into future simulations.

In the plunging jet configuration, gas has two different morphologies (see Figure 1). The gas above the water level is a continuous phase, whereas the gas below the water level is bubbly, that is, a dispersed phase. The water can be regarded as a continuous phase everywhere. For modelling this with the Euler-Euler method, two approaches are possible.

Figure 1: Morphologies of the phases near an impinging jet.

One can use two different phases for the two morphologies of gas. Then, water is treated as a third phase. Gas entrainment near the jet and degassing at the water surface has to be modelled with sources and sink terms that describe the conversion of gas from a continuous to a dispersed (bubbly) morphology and vice versa. This requires algorithms that identify the regions of entrainment and of degassing.

The other approach uses only two phases, one for water and one for gas. The different morphologies of the gas then have to be reflected by different coefficients in the closures for the momentum transfer between the gas and water phases.

The first simulations presented here are performed with water as a continuous phase and gas as a dispersed phase. Thus, the gas is assumed bubbly everywhere in the domain, and a constant drag coefficient is applied. The influence of the magnitude of the drag coefficient is investigated. Then, a more complex drag model is tested, which take into account the different morphologies of the gas phase.

2. Definition of the Test Case

2.1. Geometry and Mesh

A cylindrical tank with a diameter of 100 cm is filled through a nozzle. The water level is 50 cm below the nozzle and 150 cm above the bottom of the tank. The nozzle diameter is 19 cm (see Figure 2). To reduce the costs of computation time, only a section of five degrees is used for the simulation to be performed as 2D axisymmetric calculation (see Figure 3).

Figure 2: Geometry of a cylindrical tank filled with water through a nozzle above the free surface.

Figure 3: Section (five degrees) of the cylindrical tank; geometry and mesh.

The structured mesh has 125 uniform cells for the total height of the domain. For the radius of the water inlet, seven uniform cells are used and 30 uniform cells for the opening (see Figure 3). The tank is quite small compared to the jet diameter and the height of the nozzle above the water level is small enough for physicality of the result to be influenced significantly by the walls. This disadvantage is accepted, since we concentrate here on the gas entrainment, which takes place where the jet hits the water surface. It can be assumed that effects far away from this area do not influence the gas entrainment. The limitations of the geometry and the low mesh resolution are meant to reduce the computational costs. This is important for parametric studies. The simplicity of the geometry is accepted here since this investigation is meant to study concepts for modelling the gas entrainment. For some of the cases, the grid is refined by reducing the cell size by a factor two in each dimension.

2.2. Fluid Properties, Initial Conditions, and Boundary Conditions and Turbulence

During the calculations, the fluids are water for the continuous phase and gas for the dispersed bubbly phase. The main properties (25^°C and atmospheric pressure) for water and gas are summarized in Table 1.

Table 1: Fluid properties.

The domain is partially filled with water up to a level of 150 cm above the bottom. The distance between the free surface and the water nozzle is then equal to 50 cm. The initial velocity for both water and gas in the computational domain is taken equal to 0 m/s in each direction. The hydrostatic pressure is initialized accordingly to the water level in the domain.

Liquid Inlet
The jet is injected through the nozzle with a velocity of 3 m/s. The volume fraction is 1 for water and 0 for gas.

Gas Outlet
For the gas outlet, an opening condition is used. The volume fraction is 1 for gas and 0 for water. Aconstant relative pressure equal to 0 Pa is assumed. For the fluid, a velocity normal to the boundary condition is considered.

Liquid Outlet
For the liquid outlet, an outlet condition is used. The volume fraction is 1 for water and 0 for gas. Therefore, the gas mass flow rate is equal to 0 kg/s at this boundary condition. For the maintenance of a constant liquid level, the liquid mass flow rate leaving the domain is defined equal to the liquid mass flow rate introduced by the injector.

Walls
Outer walls are adiabatic walls and are defined using a no slip boundary condition. For the “inner walls” caused by limiting the domain to a section, a symmetry boundary condition is applied. In the case of stratified flows, the buoyancy force causes a separation of gas and water.

Turbulence Model
The homogeneous shear stress turbulence (SST) model is applied (i.e., no separate calculation of the turbulence for both phases). In the ANSYS CFX-Solver modelling guide [1], a homogeneous turbulence model is recommended for separate flow and stratified flow, whereas for dilute dispersed two-phase flow (e.g., bubbly flow), the manual recommends using separate turbulence model for each phase. In the plunging jet, separate flow and bubbly flow coexist in one domain, so none of the turbulence approaches is suitable everywhere in the domain. The calculations presented below are calculated with a homogeneous SST model by default. For comparison, some calculations are repeated with an inhomogeneous turbulence model, which is the SST model for the liquid phase and a laminar assumption for the gaseous phase (see Section 3.3).

3. Dispersed Phase Model for Gas

The simplest approach for modelling the plunging jet is achieved if the water is treated as continuous phase and the gas is a dispersed phase with a constant particle diameter ( mm). This approach neglects the fact that the particle model is not appropriate for the gas above the water level.

3.1. Drag Model

For the bubbles, a constant drag coefficient is used. The default value used here is , which is the drag coefficient for solid spheres in the Newton range. To study the effect of the particle drag coefficient on the gas entrainment, the simulations are performed with the drag coefficient and with a reduced value for comparison.

3.2. Nondrag Forces

Since the focus of this study is on the gas entrainment at the surface, nondrag forces (lift force and turbulent dispersion force) are neglected here. It is expected that the turbulent dispersion force causes an increase of the horizontal extension of the bubble plume. Nevertheless, an application of nondrag forces above the water level is meaningless. Therefore, nondrag forces are not modelled here.

3.3. Results for the Simulations with Gas as Dispersed Phase

A few seconds after the jet release from the nozzle the interface becomes stable and the gas void fraction field also becomes steady (see Figure 4). A reduction of the drag coefficient from to has no significant effect on the gas void fraction field. Thus, the drag coefficient cannot be used as a parameter that influences the gas entrainment in the simulation. If the SST model is applied for the liquid phase and the turbulence of the gas is neglected (laminar assumption), the gas void fractions are similar to those in Figure 4 which have been calculated with a homogenous SST model. Therefore, the coupling of both phases by sharing the same turbulence field does not contribute to gas entrainment. In the subsequent simulations, only the homogenous SST model is used.

Figure 4: Gas void fraction fields for a simulation using a dispersed phase model for gas homogenous SST turbulence for both phases.

3.4. Vertical Gas Fluxes below the Water Level

For a characterization of the gas entrainment, it is advantageous to use integral quantities for the intensity and the geometry of the gas plume. By performing an extensive survey of experimental studies Bin [2] obtained correlations for the penetration depth and the entrainment rate. The penetration depth is the vertical extension of the gas plume below the water level. Bin's correlation [2] for the penetration depth in meter is where is the nozzle diameter in meter and is the vertical jet velocity at the water level (in m/s). Due to gravitational acceleration, the velocity of a free falling jet is increasing until it hits the surface. If is the height of the nozzle above the surface and is the liquid velocity at the nozzle, can be calculated as For the height of and  m/s, one obtains  m/s for the jet velocity at the water level and a penetration depth of 215 cm according to (1). The predicted value for the penetration depth is larger then the depth of the water in the tank. Therefore, the length of the gas plume might be restricted artificially by the geometry. In fact, according to Figure 4 the gas plumes almost reach the bottom of the tank. For a better quantification of the vertical distribution of the gas, the gas void fraction is integrated on horizontal planes: Since the gas void fraction is dimensionless, the integral (3) yields the dimension of an area for . This can be interpreted as the area occupied by gas on the horizontal plane . In Figure 5, is plotted versus the depth below the water level. Here, is normalized by the area of the jet cross-section at the inlet. There is only a little difference between the values for and (see also Figure 4).

Figure 5: Normalized gas area as a function of the depth below the water surface.

The depth at which the normalized gas area is zero can be used to define a penetration depth for jets. According to Figure 5, the penetration depth is 150 m which means that the plumes reach the tank bottom. This is in accordance with the prediction of (1).

The entrainment rate is the ratio of the gas flux entrained below the water by the impinging jet and the water flux of the jet. The correlation for the entrainment rate suggested by Bin [2] is where is the jet height above the water level and is the gravity. For the boundary conditions used in the simulations, this correlation yields an entrainment rate of 0.08. Another correlation was obtained by Ohkawa et al. [3]: This correlation yields an entrainment rate of 0.036 for the boundary conditions used here.

To compare the simulation results in terms of entrainment rate with the predictions given by correlations (4) and (5), the gas fluxes below the water have to be investigated more closely. The product of the gas void fraction and the vertical velocity of the gas defines a vertical gas flux density : The upward and the downward fluxes can be distinguished by the definition of So the total downward flux at a certain level below the surface is where is the horizontal cross-section of the domain at a certain level below the surface. The total upward flux is calculated in the same way. Figures 6–8 show the gas void faction, the vertical velocities, and the vertical gas flux density at a depth of 30 cm below the water level.

Figure 6: Gas void fraction as a function of the horizontal distance to the jet axis, 30 cm below the water level. The values for (blue) and (mangenta) differ hardly.

Figure 7: Vertical velocities as a function of the horizontal distance to the jet axis, 30 cm below the water level ().

Figure 8: Vertical gas flux density as a function of the horizontal distance to the jet axis, 30 cm below the water level ().

In Figure 9, the total upward and downward gas fluxes are shown for the two jets modelled with the drag coefficients and . The gas fluxes below the water level are normalized by the water flux of the jet at the nozzle and plotted as function of the depth below the water level. For , the upward and downward gas fluxes are similar which means that the solution is steady. At each level, the same amount is transported upwards and downwards. For the different drag coefficients, the gas downward fluxes are also similar. Thus, the drag hardly contributes to gas entrainment. The carry-under of gas therefore seems to be mainly caused by numerical effects within the solver. The curves for the normalized gas fluxes show a local minimum ca. 10 cm below the water level. This is the depth of the deformed water surface (“trumpet”) near the jet. At a depth of 60 cm, all the curves have a local maximum. This can be explained by the re-entrainment of bubbles, which are trapped in the vortex caused by the jet. At the depth of 20 cm, the normalized gas fluxes are about 0.06 which is just between the predictions for the entrainment rate by Bin (see (4)) and by Ohkawa (see (5)). Since the entrainment is mainly caused by numerical effects in this setup, we can expect the value to be sensitive to the geometry and resolution of the grid. This must be studied in future. However, it seems to be a coincidence that the simulated entrainment rate in this simulation is in the range predicted by empirical correlations.

Figure 9: Normalized gas fluxes as function of the depth below the water level. Blue symbols: upward gas flux (). Red symbols: downward gas flux (). Yellow: downward gas flux ().

Figure 10: The continuous phase density as a function of the gas void fraction.

Figure 11: The area density as a function of the gas void fraction.

4. The Simmer-Drag Model

In the previous simulations, the gas was treated as a dispersed phase everywhere in the domain. However, the SIMMER model, first introduced into the SIMMER-code [4], takes into account the distinction in morphology that phases can have in the domain. The morphology of the phases has to be reflected by appropriate parameters in the drag force. The magnitude of the force density for the drag is where is the drag coefficient, is the interfacial area density, and ρ is the density of the continuous phase (if the other phase is a dispersed phase). is the relative velocity between the two phases.

In the SIMMER model, the drag force depends on the gas void fraction . The gas is assumed to have the morphology of bubbles where the gas void fraction is low, that is, . Where the gas void fraction is high (), droplets are supposed to be present in gas. In the intermediate range (), a linear interpolation between bubble drag and droplet drag is performed. This means that the interfacial area density and the continuous phase density depend on the gas void fraction, which is used as indicator for the morphology of both phases.

4.1. The Continuous Phase Density

If the gas void fraction is low, the liquid phase is the continuous phase (). For high gas void fractions, the gas is the continuous phase (). In the intermediate range, the density is interpolated:

4.2. The Area Density

The total area density for a spherical particle is where is the particle void fraction. The drag coefficients for particles are related to projected areas. The projected area of a sphere is 1/4th of its total area. Therefore, the so-called projected area densities for spherical bubbles and droplets are calculated as where und are the bubble diameter and the droplet diameter, respectively. In the simulations for simplicity, the same particle diameters are applied for droplets (). Similar to the continuous phase density, the global area density a is defined as

4.3. The Drag Coefficient

Bubbles are assumed spherical, where the drag coefficient that has a constant value of 0.44 is applied. As an alternative, the Schiller-Naumann drag correlation is used which reads Since the material properties of the continuous phase are included in the Reynolds number, this equation yields two drag coefficients: for bubbles and for droplets. The drag coefficient at a position in the domain is calculated according to the local gas void fraction in the same way as it is done for the continuous phase density and the area density: If a constant drag is used, a case differentiation is not necessary. Then is applied everywhere.

4.4. Results

The development of the gas void fraction near the jet is studied using transient calculations. The drag is modified by applying either a constant drag coefficient or the Schiller-Naumann drag correlation. The calculations show a steady behaviour after a few seconds (see Figure 12).

Figure 12: Gas void fraction for a constant drag () and particle diameter  mm.

The influence of the drag model (constant drag versus Schiller-Naumann drag) and of the particle diameter is very low. The gas entrainment seems to be always overestimated, since gas void fractions higher than 60% occur below the surface in all simulations.

There is no free parameter inside the SIMMER drag model, which could be used to adjust the entrainment according to an empirical correlation or another physical entrainment model. The effect of modified drag coefficients has not been studied yet. However, using arbitrary drag coefficients causes unphysical velocities for buoyant particles (e.g., bubbles) and it is therefore meaningless.

5. The Algebraic Interfacial Area Density (AIAD) Model

5.1. Drag Model

The algebraic interfacial area density model applies two different drag coefficients, for bubbles and for free surface. The interfacial area density also depends on the morphology of the phases. For bubbles, the projected interfacial area density is where is the bubble diameter and is the gas void fraction. For a free surface, the interfacial area density is Since the concept of a continuous phase is not meaningful in the range of medium gas void fractions, instead of a continuous phase density, an average density is applied in (9). The average density is defined as where and are the liquid and the gas phase densities, respectively. In the bubbly regime, where is low, the average density according to (18) is close to the liquid phase density , which is the continuous phase density in this case. According to the flow regime (bubbly flow or stratified flow with a free surface), the corresponding drag coefficients and interfacial area densities have to be applied. This can be done by introducing a blending function which is 1 for bubbly flow and 0 for stratified flow. Then, the area density and the drag coefficient are well defined everywhere in the domain by It is not easy to find an algorithm that recognizes the flow regime of course. A very simple approach identifies the flow regime by using a gas void fraction limit . Bubbly flow is assumed, where , and stratified flow everywhere else. This would mean that blending function is a step function. To avoid numerical problems, a continuous blending function is preferred (see Figure 13): For a first judgement, the gas entrainment is quantified by the gas void fractions just below the liquid interface. These are investigated for various values of the free surface drag coefficient and gas void fraction limits . For the bubble drag coefficient, a constant value of is taken, based on the drag of rigid spheres at the medium to high Reynolds number regime. As bubble diameter  mm is chosen.

Figure 13: Blending function according to (21) for (black curve) and (red).

5.2. Variation of the Surface Drag Coefficient

It is not clear which surface drag coefficient is appropriate for the situation of the impinging jet. The value of has to include subgrid information of the free surface structure (“rough” or “smooth”), and this certainly depends on the grid resolution, since with a finer mesh more details of the surface structure are resolved. Therefore, the free surface drag coefficient is varied over several orders of magnitude. Its influence on the gas void fraction below the water surface is studied while keeping the gas void fraction limit constant at . Note that the vertical water velocity at the nozzle is kept constant at  m/s. The simulation is performed in the transient mode, but the result is almost in a steady state 10 seconds after the start when the jet is released from the nozzle.

Figure 14 shows the gas void fraction for . The gas entrainment seems to be overestimated here, since even at a depth of 50 cm below the water surface gas void fractions of 0.5 appear. In Figure 15(a), the corresponding bubble area density is displayed. Note that is proportional to the gas void fraction in (16), where it is greatest at . Of course bubbles are not assumed to be present where . According to (19) and (20), the blending function switches to free surface area density at high gas void fractions. The free surface area density for this case is shown in Figure 15(b). In Figure 16, the total area density according to (19), the total drag coefficient according to (20), and the product of and are shown.

Figure 14: Gas void fraction for the free surface drag coefficient .

Figure 15: (a) Bubble projected area density . (b) Free surface area density .

Figure 16: (a) Total area density . (b) Total drag coefficient . (c) Product .

As one can see from Figure 17 that the gas entrainment below the surface decreases if the surface drag coefficient is reduced. Note that the solver does not converge when . Since the maximal gas void fraction below the water surface is similar for and , it is obvious that entrainment cannot be suppressed by further reducing . Even for a very small bubble drag coefficient (), the gas entrainment is not negligible (see Figure 18). This indicates that numerical diffusion contributes to entrainment.

Figure 17: Gas void fraction for various surface drag coefficients. Representative plots at time seconds.

Figure 18: Gas void fraction for a very small bubble drag coefficient.

5.3. Variation of the Blending Function

By changing the gas void fraction limit the blending function can be modified. The value of has a significant influence on the gas entrainment (see Figure 19). It is not clear which value is appropriate since has no physical meaning. The definition of the blending function in general and of in particular is arbitrary to some extent. Note that in this case, the gas void fraction is used as criterion to identify the location of the surface. This is a quite simple approach of course and a more sophisticated blending function could use the gradient of to identify the surface since this gradient is high near the surface.

Figure 19: Gas void fraction for and various .

5.4. Grid Resolution

To check the influence of the grid resolution on the numerical solution, one calculation () is repeated with a doubled spatial resolution. The gas void fraction fields are similar (see Figure 20).

Figure 20: Gas void fraction for two different grid resolutions, , representative plots.

Figure 20 shows that the gas plume is narrower in the calculation with the higher resolution. By an integration of the vertical flux density at this level according to (8) and by normalizing the result with the water flux at the nozzle , the dimensionless entrainment rate is obtained. For the coarse mesh, the entrainment rate is 3.5% and for the fine mesh it is about 6.4%. Of course the results of a CFD model should not depend on the grid resolution (see [5]). However, in a simulation of an impinging jet with increasing resolution more details of the complex surface geometry at the impinging zone is resolved. In the borderline case of an infinite resolution, the real bubble generation process could be captured. With a decreasing resolution, the geometry of the impinging zone is further simplified. This is the reason why the resolution has an effect on the simulated entrainment.

6. Conclusion

Generally for the CFD modelling of large hydrodynamic configurations with multiphase flow, the Euler-Euler approach is used. This is the reason why the capabilities of the Euler-Euler approach for the modelling of the impinging jet are investigated in this contribution. The physical process of bubble generation near the jet occurs on a very small scale, which cannot be resolved in a large-scale simulation. Therefore, the gas entrainment has to be described by a model, which represents the physics of the entrainment (e.g., a correlation). For the implementation of such a physical model in the frame of a CFD code, a mechanism is required that allows the adjustment of the gas entrained in the simulation according to the correlation. Since the gas and liquid phases tend to separate due to the buoyancy force, it is an obvious choice to use its counterpart—the drag force—to obtain gas entrainment below the surface. However, if the gas is modelled as dispersed phase in an Euler-Euler simulation, the entrainment barely depends on the magnitude of drag coefficient, and it is obviously caused by numerical effects. Thus, this approach is not suitable for the implementation of a physical model for gas entrainment.

The SIMMER model assumes bubbly flow, where the gas void fraction is low and it assumes droplet flow, and droplets in gas, where the gas void fraction is high. A variation of the drag force—either by modification of the assumed particle diameter for bubbles and drops or by using different correlations for the drag coefficients of spherical particles—does not have a significant effect on the gas entrainment in the simulations performed with the SIMMER model. The gas entrainment is overestimated in all simulations, and there is no free parameter inside the SIMMER model, which allows the modification of the amount of gas entrainment in the simulation.

The algebraic interfacial area density (AIAD) model was found to be a suitable approach to adjust the entrainment. There are two free parameters inside the AIAD model that have a strong influence on the suction of gas across the liquid interface. These parameters are a drag coefficient for the free surface and the shape of the blending function. The blending function is used to identify regions of stratified flow (free surface flow) and regions of dispersed phase flow (bubbly flow) in order to apply the appropriate drag model.

Nevertheless, the gas entrainment calculated with the AIAD model is arbitrary, as the model does not realistically reflect physics of the bubble entrainment. Further investigations will be performed to improve the parameterization in terms of the AIAD model.(i)In the AIAD model, only the magnitude of the gas void fraction is evaluated by the blending function to identify the location of the free surface. In a more sophisticated approach, more criteria could be evaluated such as the gradient of the gas void fraction.(ii)It is not clear which drag coefficients for the free surface are appropriate. The drag coefficient should also reflect the roughness of the jet surface, for example.(iii)Up to now, the blending function is meant to identify the location of the free surface. With a more complex algorithm, it might be possible to identify the region where the jet entrains gas. This would allow applying special closure models (e.g., drag forces) to obtain a more controlled gas entrainment there.(iv)The parameters for a realistic entrainment probably depend on the grid resolution.(v)Up to now, the literature about gas entrainment near impinging jets is rather fragmentary. More experimental data are necessary to adjust the CFD models and obtain realistic entrainment in simulations. The Forschungszentrum Dresden-Rossendorf (FZD) is planning to perform new experiments with impinging jets. New sensors for multiphase flow measurements, which have been developed at the FZD, will also be used.

Acknowledgment

The NURESIM project is partly funded by the European Commission in the framework of the Sixth Framework Program (2004–2006).