MODELLING AIR AND WATER TWO-PHASE ANNULAR FLOW IN A SMALL HORIZONTAL PIPE

Numerical simulation using computational fluid dynamics (CFD) has been carried out to study air and water two-phase flow in a small horizontal pipe of an inner diameter of 8.8mm, in order to investigate unsteady flow pattern transition behaviours and underlying physical mechanisms. The surface liquid film thickness distributions, determined by either wavy or full annular flow regime, are shown in reasonable good agreement with available experimental data. It was demonstrated that CFD simulation was able to predict wavy flow structures accurately using two-phase flow sub-models embedded in ANSYS-Fluent solver of Eulerian–Eulerian framework, together with a user defined function subroutine ANWAVER-UDF. The flow transient behaviours from bubbly to annular flow patterns and the liquid film distributions revealed the presence of gas/liquid interferences between air and water film interface. An increase of upper wall liquid film thickness along the pipe was observed for both wavy annular and full annular scenarios. It was found that the liquid wavy front can be further broken down to form the water moisture with liquid droplets penetrating upwards. There are discrepancies between CFD predictions and experimental data on the liquid film thickness determined at the bottom and the upper wall surfaces, and the obtained modelling information can be used to assist further 3D user defined function subroutine development, especially when CFD simulation becomes much more expense to model full 3D two-phase flow transient performance from a wavy annular to a fully developed annular type.


Introduction
Two-phase annular flow in a horizontal pipe is of great importance in a wide range of industrial applications, including chemical process, heat transfer system, steam generator and pipeline transportation 1,2 .In the heat exchanger in particular, the state-of-the-art application is the evaporation and the condensation processes, where the void fraction of air and water changes during the heat transfer process.Hence, the two-phase flow behaviors can have profound effects on overall system device performance; because non-uniform wall heat transfer can produce uneven temperature distributions along the pipe wall surfaces, where water film distributions in annular flow region would be essential to improve local heat transfer efficiency 2 .It is widely known that two-phase flow in a horizontal circular pipe can produce different flow patterns, which is required to be clearly distinguished.This gives arise an attention for researchers/engineers to develop effectively numerical and experimental methods that are able to accurately predict twophase flow transient regimes.Hence, whole system performance can be improved via detailed design and optimization processes.
A crucial step for two-phase flow numerical analysis is to determine flow characteristics by using the most appropriate computational model that can correlate different flow regimes.Taitel and Dukler 3 developed a mechanistic model enabling to predict the transient of two-phase horizontal pipe flow using a flow regime map based on physical concepts.Mandhane et al. 4 developed their model using experimental data.The two-phase flow pattern and transition were also studied by Hewitt and Ishii 5 .They have made considerable efforts to develop theoretical models in order to correlate specific twophase flow regimes, in which various inflow conditions are encountered, such as flow rate, fluid properties and pipe dimensions, etc 6 .
Two-phase flow transition can produce a range of flow patterns from bubbly flow, slug flow and wavy flow to full annular flow.The flow transient process has a greatly influence on subsequentflow behaviour even at far downstream, as the flow entrainment and mixing are strongly dependent on the device configuration due to the slow exchange process of mass between vapour core and liquid film 7 .Figure 1 depicts a sketch of twophase annular flow pattern where the wavy liquid film is climbed onto the wetted wall surface with a gas core flowing along the center of the pipe, entrained with liquid droplets.The gravitation force has a major influence on the distribution of the annular film around the pipe walls where the liquid film is usually thicker at the bottom wall than that at the top wall.While the annular flow may wet the entire pipe walls, a dryout can occur at the upper pipe wall surface, where no liquid film coverage presents 4 .The dry out of the liquid film will significantly reduce the heat transfer efficiency.The present CFD studies seen in most literatures are mostly focused on wave annular or full annular flow region where the liquid will wet the entire pipe.increase of gas phase wave-front velocity and unsteady oscillating frequency with the increase of gas/liquid flow rates, etc. Mori et al. 9 also identified some large wave structures in a vertical pipe flow, noted as large disturbance waves that were found near the annular-churn flow boundary where higher liquid or lower gas superficial velocity is predominated.As a result, the flow kinematic viscosity was found inversely proportional to the wave frequency.Martin 10 suggested that wave velocity was not strongly dependent on the tube diameter at same gas velocity, and the wave frequency decreases with the increase of tube diameter.Unlike the vertical flow, the film distributions in a horizontal pipe are non-uniformly distributed along the pipe walls in nature 2 .
While numerous studies have been carried out to investigate two-phase flow in horizontal and/or vertical pipes, there is little understanding on the two-phase flow pattern and transition, particularly wavy annular flow behaviour in a horizontal pipe, owing to very limited empirical data.The experiments carried out by Schubring and Shedd 11 for a horizontal annular flow has led to the identification of two types of wave annular flow structures adhered to annular film:-the first being small wave ripples on the base film, moving at low velocity without carrying mass, and the second disturbing wave carries almost all the mass and travels along the tube at relatively high velocity, where larger film thickness characteristics were presented resulting in the interfacial shear related to film roughness.The scaling study 11 indicated that two-phase flow model is strongly correlated to the pipe diameter and the gas superficial velocity.Moreover, the study of Hewitt et al. 12 on investigation of liquid film structure in a 32mm diameter pipe has shown that a large number of air bubbles of all sizes are entrained into the water film, and they can be entrained and released continuously.This can lead to the formation of different film thickness characteristics under gravity effect.Subsequently, the liquid film thickness can have strongly influences on heat transfer and pressure drop.Styrikovich and Miropolski [13] found that there was a temperature difference between the top and the bottom area of the pipe and it varies with gas and liquid superficial velocities.
Jayanti and Hewitt 14 developed a 2D numerical model using computational fluid dynamics (CFD) to investigate a vertical varying film induced by an upward liquid and/or gas flow in annular flow regime.Portela et al. 15 studied the second flow affecting the liquid film distributions around a horizontal pipe using Eulerian-Lagrangian computational method, in which they proposed the secondary flow influence due to the droplets dragged along with its counter rotating cells, so that droplets were transported to the top wall of the pipe.An increase of droplet concentration presents in the gas core and the top wall regions of the pipe.Their model was used to simulate the particle-laden turbulence pipe flow where droplets were considered as small solid sphere to encounter the wall roughness, that inducing secondary flow patterns.Unfortunately, the liquid film characteristics were not fully provided in those studies.So far, very few authors have studied two-phase flow behaviors and liquid film distribution characteristics in a horizontal pipe using Eulerian-Eulerian CFD framework.This technical approach however has important implications to facilitate the development of an efficient flow system containing small diameter pipes, e.g. industry or domestic heat exchangers in which the ultimate goal is to carry out optimal design at a relative low cost and a short turn-around process time.
The primary objective of this study is to model flow pattern and transient process experimented by Schubring and Shedd 6,11 by using modern CFD approach to investigate air and water two-phase flow behaviours particularly the annular wavy flow distribution characteristics and the base film thickness.A user defined function (UDF) subroutine (ANWAVER-UDF) developed and implemented in ANSYS-Fluent solver under Eulerian-Eulerian framework will be applied to model the wave annular flow pattern.The results will be verified against available experimental data of Schubring and Shedd 6,11 .While proposed CFD method is general, only 2D case study will be carried out for parameter studies over a wide range conditions hereby in this paper.

Computational domain and mesh
Figure 2 presents a schematic view of 2D cross section through a straight horizontal pipe with all dimensions illustrated.The pipe is 1.0m long and 8.8mm in diameter, same as that adopted in the experiments of Schubring and Shedd 6,11 .Figure 3 displays sample structured mesh with a total of 100×88 grid points, uniformly distributed in the streamwise x-direction and highly stretched (with a stretch factor of 6) in the wall normal y-direction, respectively.The mesh quality is measured by max skewness value of 2.4×10 -4 , which is much lower than a criterion of 0.9 suggested by ANSYS mesh tool 16 .The Reynolds Number (Re) based on the airflow velocity at the pipe centre core is about 7700 in turbulent flow regime, while the water liquid film Re number is about 2785, in laminar transition to turbulence region.

Governing Equations of two-phase flow under Eulerian-Eulerian framework
The two-phase flow will be solved under Eulerian-Eulerian framework using a control volume (CV) method for mass, momentum and energy conservation equations implemented along with two-phase flow sub-models 16 .Mass conservation equation:

Inflow
where q v is the velocity of phase q and pq m  characterizes the mass transfer from th p to th q phase, and represents the mass transfer from phase 'q' to phase 'p'.'n' indicates a number of 'p' phase.In this simulation, 'p' represents primary phase of air and 'q' presents for the secondary phase of water, respectively.is source term of 'q' phase. is volume of faction of 'q' phase. is density of 'q' phase.The volume of fraction for all phases must equal to 1 17,18 , so that The concentration of vapour/liquid in two-phase flow is presented on the same coallocated meshes using vapour/liquid volume of fraction 17,18 .
where V q is the volume of the phase 'q' in a cell/domain, V is volume of the cell/domain.

Conservation of Momentum:
The momentum balance for 'q' phase that presents liquid secondary phase in this simulation: q q q q q q q q q q q n pq q lift q vm q pq pq qp qp where q  is stress-strain tensor the q th phase, q  and q  are the shear and bulk viscosity of phase 'q', q F  is an external body force, , lift q F  is lift force, and , vm q F  is virtual mass force if there is droplet entrained in continuous 'p' (gas) and 'q' (liquid) phases, pq R  is an interaction force between gas and liquid phases, and 'p' is the pressure shared by gas and liquid phases.pq v is inter-phase velocity, and its value is dependent on mass flowrate when 0 indicating of mass in phase 'q' is transferred to phase 'p'.

Energy conservation equation:
    1 : q q q q q q q q q q q q q n pq qp pq pq qp p p a h a u h a u q S t t where h q represents the specific enthalpy of the q th phase, q q  is the heat flux, q S is source term, contributed by sources of enthalpy, Q pq is intensity of heat exchange between p th and q th , and h pq is the inter-phase enthalpy.q u  is velocity magnitude of 'q' phase in x-direction.q q  is heat flux of 'q' phase.

Flow boundary conditions
The superficial velocities of air ( sg U ) and water ( sl U ) entering the pipe at the inlet are calculated via following equations.The G represents the total mass flux of all phases, whereas m  and x are the mass flow rate and air quantity, respectively 6,11 as where A is the pipe cross section surface area, the subscripts 'sg' and 'sl' represent superficial gas and superficial liquid, subscripts g and l stand for gas and liquid phases, respectively.The liquid droplet Reynolds number is calculated by where l  is the liquid dynamic viscosity.

Liquid annular film distribution
The initial liquid (water) wave structure and the base film thickness at pipe inlet was defined using a user defined function (UDF) subroutine in Fluent called ANWAVER model, where the maximum bottom and top base film thicknesses were set to 0.1938mm and 0.125mm, respectively.The model was coupled collaboratively with Eulerian-Eulerian two-phase flow model.Present simulation attempts to predict an annular flow pattern with the water film coated around the pipe circumferentially.The base flow films were then determined at the top and the bottom pipe walls.A time dependent wave function was introduced to create the liquid (water) wave structure and the base film developing mechanism as where t is iteration time step,  and b is wave function parameters associated with the time step and the base film thickness (Trio), respectively, and the latter is based on water/air superficial velocity.The averaged base film thickness is calculated using separate base film thickness measured at the top and the bottom pipe walls 6,11 as where b represents the bottom pipe wall; t the top pipe wall of tube and δ the base film thickness measured below the wave baseline at a distance to the pipe inlet.

Boundary conditions and simulation methods
The inlet air and water fluid superficial velocities are 13 m/s and 0.301 m/s, respectively with uniform distributions at a given temperature of 15 o C. The averaged relative pressure at the outlet is specified to 0, and the wall is adiabatic and no slip.The air and water at the inlet has a volume of fraction (VOF) of 0.104.
Simulations were carried out in transient mode to account for the buoyancy force effect on two-phase flow annular film thickness distributions along the pipe walls, as the air and water interfacial shear stress is strong.A user defined function (ANWAVER-UDF) model was adopted to model the transient process of wavy annular flow to full annular flow state, and the characteristics of the top and the bottom water film thickness being determined.A second-order phase-coupled numerical scheme was used to correlate pressure and velocity in an implicit manner.The turbulence model was standard k-ε model for both air and water phases.The air/water interface was modelled using high resolution interface capturing (HRIC) scheme to enhance the resolution of air and water interface.The flow was resolved iteratively using a small time step of 10 -5 with a total of 300,000 time steps to reach a residual of 10 -5 .Both pressure and water volume of fraction monitor points was defined at the pipe outlet to verify the convergence of the solution.

Results and Discussion
CFD can reveal strong dependence of water film thickness distributions with air and water superficial velocity, e.g. a wave annular flow feature can be observed at relative lower air and water superficial velocity, in contrary to a full annular flow scenario.
Figure 4 gives contours of air/water two-phase flow pattern at seven time instances.It can be seen that the flow transient occurs from plug/slug flow at time t = 0.434s (see Fig. 4a) to slug/wavy annular patterns at time t = 0.705s (see Fig. 4b).A wavy annular flow pattern is also observed at time t = 0.865s, 1.152s, 1.32s, 1.792s (see Figs. 4c-4f), respectively.Finally, a full annular flow pattern is developed at time t = 2.036s (see Fig. 4g).The transition of two-phase flow is also clearly presented.The asymmetrical wave flow pattern around the pipe top/bottom walls (see Figs. 4c-4f) may be attributed by the buoyancy induced flow instability and strong interfacial shear stress at the air/water interface.Zero water volume of fraction indicates a pure air occupation, which is mainly located at the centre of the pipe.The entrainment of air into water is predicted for the transient from plug/slug flow to wave annular flow patterns (see Figs. 4a-4b).It can be seen that the air 'pocket' is entrained at the front of water slug followed by the development of annular water film around the pipe top/bottom walls.The predicted tendency of air and water two-phase flow transition is consistent with the findings obtained by another similar study 13 .The initial water film at the pipe inlet was defined using wavy annular flow function (ANWAVER-UDF), where the bottom and the top base film thicknesses are set to 0.1938mm and 0.125mm, respectively.It was found that this film thickness changes with transient time and at full annular flow region (see Fig. 4g), it is thicker at the bottom wall than that at the top wall.Figure 7 displays CFD predictions compared to that of experimental data [6, 11] in terms of the liquid (water) film thickness ratio on the bottom/top walls (b/t) as a function of air superficial velocity.The water base film thickness is measured between the top/bottom walls and the air/water interface.For wave annular scenario, the water base film thickness is estimated between the top/bottom walls and the wave base line similar to the sine wave.It can be seen that the water base film thickness is higher at wavy annular flow region where the maximum bottom/top (b/t) film thickness ratio is determined at 3.6 compared to a value of 2.2 for a fully developed annular flow region.Figure 8 illustrates the water distribution spectra represented by the water volume of fraction for a wavy annular flow transient to a full annular flow.Figure 8a shows a full coverage of water film on the pipe bottom wall where the water VOF is equal to 1, that is consistent to those contours as seen in Figs.4c-4g.A gradually development of water film coverage on the top wall is also observed, where a full coverage of water film is predicted at the pipe outlet.Although the magnitude of water VOF at the centre of the pipe is low in the wavy annular flow region (about 3×10 -7 as seen in Fig. 8b), there is an implication of water moisture presented in the air gas core located circumferentially around the centre of the pipe.These moistures formed by small water drops can travel further upwards to wet the top wall surface, and thus a full coverage of water film can be predicted along the top wall of the pipe.

Conclusions
A two-phase flow CFD model has been employed to investigate water wavy annular flow transient to full annular flow in a horizontal small pipe using a commercial software ANSYS-Fluent, in which a specific user defined wave function (ANWAVER-UDF) is developed and implemented under Eulerian-Eulerian framework.CFD simulation has successfully modelled wavy annular flow structure and its transient process to a fully developed annular flow.The CFD predicted water film thickness distributions are in generally good agreement with the experimental data.The air/water turbulence kinetic energy and turbulence eddy dissipation rate are found higher in a full annular flow, compared to that of a wavy annular flow.The high turbulence kinetic energy may have contributed to high liquid (water) film coverage on the top pipe wall due to the wave induced droplets travelling further towards there.Consequently, this will lead to the top wall water film thickness increase for a full annular flow.Nevertheless, the water volume of fraction in the center of the pipe is found low for a full annular flow, compared to a wavy annular flow.The discrepancies of liquid (water) film thickness determined at the bottom and the top wall surfaces from CFD predictions in comparison to that of experimental data can be used as the baseline for further user defined function (UDF) subroutine development.This can eventually be used in conjunction with discretized particle model (DPM) to investigate the underlying mechanisms of the influence of air bubble and water droplet in wetting the pipe top wall surface.

Fig. 1 .
Fig. 1.A sketch of annular flow in a horizontal pipe.Azzopardi et al.8 studied two-phase flow liquid entrainment in a vertical pipe using experimental method.Their experiments have revealed some general trends, such as the

Fig. 2 .
Fig. 2. A sketch of 2D cross section of 1m long horizontal straight pipe of 8.8mm I.D.

Figure 5
Figure5illustrates contours of turbulence eddy dissipation rate of the air.It shows that the air turbulence eddy is dissipated significantly at the pipe top/bottom walls, where air/water interface appears to form a closed air core around the centre of the pipe.Near the top/bottom walls the transparent region has low air turbulence eddy dissipation rate, indicating that the water film mainly presents there with little air penetrating the interface.A wavy flow structure is predicted, and the water film distribution is more uniform and straight distributed at the top wall than that at the bottom, as depicted in Figs.5a, 5b.The high turbulence eddy dissipation implies a strong interfacial activity at air/water interface to balance high turbulence kinetic energy in the same region.A wavy annular flow to a full annular flow transient (see Fig.5c) occurs while turbulence kinetic energy increases with the increase of air/water superficial velocity (see Figs.6a-6b).In general, the base water film velocity is much slow compared to the air velocity.

Fig. 7 .
Fig. 7. Comparison of CFD predicted results and experimental data of Schubring and Shedd [6, 11] in terms of the ratio of the bottom (b) and the top (t) water film thickness as a function of air superficial velocity.

Fig. 8 .
Fig. 8. Water volume of fraction of wave annular flow: (a) at the bottom/top wall surfaces; (b) at the centre of the pipe.