Computational fluid dynamic simulations on liquid film behaviors at flooding in an inclined pipe☆

Jianye Chen ,Yuan Tang ,Wei Zhang ,Yuchen Wang ,Limin Qiu ,Xiaobin Zhang ,*

1 Institute of Refrigeration and Cryogenics,Key Laboratory of Refrigeration and Cryogenic of Zhejiang Province,Zhejiang University,Hangzhou 310027,China

2 China Aerodynamics Research and Development Center,Mianyang 621000,China

Keywords:Two phase flow Flooding Countercurrent flow limitation Computational fluid dynamic Liquid film Inclined pipe

A B S T R A C T The complex liquid film behaviors at flooding in an inclined pipe were investigated with computational fulid dynamic(CFD)approaches.The liquid film behaviors included the dynamic wave characteristics before flooding and the transition of flow pattern when flooding happened.The influences of the surface tension and liquid viscosity were specially analyzed.Comparisons of the calculated velocity at the onset of flooding with the available experimental results showed a good agreement.The calculations verify that the fluctuation frequency and the liquid film thickness are almost unaffected by the superficial gas velocity until the flooding is triggered due to the Kelvin–Helmholtz instability.When flooding triggered at the superficial liquid velocity larger than 0.15 m·s?1,the interfacial wave developed to slug flow,while it developed to entrainment flow when it was smaller than 0.08 m·s?1.The interfacial waves were more easily torn into tiny droplets with smaller surface tension,eventually evolving into the mist flow.When the liquid viscosity increases,the liquid film has a thicker holdup with more intensive fluctuations,and more likely developed to the slug flow.

1.Introduction

In counter-current two phase flow,with the increase of gas flow rates,the flow becomes unstable and finally part or whole of the liquid film reverses flow direction,defined as the onset of flooding.Flooding,as the counter-current flow limitation,is encountered in many industrial devices,such as heat pipe,reflux condenser,packed column and some nuclear reactor accidental scenarios.It often deteriorates the regular operation of devices and has been the subject of numerous investigations in the past decades.

The computational fluid dynamic(CFD)method is a potent tool as it provides more insight into the physics of complex two-phase stratified flows.Volume of fluid(VOF)method with the surface tracking technique and the Eulerian model are widely used to model the flows[1].Murase et al.[2]conducted numerical calculations of the 1/15th scale of pressurized water reactor hot leg with the VOF method.The results underestimated the water flow rates at the upper end of the inclined pipe and overestimated that in the horizontal segment at flooding.Later,they[3]improved the computational grids and schemes to model steam-water flows at 1.5 MPa under pressurized water reactor(PWR)full-scale conditions,and obtained results consistent with the Upper Plenum Test Facility(UPTF)data except in the cases of large steam volumetric flux.Overall,VOF method obtains limited success in modeling flooding because one set of N–S equations are shared by the two phases and the momentum exchange between them ignored.

The Eulerian model solves a set of Navier–Stokes equations for each phase and coupling is achieved through the shared pressure,interphase momentum exchange and energy exchange.The interphase drag force modeling mostly affects the calculation precision when flooding occurs.Wang and Mayinger[4]applied the interfacial friction factor model proposed by Lee and Bankoff[5]to model the UPTF and got a satisfactory result in spite of a little difference of the flow patterns with experimental results in the horizontal leg.Minami et al.[6]and Utanohara et al.[7]conducted CFD simulations on countercurrent flow in a 1/15th scale of PWR hot-leg model.The required interfacial friction correlations were selected from a combination of the available onedimensional experimental correlation for annular and slug flows that gave the best agreement with the experimental data.

However,a general geometry-independent model closer to physics and less empirical is a long-term objective for the drag force modeling for flooding simulation.The Algebraic Interfacial Area Density(AIAD)method was adopted by H?hne et al.[8–11]to model flooding phenomena in PWR.The results show its success in predicting the transition among different flow patterns.To validate the general usefulness,more work is still required in modeling the flooding in the complex channels where effects of gravity are important.

In this paper,a comprehensive numerical investigation on flooding in inclined pipe is carried out.To verify the accuracy of the numerical method,calculated superficial gas velocities at flooding are compared with the published experimental data.Emphasis is on the wave characteristics and the flow patterns during flooding.Moreover,impacts of surface tension and liquid viscosity on flooding are investigated.

2.Numerical Models

2.1.Control equations

In the calculations,we solve the independent mass and momentum conservation equations for the two phases,which have the following form:

where k refers to gas(G)or liquid phase(L).The drag force τD=which is derived from the interfacial shear stress,is most conveniently expressed in terms of the drag coefficient CD:

where ρLGis the average density,u is the relative velocity and A is the projected area of the control volume in the flow direction.In the AIAD model,CDhas different correlations in the full range of volume fraction of gas phase,and it allows the detection of the morphological form and the corresponding switching for each correlation from one object pair to another[12].The asymptotic limits of bubbly and droplet flows are improved by comparing different coefficients.Details can be found in Refs.[8–11].

CFD software Ansys Fluent 14.5 including the two-fluid model was used.Multi-Fluid VOF model for the Eulerian multiphase allows using the sharpening schemes Geo-Reconstruct,compressive,and CICSAM with the Explicit VOF option.This model overcomes some limitations of the VOF model due to the shared velocity and temperature formulation,and is often used for the cases requiring sharp interface treatment between phases.More details can be referred to[13].The AIAD model,which has successfully predicted the interfacial drag force in the countercurrent flow during flooding[8–11],is adopted in this work.

2.2.Turbulence closure

In the counter-current free surface flows,the high velocity gradient at the phase interface will generate high turbulence disturbance in both phases when using differential eddy viscosity models.Hence,turbulence damping is required in the interfacial area to correctly model such flows.For the two-fluid formulation,Egorov et al.[12]proposed a symmetric damping procedure based on the standard ωequation,which provides a solid wall as damping of turbulence in both phases,and is formulated by Wilcox[14]as follows:

which is added as a source to the ω-equation below:

2.3.Boundary conditions

Deendarlianto et al.[15–17]have carried out lots of flooding experiments in inclined pipes and got abundant experimental data.In this paper,the core computational domain of 1 m length and 16 mm I.D.is modeled as a two-dimensional(2D)structure by reference to the experimental setup[15–17],as shown in Fig.1.The same 2D model for two-phase pipe flow in a pipe has also been adopted by other researchers[18,19].Water inlet is conical in order to decrease the disturbance.The inclination is regulated by changing the gravitational direction.

The original grid in the whole domain consists of 95200 cells,as shown in Fig.2.The grid is doubled to demonstrate the grid independence.The calculated dynamic liquid film thicknesses(h)with a void fraction of 50%as the interface are compared with the experimental counterpart when the flow is simulated to the steady state.It indicates that the liquid film wave fluctuates regularly with almost the same amplitudes and frequencies on both computational grids,as depicted in Fig.3.Therefore,the original grid of 95200 cells is chosen for the following simulations.

The inlets and outlets are velocity inlets and pressure outlets,respectively.The two phases are set as adiabatic and incompressible.At the beginning,an initial liquid film is maintained at a constant superficial liquid velocity(UL)with a very small gas flow rate UG.Here,the superficial velocity is

where V is the volume flow rate and D is the inner diameter of the tube.Then,UGincreases with a small gap until partial liquid film reverses its flow,the moment is defined as the occurrence of flooding,and UGat this moment is called critical gas velocity(UGF).

3.Result and Discussions

3.1.Critical flooding velocities

Fig.4 illustrates the comparison of calculated UGFwith the experimental data[17]at different ULwhich as listed in Table 1.Overall,the simulated UGFpredicted the experimental UGFwith the accuracy of±25%,which shows the reliability of the simulation method.The deviation between the simulated UGFand measured UGFis attributed to the different definitions of the onset of flooding.In the experiments[17],the onset of flooding was identified by the maximum airflow rate at which the discharged liquid flow rate is equal to the inlet liquid flow rate.In the numerical simulations,the wave reversal in the close inspection of flow pattern variation was marked as the onset of flooding,the same criterion being also adopted by most of researchers[20–22].

3.2.Interfacial characteristics before flooding

Karimi and Kawaji[23]and Vijayan et al.[24]performed flooding experiments and reported that the liquid film thickness tended to be little affected by the gas flow before flooding,while it increased sharply when close to flooding.Luo et al.[25]also concluded that the effects of gas phase could be neglected at low flow rate through the experiments on inclined plates.In this simulation,a similar trend has been observed.As in Fig.5,at the upper part of the tube(x=0.5–0.8 m),the liquid accelerates downwards after entering the tube and the drag force does little effect on h at UG=0.08 and 0.75 m·s?1.When almost approaching flooding(UG=1.5 m·s?1),the drag force on the liquid film exerted by the gas flow increases significantly,which slows down the liquid film.Therefore,h increases at the upper part(x=0.5–0.8 m).It also indicates that all the cases with different gas velocity have a smoother film at x=0.5–0.8 m than the lower part,indicating that the interfacial waves progressively slow down when propagating to the liquid outlet,as a result,the amplitude becomes larger.Meanwhile,when close to flooding,as being impeded by large interface waves at the lower part,both the gas and liquid velocities drop at the upper part,making the interface relatively smooth.

Fig.1.2D calculation model of a pipe.

Fig.2.Local grid of the water inlet region.

Fig.3.The variation of film thickness with different meshes.

Fig.4.Comparison of flooding velocities for simulations and experiments[17].

The wall shear stress by the liquid phase is measured as an important index of turbulence effects in two-phase flow experiments.Pantzali et al.[20]reported that the measured fluctuation of τwclosely reflected the fluctuation of thin liquid layers in an inclined small diameter tube.The literature[26,27]indicated that the dominant frequency of τwis less than 10 Hz in vertical pipes.Drosos et al.[22]experimentally measured τwbefore flooding in vertical tubes and found the dominant frequency is between 7 and 10 Hz for the free falling film and approximate 2–4 Hz at larger UGnear flooding.While in the present calculations,the wall shear stress by the liquid phase was monitored at the location near the liquid exit(x=0.1 m).The dominant frequency seemsto be independent with UGand the value is about 3 Hz.An additional regular frequency peak occurs at approximately 5 Hz,which is probably caused by the significant growth of the smaller waves.This divergence between the experimental and simulation results is possibly due to the different boundary conditions.In the experiments,ReLof the tested liquid film is between 135 and 265,which is classified to laminar flow,while the liquid flow in Fig.6 is turbulent with ReLof 3976.

Table 1 Calculating conditions

Fig.5.The variation of the liquid film thickness with the increase of UG before flooding.

Fig.6.The FFT analyses of τw with different UG before flooding.

3.3.Interfacial characteristics at flooding

As UGreaches to about 2.9 m·s?1while keeping ULconstant of 0.25 m·s?1,an interfacial wave suddenly blocks the whole cross section of the pipe and then flows backwards under the drag force of gas,which is the onset of flooding.Fig.7 presents the formed slug from simulation and experiment[15]at the same condition.The streamline shows that only the wavy portion reverses its flow direction while the main stream keeps unchanged.At this time,the inlet pressure increases sharply because of the block of the pipe.Soon afterwards the slug is swept off by the great differential pressure,resulting in that partial liquid falls back to the liquid film and the left is blown off the air outlet.The pressure gradient falls down as the gas space becomes continuous again,and the corresponding pressure variation at the gas inlet is shown in Fig.8.It is found that the simulation pressure is generally lower than the experimental data,regardless that they have the same variation trend.This may be attributed to the different liquid inlet conditions of the experimental apparatus and the calculation model,as referred in Section 2.3.This difference may lead to the different interfacial behaviors and gas spaces near the liquid inlet.

Fig.7.Flow structures of the counter-current gas liquid flow for simulation and experiment[15]at the moment of flooding onset(UL=0.25 m·s?1,UG=2.9 m·s?1,θ=30°).

Fig.8.Air inlet pressures for simulations and experiments during flooding[15](UL=0.25 m·s?1,UG=2.9 m·s?1,θ=30°).

Fig.9.Flow structures of the counter-current gas liquid flow for simulations and experiments[15]at flooding(UL=0.03 m·s?1,UG=9.1 m·s?1,θ=45°).

For the case of UL=0.03 m·s?1,θ=45°,the liquid film is much thinner than UL=0.25 m·s?1,θ=30°.Flooding is triggered when UGincreases to 9.1 m·s?1.Fig.9 illustrates comparatively the flow patterns at flooding obtained from experiment[15]and corresponding CFD calculation.Different from the phenomenon in Fig.7,the interfacial wave typically flows backwards to the liquid inlet before it grows to block the whole cross section.The local gas velocity at the wave crest reaches maximum for the decreased space of gas,leading to the biggest shear stress on the wave crest.The gas flow tears the wave crest into numerous droplet entrainments.The thin liquid film is easy to entrap bubbles,which is observed in both the experiments[15]and simulations.Figs.7 and 9 show the different flooding mechanisms.Actually,the simulation results of all the conditions in Table 1 indicate that flooding acts as the local bridge when ULis bigger than about 0.15 m·s?1,while it acts as isolated waves with numerous entrainments when ULis smaller than 0.08 m·s?1.The conditions depend on θ when ULis between 0.08 and 0.15 m·s?1with water–air:it seems that at a lower inclination angle,flooding tends to act as the local bridge.The comparisons with the experimental results in Figs.7 and 9 illustrate that the simulation results are convincing.

Fig.10.Simulated wall shear stress with gas flow rate.

The temporal variations of τwat the same liquid Froude number(FrL)are qualitatively compared between the calculated data and the experimental data[20]in Fig.10.The dimensionless FrLis used to guarantee the same liquid inertial force effects for the different gas conditions,which is characterized by the ratio of inertia force to gravity of fluids is defined as

The simulation results show that with the increase of UG,and the oscillation of τwbecomes intensified.When UGreaches 6 m·s?1,τwfluctuates more drastically.At about 3.3 s when flooding occurs,τwintermittently falls to zero,then turns back to about 4 kg·m–1·s?2soon afterward.This fluctuation trend of τwis due to the changes of liquid velocity during flooding.As a whole,the simulated fluctuation of τwcoincides qualitatively with the experimental data.

3.4.Influence of surface tension on flooding

The influence of different surface tensions(σ)on flooding is investigated in the simulations.Deendarlianto et al.[16]investigated experimentally the effects on flooding in the aforementioned apparatus.A surfactant was added to vary the liquid surface tension without altering other liquid properties.In the simulations,liquids with surface tension of 0.072,0.054,and 0.034 N·m?1were investigated corresponding to the experiments.And they are labeled as S72,S51,and S34,respectively,as listed in Table 2.

Table 2 Liquid properties employed to study the influence of surface tension on flooding

3.4.1.Interfacial characteristics before flooding

The FFT analysis of τwwith different σ before flooding is shown in Fig.11.It indicates that the dominant frequency maintains approximately 3 Hz with different σ and UG.At the case of S34,much more smaller waves occur with higher frequencies,illustrating that smaller σ brings about more small disturbance to the liquid film.The amplitude of the dominant wave is bigger with smaller σ.When UGis near UGF,the dominant amplitude of S72 is much smaller than that of S51 and S34[Fig.11(a)].UGhas great influence on the amplitude in S51 case,while it affects little on the amplitude in S72 case.

Fig.11.FFT analysis of τW with different σ.

Fig.12.Flow patterns at flooding(S34,UL=0.25 m·s?1,UG=2.6 m·s?1,θ=30°).

Fig.13.Flow patterns at flooding(S72,UL=0.25 m·s?1,UG=2.9 m·s?1,θ=30°).

3.4.2.Flow pattern transition at flooding

Jepson et al.[28]held the opinion that surface tension is a stabilizing force for the interface.They found by experiments that the decrease of σ will bring about more entrainments during flooding.Sacramento andHeggs[29]found that the emergence of entrainments helps to trigger flooding.The simulation results of flow pattern variation of the case S34 at flooding are shown in Fig.12.It illustrates that droplets are easily formed from the liquid film,and then flow towards the gas outlet entrained by the gas flow.Some of the droplets converge to liquid lumps,which then fall back to the liquid film,disturb the liquid film and lead to a big wave[Fig.12(b)].Under this circumstance,flooding is triggered by the entrainments.As a comparison,the flow pattern variation of the case S72 presented in Fig.13 reveals that no liquid drops are torn off the liquid film at flooding.

Table 3 Liquid properties employed to study the influence of liquid viscosity on flooding

Fig.14.The variations of liquid film thickness with different viscosities.

Fig.15.FFT analysis of τw with different viscosities.

3.5.Influence of liquid viscosity on flooding

The influence of liquid viscosity(μL)on flooding is also investigated in the simulations,and the modeled cases are listed in Table 3.

3.5.1.Interfacial characteristics before flooding

Fukano and Furukawa[30]conducted experimental researches to study the influence of μ L on the liquid film in vertical tubes and found that μ L has a great effect on the liquid hold up,namely liquid film thickness(h).With the same UL,the bigger μ L leads to a higher h.In the simulations,the same conclusion has been drawn in the inclined tube as illustrated in Fig.14.The effect of μ L on the fluctuation of the liquid film is smaller with lower UG,compared with the case of V2 with higher UG.Literature survey reveals that the effect of μ L on the stability of the liquid film has not been completely clarified[31,32].On one side,large μ L helps the stability of the interface,therefore delays the occurrence of flooding;on the other side,large μ L will bring about larger h,correspondingly,narrow space for gas flow,leading to a larger real gas velocity.Consequently,the shear force on the wave crest by the gas flow increases,which means the critical unstable conditions caused by the K–H instability easier to reach.It seems that the ultimate results depend on the compromise of the two effects,and our simulations support the later opinion.

In regard to the effects of μ L on the frequency,as shown in Fig.15,with higher μL,τWhas bigger dominant wave amplitudes.μ L has a slight effect on the dominant frequency which is among 3–5 Hz.

3.5.2.Flow pattern transition at flooding

When flooding occurs,μ L plays an important role in the transition of flow patterns.For the case of V0.5,even at a high film thickness,the growing waves can hardly touch the top of the pipe.The shear stress on the wave crest by the gas flow overcomes the small liquid interior viscous force.The wave reverses and entraps some gas,leading to a big bubble enveloped by the wave,as shown in Fig.16(a).In contrast,considering the case of V2(Fig.16(b)),the wave can stably keep its form compared to the case of V0.5.The growing wave bridges the top of the pipe and blocks the whole cross section,then a large slug forms.

4.Conclusions

The CFD calculations based on coupled Eulerian model with Multi-Fluid VOF model was performed to investigate the flooding phenomenon in an inclined pipe.The simulation results validated that this model predicted well the flooding velocities under different inclinations and can give more insight to the film structures in gas-liquid two phase flow than experiments.

Two main flow patterns occurred at flooding:with a big superficial liquid velocity(0.15 m·s?1or larger),the interfacial wave developed to slug flow while it developed to entrainments with a small superficial liquid velocity(0.08 m·s?1or smaller).

Surface tension and liquid viscosity have important influence on the flow pattern variation:(1)surface tension acts as an stabilizing force on the liquid film and with the decrease of surface tension,the liquid film becomes more unstable and the secondary waves increase.While at flooding,the liquid film with lower surface tension brings about much more liquid entrainments,which disturb the liquid film in turn and trigger the occurrence of flooding.(2)With a higher liquid viscosity,the liquid film has a higher thickness and more drastical fluctuation at big superficial gas velocities close to flooding.Higher liquid viscosity helps the interfacial wave to keep its shape at flooding and form the slug flow during flooding.By contrast,with lower liquid viscosity the wave crest prefer to be scraped by the gas flow.

Nomenclature

Akinterface area density of phase k,m?1

B damping factor

CDdrag coefficient

D diameter of the pipe,mm

FrLliquid Froude number

Gωgeneration of ω,kg·m?3·s?2

h liquid hold up or liquid film thickness,mm

Δn typical grid cell size across the interface,m

Re superficial Reynolds number

Sωsource term for turbulence damping,kg·m?3·s?2

UGsuperficial gas velocity,m·s?1

UGFsuperficial gas velocity at flooding,m·s?1

ULsuperficial liquid velocity,m·s?1

ukvelocity vector,m·s?1

V volume Flow rate,m3·s?1

Yωdissipation of ω due to turbulence,kg·m?3·s?2

Fig.16.Comparison of flow patterns at the incipient flooding with different μ L.

α volume fraction

β κ–ω model closure coefficient of destruction term

Γωeffective diffusivity of ω,kg·m?1·s?1

θ inclination angle relative to horizontal,(°)

μ viscosity,kg·m?1·s?1

ρ density,kg·m?3

ρLGaverage density of liquid and gas

τDinterfacial shear stress,kg·m?1·s?2

τtturbulent shear stress,kg·m?1·s?2

τvviscous shear stress,kg·m?1·s?2

τwwall stress,kg·m?1·s?2