Predictability and benefits of coupled Eulerian-Lagrangian approach over Eulerian characterization of droplet annular flow

In gas–liquid two-phase flow, the coexistence of a well-defined interface along with a random dispersion of phases possesses not only fundamental fluidic interest but also finds novel application in process, power and refrigeration industry. Out of different flow patterns in two-phase flow, a wavy-annular association of phases in macro-scale is often linked up with the generation of micro-scale droplet swarm in either adiabatic situation or dry-out related to heat transfer. Various researchers tried to address such multi-scale features in wavy annular and droplet flow using experimental observation (Hewitt and Nicholls 1969, Ishii and Grolmes 1975, Azzopardi 1986, Jepson et al 1989, Hazuku et al 2008, Wang et al 2013) and numerical simulations (Tryggvason et al 2006, Kumar et al 2016). The information, which came into pieces through unified scale experimental or numerical investigations, is reported next. Disturbance waves are one of the essential macro-scale characteristics often visually perceived in wavy-annular flows. These are made distinctive based on a magnitude of the relative velocity of phases to form waves. Many researchers (Taylor and Hewitt 1963, Azzopardi 1986, Sawant et al 2008, Zhao et al 2013, Alekseenko et al 2015) have conducted abundant experimental studies over the years to delineate its formation, departure and coalescence. The interaction of initial waves results into the development of disturbance waves, a fast-moving ripple (Alekseenko et al 2015). The frequency of disturbance waves diminishes due to coalescence (Azzopardi 1986). Gaseous phase Reynolds number and relative interfacial roughness are dominant parameters for evaluation of wave height in the tube. Wave height degrades with the increase in gas phase Reynolds number (Wang et al 2004).

The different set of momentum flux of liquid over the core gas eventually results in the development of a series of waves in annular flow configuration. At a high value of relative momentum flux, film thickness begins to deplete by entrainment of tiny droplets into the gaseous core. Mass, momentum and energy transfer in wavy annular flow is predominantly influenced by droplet entrainment rate. In diabatic annular flow, critical heat flux is significantly affected by entrainment rate and droplet size (Kataoka et al 2000). The height measured from the substrate surface termed as interfacial roughness along with the interfacial friction is the indispensable parameters for the entrainment rate. An increase in gas Reynolds numbers diminishes the interfacial roughness (Wang et al 2004). Primarily two different regions of entrainment in annular flow are the entrance and quasi-equilibrium regions (Ishii and Mishima 1989). For high film Reynolds number, droplet entrainment phenomenon occurs due to shearing off of undercutting wave. Experiments (Azzopardi and Hewitt 1997, de Bertodano et al 1997, Descamps et al 2008) have been done on the calculation of entrainment rate of droplets and proposed correlation (Kataoka and Ishii 1982, Ishii and Mishima 1989, Nigmatulin et al 1995, Zaichik et al 1998). Droplets hit the liquid film and increase the film thickness or rebound from the film depends upon the sticking efficiency and relative interfacial roughness (Šefko and Edin 2015). Droplet deposition can be broadly divided into two steps; i.e. transport of droplets and adhesion to the liquid film or wall. Study of droplet deposition and size distribution of droplet is vital for the calculation of the spatial location of dry out or burn out in dispersed droplet annular flow. An experiment conducted by Okawa et al (2005) for the mass transfer rate of a liquid phase in a vertical tube reported that mass transfer coefficient of deposition rate tends to decrease with increase in droplet concentration in the gaseous core. Later, Okawa and Kataoka (2005) have presented a correlation of deposition rate and found that at high concentration, the deposition rate is a function of droplet concentration and superficial gas velocity, whereas, at low concentration, the only superficial gas velocity is significant.

First, the idea behind the transition mechanism of scales is reported in the work of Ishii and Grolmes (1975). They described the transformation as a sequence of roll wave, wave undercut, liquid impingement, bubble bursting, and liquid bulge disintegration. Later, numerous studies (Ishii and Mishima 1989, Kuo and Cheung 1995, Kataoka et al 2000, Okawa and Kataoka 2005, Okawa et al 2005, Dasgupta et al 2017) have been done and reported in search of entrainment and deposition rates. But minimal studies have been focused in the direction of understanding such multiscale mechanism. Azzopardi (1997) has schematically described bag break up at low phase flow rate, and ligament breaks up at high phase flow rate. Liquid protrusion from the wavy film due to excessive liquid inertia or gaseous shear sheds of droplets which will contribute in the entrainment of liquid in the gaseous core. Kumar et al (2016) have performed the numerical simulation to explain the transition mechanism via the route of undercutting, rolling and orificing, in a vertically upward annular flow. Continuous entrainment of such droplets at different azimuthal and axial locations leads towards liquid dispersion at the core and thinning of a film at the wall. On the other hand, depending on post-birth dynamics, droplets can also deposit back to the liquid film and act against film thinning. Algebraic summation of entrainment and deposition dictates film thickness in the downstream section of a conduit carrying gas–liquid phases. In extreme condition, a film can extinct and give rise to dry-out. To avoid dry-out, knowledge about entrainment and deposition rates are of utmost importance. However, previous studies revealed the initiation of annular to droplet transition but have not targeted understanding of overall conversion of a well-defined interface to dispersed one. The major problem lies in the simultaneous and cumbersome capture of micro-scale droplet swarm along with the macro-scale wavy interface. The phenomenon can be termed as multi-scale signifying tracking of small-scale droplets as compared to wavy-pattern of a diametrical scale of the conduit. A proper understanding of both these interfaces requires either tremendous computational efforts or synchronized and seamless multi-scale model. The volume of Fluid (Hirt and Nichols 1981, Popinet 2009) and Lagrangian point particle (Tomar et al 2010, Ling et al 2015) based multi-scaling has been tried for predicting applications like gas assisted jet atomization (Ling et al 2015). Still, similar Eulerian-Lagrangian models are not yet proposed for droplet annular flow, till date.

This work described the complete association of droplets in the entire gas core along with the wavy annular interface and highlighting the advantage of a multi-scale model over the conventional volume of fluid (VOF) technique. A representative section of entrained droplets along with an annular film in the present study is shown in the left wing of figure 1. Here, though macro-scale wavy film co-exists with micro-scale isolated droplets, use of the unified model (volume of fluid; described in section 2) may adequately capture multi-physics of widely varying length scale in the domain. To do so, computational efforts also sore up. Present manuscript proposes the use of a coupled multi-scale model for efficient capture of the same flow physics, as described in the right wing of figure 1. In the next section, briefly, we have presented the coupled multi-scale model before discussion of the simulation results on droplets swarm generated from annular interfaces.

Zoom In Zoom Out Reset image size

A typical fluidic mechanism representing the coexistence of annular liquid film and dispersed droplets is analyzed numerically inside a vertically upward cylindrical tube. Adaptive mesh refinement technique based on gradient of VOF tracer is adopted to capture the interfacial dynamics meaningfully in the simulations using Gerris (Popinet 2003), a computational tool. Three-dimensional (3D) finite volume based Eulerian VOF, as well as seamless integration of VOF with Lagrangian point particle (LPP) approaches, are tried to capture such multiscale process. In this work, first, we have demonstrated the multi-scale mechanisms using widely varying mesh sizes in VOF framework and then mimicked the similar fluidic understanding using coupled VOF-LPP method. It is needless to mention that the latter approach is computationally efficient and requires lesser storage space. The ability of VOF based Gerris solver for prediction of smaller scale interfacial entity is well proven in the literature (Kolobov et al 2007, Fuster et al 2009, Agbaglah et al 2011, Aristov et al 2012, Tripathi et al 2015, Kumar et al 2016). In contrast, an isolated effort (Tomar et al 2010) has noticed to prove coupled VOF-LPP module of Gerris as an efficient replacement of VOF scheme, specifically for atomization issues.

In this section, we have briefly introduced both the algorithms before proceeding towards the understanding of the fluidic features of the phenomena. The numerical test setup is made in a cylindrical tube with diameter D and length L (schematically shown in figure 2(a)) to begin with annular liquid entry adjacent to wall and gas co-flow in the core. A typical 24 mm diameter and 144 mm length of the cylindrical domain containing air-water pair has been considered for the simulations. Results of the simulations are presented in non-dimensional length and time scales with mention of velocity and momentum flux ratio. Temporally adaptive meshes are used to discretize the domain, where the quadtree structure of the discretization is used to adaptively follow the small structures of the liquid film and surface waves, thus concentrating the computational effort on the area where it is most needed. This is done using an object class, GfsAdaptGradient in Gerris. Gradient of tracer criteria has been used to determine where refinement is needed. It uses a 'cell cost' defined as the norm of the local gradient of a given variable or function multiplied by the cell size. If the cell cost is more than a set value near the waves of liquid surface the mesh is refined to a 'minlevel' using quadtree structure. Flow configuration is set from bottom to top against gravity. As a representation, the thickness of adhered film with the wall has been set to D/8. In the present study, we have taken the film thickness of the wall adhered liquid film thicker than reported experimental observations (Azzopardi and Hewitt 1997, Cherdantsev et al 2014). It has been considered to deal with the numerical complexities of resolving the thin film flows and the computational expenses. The thin films appear as dry patches and break into several rivulets due to poor meshing if the disproportionately thin film is considered within the computational limit. Such rivulets' dynamics are entirely distinct from the continuous film flows due to the contact angle issues. These thick films would have resulted in flooding towards the gas core with less liquid velocities. Therefore, the liquid velocities are taken one order higher to have a smooth interface with the upward liquid flow when gas velocities are less.

Zoom In Zoom Out Reset image size

Following dimensionless numbers are involved in conducting the numerical study; velocity ratio ${({v_l}/{v_g})}$, non-dimensional time (t* = t√(σ/Δρδ³)), non-dimensional location (i* = i/D; $\forall i \to x,y,z$) momentum flux ratio ${ }\left( {J = { }{\rho _l}{v_l}^2/{\rho _g}{v_g}^2} \right)$, Reynolds number $({\text{R}}{{\text{e}}_t} = {{{{\rho _t}{j_t}D}}/ {{{\mu _t}}}};\forall t \to l,g){ }$, and Weber number $\left( {{\text{W}}{{\text{e}}_t} = { }\frac{{{\rho _t}{j_t}^2D}}{\sigma };\forall t \to l,g} \right)$. Table 1 shows the ranges of values tested in the present context.

Table 1. An account of non-dimensional parameters and their range, used in the present numerical study.

${v_l}/{v_g}$	$J$ × 10−5	${\text{R}}{{\text{e}}_l}$ × 10−5	${\text{R}}{{\text{e}}_g}$	${\text{W}}{{\text{e}}_l}$ × 10−3	${\text{W}}{{\text{e}}_g}$
7–40	0.4–13.3	0.35–2.1	500–1000	4–25	0.02–0.04

Momentum flux is one of the dominant factors along with geometrical configuration, causing the coexistence of annular film and droplet swarm at the core. Momentum fluxes for both the phases are tabulated in Kumar et al (2016) to show rolling (${v_l}/{v_g}$ ≪ 1), orificing (${v_l}/{v_g}\,$ ≈ 1) and undercutting (${v_l}/{v_g}$ ≫ 1) behavior. They have elaborated short-time dynamics of localized interfacial deformation to describe these fluidic features. It has been well established by Azzopardi (1997) that the drops are not created from the entire film interface, but specifically arise from the interfacial features, referred to as disturbance waves. It also needs to be mentioned that two main mechanisms (bag breakup and ligament breakup) have been identified as the cause of droplet entrainment depending upon the inlet flow parameters entrainment (Ishii and Grolmes 1975, Azzopardi 1983, 1997, Azzopardi and Hewitt 1997, Badie et al 2001, Lecoeur et al 2011, Cherdantsev et al 2014, Pham et al 2014). However, Cherdantsev et al (2014) have experimentally obtained that both the mechanisms may occur under the same flow conditions and on the same disturbance waves. Furthermore, the reasons or parameters separating these two ways of droplet integrations, as obtained by Cherdantsev et al (2014) is also not consistent with Azzopardi (1983).

It can be noted that always there should be extreme deformation of the interface for the entrainment to occur. Such deformations of the interface are dependent upon the forces acting on the crests of the interfacial waves, localized flow conditions around them and the shape of the interface. The type of mechanism of droplet disintegration will also be dependent on all these features. Hence, for resolving the dynamics of entrainment from the fundamental aspect, one needs to consider the effect of the spatiotemporally varying local flow conditions on the interface deformations. The present analysis has been reported for different velocity ratios (${v_l}/{v_g}$) of phases in the undercutting regime. Fluidic situations are shown at moderate ratio of liquid to gas momentum flux and adequate film thickness so that with the existing computational resources simulations can be completed within allowable time without missing any fluidic physics in the process. In actual industrial annular flow situation, with much thinner wall adhered liquid film and real flow situations, finer droplets will be produced at a faster rate from many protrusions of the liquid film. We have chosen the domain in such a manner so that the effect of flow conditions supplied at the inlet should be prominent throughout the domain. Then, at the inlet, various flow conditions are supplied to include the relative importance of the momentum exchange between the phases. The high liquid velocities are also supplied because Ishii and Grolmes (1975) have obtained a region on the transition plot. They indicated that liquid flowing as a film also contributes to the momentum exchange between phases. It is to be noted that considered flow velocities may not sustain annular regime but will be a perfect case study for establishing the benefit of the hybrid method over the VOF based simulations. In the early stage, the inlet momentum of liquid due to its high velocity will continue climbing up the wall instead of filling the pipe. Soon annular structure will be re-organizing the artificial interface configuration to an area filling regime based on superficial velocities of fluids. In the subsequent section, we have discussed the direct numerical simulation approach used for spatial and temporal fully resolved simulation of the present study. It needs to notice that we have taken care to accommodate at least six cells across the least diametral plane of the entities for accurate capturing of interfacial dynamics. It needs to be noted that present numerical methodology does not consider turbulence models for reported phenomenon of film to droplet transition. A better model including turbulent calculations may extract more correct interfacial configurations. However, higher order numerical simulations will be complex and costly than the reported one.

Incompressible mass conservation and two-phase momentum balance equations incorporating surface tension are used for modeling the fluidic interactions. Towards comprehension to the readers, a brief outlines of the governing equations are given in equations (1) and (2) as:

$$\begin{equation}\nabla .{{\boldsymbol{u}}} = 0\end{equation} \tag{ 1 }$$

$$\begin{equation}\rho \left[ {{\partial _t}{{\boldsymbol{u}}} + {{\boldsymbol{u}}}.\nabla {{\boldsymbol{u}}}} \right] = - \nabla p + { }\nabla .\left( {2\mu {\mathbf{D}}} \right) + \sigma k{\delta _s}{\mathbf{n}} + {{\mathbf{F}}_{\mathbf{b}}} + {{\mathbf{F}}_{\mathbf{p}}}\end{equation} \tag{ 2 }$$

where, ${{\boldsymbol{u}}} = u\hat i + v\hat j + w\hat k$, $\rho$ and $\mu$ represents fluid velocity, density and viscosity, respectively. D represents deformation tensor, which can be written as ${{\boldsymbol{\,}}}{{{\boldsymbol{D}}}_{ij}} = \,\frac{1}{2}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)$, and $\,\sigma$, k, ${\delta _s}$ and n are representing surface tension coefficient, surface curvature, Dirac delta function and unit normal vector respectively. ${{\mathbf{F}}_{\mathbf{b}}}$ is the body force and ${{\mathbf{F}}_{\mathbf{p}}}$ is the source term in momentum equation of a mesh due to particle interaction. For VOF-LPP method only, ${{\mathbf{F}}_{\mathbf{p}}}$ attains a nonzero value, discussed later on.

Density and viscosity are calculated by volume-fraction weighted average equations (Popinet 2003) such as,

$$\begin{equation}\rho \left( {\text{c}} \right) = {\rho _1}{\text{c}} + {\rho _2}\left( {1 - {\text{c}}} \right)\end{equation} \tag{ 3 }$$

$$\begin{equation}\mu \left( {\text{c}} \right) = {\mu _1}{\text{c}} + {\mu _2}\left( {1 - {\text{c}}} \right).\end{equation} \tag{ 4 }$$

Phase fraction (c) is advected following VOF equation as:

$$\begin{equation}{\partial _t}{\text{c}} + {{\boldsymbol{u}}}.\nabla {\text{c}} = 0\end{equation} \tag{ 5 }$$

where the volume fraction, c defined as:

$$\begin{equation}{\text{c}} = \left\{ {\begin{array}{*{20}{l}} {0},& {{\text{ cell}}\,{\text{filled}}\,{\text{with}}\,{\text{air}}} \\ {1},& {{\text{ cell}}\,{\text{filled}}\,{\text{with}}\,{\text{water}}} \\ {0 \lt c \lt 1}, & {{\text{ cell}}{\text{contained}}\,{\text{interface}}} \end{array}} \right..\end{equation} \tag{ 6 }$$

A finite volume based staggered grid discretization (Popinet 2009) has been applied over the domain for numerical simulation. The Classical time-splitting projection method is used to frame modified momentum equation by decoupling pressure field from equation (2). Velocity advection (${{\mathbf{u}}_{n + 1}}.\nabla {{\mathbf{u}}_{n + 1}}$) and viscous term (${{\boldsymbol{D}}}$) are discretized by Bell et al (1989) second order unsplit upwind and by Crank-Nicholson second order accurate scheme (Popinet 2003), respectively. Scheme to calculate ${{\mathbf{F}}_{\mathbf{p}}}$ is discussed later. The numerical route can be broadly categorized in the following steps:

• (a)  
  In the first step, the auxiliary cell centered velocity field (${\mathbf{u}}_*^c$) is calculated from the momentum equation by decoupling the pressure field.
• (b)  
  Then, the auxiliary face centered velocity field (${{\boldsymbol{u}}}_*^{\;f}$) is computed, from quad/octree mesh domain, by taking an average of cell centered values linked with every faces to achieve consistency of fluxes (Popinet 2003).
• (c)  
  The divergence of the auxiliary face velocity field ($\nabla .{\mathbf{u}}_*^{\;f}$) is computed by finite volume approximation as shown in equation (7):

  $$\begin{equation}\nabla .{\mathbf{u}}_*^{\;f} = \,\frac{1}{\Delta }\mathop \sum \limits_f {\mathbf{u}}_*^{\;f}.{{{\boldsymbol{n}}}^{\;f}}\end{equation} \tag{ 7 }$$

  where, ${{{\boldsymbol{n}}}^{\;f}}$ is a unit normal vector to the face and $\Delta$ is cell size of the control volume.
• (d)  
  The divergence of the auxiliary face velocity field is substituted in pressure-Poisson equation to compute face-centered divergence of auxiliary pressure correction field (${\nabla ^{\;f}}p$).
• (e)  
  Finally, updated value at consecutive time step for face and cell centered velocity fields are being calculated using corrected pressure.
• (f)  
  An updated field of advected density and viscosity have been calculated using volume-fraction-weighted average equations (3) and (4).
• (g)  
  Using piecewise-linear geometrical VOF scheme (Popinet 2009), updated volume fraction value is calculated.
• (h)  
  The third term in momentum equation has a curvature (k) and surface tension ($\sigma$), which is calculated by using height function model (Bell et al 1989) and parabolic reconstruction whenever radius of curvature is smaller than five cells. Surface tension is calculated by the balanced forced method (Popinet 2003), to use in pressure decoupled momentum equation.

In the coupled Eulerian-Lagrangian method, net external forces acting on the particle from the surrounding fluid is added as a source term (${{\mathbf{F}}_{\mathbf{p}}}$) in the momentum conservation equation (equation 2). The particle position and its velocity are governed by equations (8) and (9), respectively (Tomar et al 2010)

$$\begin{equation}\frac{{{\text{d}}{{{\boldsymbol{x}}}^k}}}{{{\text{d}}t}} = {{{\boldsymbol{v}}}^k}\end{equation} \tag{ 8 }$$

$$\begin{equation}\frac{{{\text{d}}{{{\boldsymbol{v}}}^k}}}{{{\text{d}}t}} = { }\frac{1}{{\rho _{\text{p}}^k{V^k}}}\left[ {{F_{{\text{drag}}}} + { }{F_{{\text{inertia}}}} + { }{F_{{\text{add}}}} + {F_{{\text{lift}}}} + { }{F_{{\text{ext}}}}} \right]\end{equation} \tag{ 9 }$$

where, ${ }\rho _{\text{p}}^k,\ {{{\boldsymbol{x}}}^k},$ ${{{\boldsymbol{v}}}^k}$ and ${V^k}$ are denoting density, position, velocity vector and volume of the kth particle, respectively. Remaining terms are the forces, given at the right-hand side of equation (9) incorporating the effect of drag, inertia, added mass, lift and external means respectively. The procedure to do the calculations of these forces are mentioned in Tomar et al (2010). As time marches, the velocity (${{{\boldsymbol{v}}}^k}$) for each particle is calculated using the first order explicit forward-Euler scheme. The source term (${{\mathbf{F}}_{\mathbf{p}}}$) is as follows:

$$\begin{equation}{{\mathbf{F}}_{\mathbf{p}}} = \mathop {\lim }\limits_{{V_{\text{f}}}} \sum\limits_{{N_{\text{p}}}}^{k = 1} {\frac{{{V^k}}}{{{V_{\text{f}}}}}} \left[ {\rho _{\text{p}}^k\left( {g - \frac{{{\text{d}}{{\boldsymbol{v}}^k}}}{{{\text{d}}t}}} \right) + \rho \left( {\frac{{{\text{d}}{\boldsymbol{u}}}}{{{\text{d}}t}} - g} \right)} \right]\end{equation} \tag{ 10 }$$

where, ${N_{\text{p}}}$ is denoting the number of particles in the fluid control volume, ${V_{\text{f}}}$. Equation (10) encompasses the reaction forces caused by the particle and the buoyancy effect. To maintain the numerical convergence, the discontinuous field of ${{\mathbf{F}}_{\mathbf{p}}}$ is smoothed using a Gaussian distribution with a standard deviation ${\sigma _{\text{p}}}$ (Tomar et al 2010). The functional form of smoothing effect, ${\overset{{\smash{\scriptscriptstyle\smile}}}{F} _{\text{p}}}$ is as follows:

$$\begin{equation}{\overset{{\smash{\scriptscriptstyle\smile}}}{F} _{\text{p}}} = \frac{{{{\mathbf{F}}_{\mathbf{p}}}{\text{exp}}\left( { - \frac{{{{\left| {x - {x_{\text{p}}}} \right|}^2}}}{{\sigma _{\text{p}}^2}}} \right)}}{{{{\left( {\sqrt {2\pi {\sigma _{\text{P}}}} } \right)}^3}}}\end{equation} \tag{ 11 }$$

where, the value of σ_p has been chosen as maximum between particle radius and cell size. Choice of σ_p helps to improve the runtime in a coarse population. Effect of σ_p becomes significant on the dynamics of the droplets in case of dense particle association.

In the present work, using coupled Eulerian-Lagrangian approach, a small droplet is replaced by LPP, which eventually places the particle at its centroid. Meshless LPP reduces the computational time by avoiding an inevitable time-taking route to capture the smaller droplet by generating more meshes in the domain. The qualifying criterion for the conversion of any detached entities/droplets into an LPP is 512 cells. The number is chosen based on a visualization study of corresponding Eulerian simulation so that only spherical or similar shaped drops are converted into particles. The selection of the threshold criterion is chosen in such a way that ensures only spherical, spheroidal entities qualify for particle conversion, and no lamellas will convert into the particle. After conversion of drop to particle, as the VOF interface is not existing surrounding the particle, mesh coarsens around it. On the other hand, in complete VOF scheme such coarsening does not happen. This coarsening reduces the presence of total number of cells in the domain and lessens runtime. Reduction in a number of cells as a threshold will disqualify a few larger drops for conversion into particle whereas an increase of the value of the threshold may promote conversion of lamellas into the particle. Net mass and momentum of the droplet are assigned to the particle. The space occupied by the droplet is replaced with matrix fluid. The velocity of the particle is assigned with the average velocity of the cells in which earlier VOF interface exist. At the same time, coarsened mesh after conversion of VOF to LPP gets the velocity depending on averaging scheme of its child cells. It needs to note that the droplet internal dynamics is not captured correctly, but its kinetics is modeled efficiently in the process.

In the next section, we have demonstrated the implication of both VOF and VOF-LPP method for droplet annular flow. In this process, we have also highlighted the computational efficiency achieved by introducing the Lagrangian particle. In this phenomenon, the macro-scale annular film is disintegrated into an entrained droplet that is tiny. For the efficient prediction of the small droplet, a proper resolution of interface structure with a sufficient number of meshes in the VOF method is utmost important. In figure 2(b), a sectorial cross-sectional view has been depicted to demonstrate mesh sufficiency in arbitrary shaped three-dimensional disconnected interfaces. Adaptive mesh structure in the figure shows at least six cells in a diametrical plane of the droplet. Entrainment protrusions obtained in the present study are a quite interesting feature which subsequently ejects a droplet. A similar mechanism has also been proposed by Wang et al (2013), in figure 2(c), a side by side comparison showing the predictability of the proposed model in the direction of droplet annular flow.

Numerical simulations are performed to observe the interfacial waviness in the annular film, which sheds of numerous droplets after subsequent pinch-off of liquid in a gas core. In a case of higher liquid velocity as compared to gas, interfacial waves protrude in the form of a bag-like structure. Though rolling is the predominant route for droplet annular flow regime, higher liquid velocity than gas is chosen in the present study to show the equivalence of VOF & coupled multiscale scheme. The choice is solely based on the reduction of calculation effort. Actual transition mechanism with multi-scale and VOF model is quite a time consuming and will be described in the future attempt. Figure 3(a) shows liquid protrusion which entraps gas beneath it. An upward moving gas stretches the lamella and thins down it to promote towards the pinching-off drop. Representative pinch-off of a droplet from a lamella at $\frac{{{v_l}}}{{{v_g}}} = 20$ is shown in figure 3(b).

Zoom In Zoom Out Reset image size

Shedding liquid droplet from the film is observed as entrainment, which is the fundamental step behind the thinning of the annular film, leading to droplet swarm. After disintegration droplet moves as per its velocity and appears as disperse entity in the gaseous core. Depending on its inception route, the droplet may join a similar population at the core or merge back in the liquid film. The latter situation is a so-called deposition. The liquid droplet journey at successive time interval is shown in figure 3(c), which depicts that droplet continues to grow a distance from its originating lamella. Perceived stages of droplet production from the annular film at ${{\text{v}}_{\text{l}}}/{{\text{v}}_{\text{g}}}{ } \gg 1$ are known as undercutting. A detailed description and flow physics around the produced droplet are mentioned in Kumar et al (2016).

Since the interface is subjected to relative momentum flux of phases, several irregular liquid protrusions can be seen frequently in well-separated annular flow. In the present study, we noticed these irregular protrusions are not axisymmetric, and a similar is reported in Dasgupta et al (2017). With time, some of the tiny irregular protrusions (ripples) merge each other and forms a big bulge of liquid which is termed as 'disturbance wave'. The amplitude and associated motion of disturbance wave are larger in comparison with ripples (Dasgupta et al 2017). The accurate prediction of entrainment solely depends on the characteristics of such waves. Figure 4(a) shows the temporal growth of amplitude of disturbance wave for various set of velocity ratios. Here, one can quickly notice the slope of the curve is least at the beginning stage, but a gradual rise in later stages. It depicts that the growth rate in the amplitude of disturbance wave is non-linear, and such non-linearity is continued for the rest of the cases. The rate of accumulation of liquid is a crucial parameter for the growth of such waves. Hence, by increasing the momentum flux of liquid film by keeping other parameters same, a sudden rise in the amplitude is observed as shown by steeper curves at high-velocity ratios.

Zoom In Zoom Out Reset image size

Further, we continued to explore the traversal of such a wave across the longitudinal direction, as shown in figure 4(b). Here, the trend, denoting the longitudinal instances of the disturbance wave is linear for each test conditions. It has established the analogy behind the velocity of such waves. In summing-up, the momentum flux of liquid associated with such waves is directly proportioning with the velocity ratios. Next, efforts are made to explain the dynamics of entrained droplets in the core.

In a pipeline carrying annular flow, produces droplets, is a multitude in nature from the different azimuthal plane, lamellas form and shed droplets at several axial heights. As a whole, if one looks at the three-dimensional pipe, an enormous number of droplets can be seen to be originated. As time progresses, the droplet population occupy the entire gaseous core creating a well-dispersed situation. Tracking of well-resolved droplet dynamics for such a large number is quite expensive in VOF framework. Moreover, holistically, droplets individual dynamics will not contribute much to the characterization of a dispersed population. Only the droplet trajectories are essential which can be well tracked by replacing it with an LPP, immediately after origination. Necessary governing equations for particle motion are shown in the previous section. In figures 3(d) and (e) we have seen the coupled VOF-LPP simulation for the same droplet, mentioned in figures 3(b) and (c). It can be comprehended from the comparative figures that the transformation of the droplet into particles is not biasing the trajectory. One may argue that the evolution of the droplet shape is not properly tracked by replacing it with a spherical particle. But these small-scale shape evolutions are having minute effects on the behavior of the droplet population.

To capture it in a better way, we have shown both the droplet trajectory and the corresponding particle in the same plane in figure 5. Dynamics of a droplet as the aftermath of entrainment from a liquid filament is shown in figure 5. Over time, a droplet instance is shown with locations connected with a line in cross-sectional (a) and longitudinal planes (b), respectively. It can be observed that droplet moves towards core and leaves behind the wavy-annular film, where origination happened. In both these figures, corresponding locations of the particle are shown after its conversion from a droplet, obtained from combined VOF-LPP method. A good match between droplet trajectory and representative particle motion justifies conversion of a droplet to a particle. Due to the conversion of a droplet to particle, the computational effort has been saved hugely. As representative savings of computational cost in VOF-LPP than VOF simulation mentioned in figure 5 is ∼24%. Speed up increases while more droplets join the cluster.

Zoom In Zoom Out Reset image size

Droplet creation through the stages of undercutting happens at different azimuthal axial locations from the film continuously, as time progresses. It makes a vast number of a droplet, which involves in own dynamics and proceeds to fill up entire gaseous core volume. In figure 5, it has been already shown that droplets move towards the core and fall back compared to film. Hence, the whole cross-section of the gaseous core fills up with many droplets, as shown in figure 6. The typical set of values for different flow situations and physical dimensions for the situations shown in figure 6 are ${v_l}$ = 0.5 m s⁻¹, ${v_g}$ = 10.0 m s⁻¹, ${\rho _l}/{\rho _g}$ = 1000, L = 6 D = 0.144 m. One can clearly understand with the progress of time (figures 6(a)–(e)) droplet swarm span increases along with population density. Variation of droplet population density at different axial planes (PP, QQ and RR) is also shown in three additional insets of figure 6. The droplet population density is found highest at a higher axial plane attached to the traversing ligament for origination.

Zoom In Zoom Out Reset image size

A longitudinal window (ABCD) is depicted in figure 6, to establish a multitude of droplet formation through the route of undercutting, as shown in figure 3. Observations of VOF based dispersion (figures 6(a)–(e)) can also be replicated using hybrid LPP-VOF simulation. This multi-scale model ensures speed up in calculations keeping the primary feature intact. A representative result of VOF-LPP simulations for the same gas–liquid inflow parameter is shown in the last row of figure 6. LPPs replace the produced dynamic droplets by the spherical liquid particles (shown in dark green in insets). Reported results for the area filled droplet at the different azimuthal axial location by using coupled VOF-LPP method have a close resemblance to VOF. An account of run time has been presented in table 2 to compare the VOF and VOF-LPP method for predicting droplet population. Initiation for this comparison has been set up as the time frame at which the first droplet is generated. Before the first droplet generation, simulation using VOF-LPP method performs proportionally to VOF method. Once droplet starts shredding, meshing adapted around the smaller droplets to resolve the features that eventually retards the VOF based simulation. On the other hand, as in VOF-LPP, smaller droplets are replaced by the particles, resulting in faster simulation. Table 2 reports ∼45% saving in simulation runtime, which involves just a short physical time (∼13 ms). The total savings cannot go proportionally; some computational time is incurred in the calculation to update more and more particles generated with time. Overall, table 2 hints the better efficiency of VOF-LPP method compared to the Eulerian VOF technique.

Figure 7 depicts the evolution process from annular to droplet annular through the route of entrainment. One can observe the number of droplets continuously increasing in the gaseous core as time progresses. Slowly and steadily whole cross-section of the tube is being occupied by droplets swarm and subsequently wall adhered liquid film depletes substantially. This figure has also shown the one-to-one correspondence between simulation made by VOF only and VOF-LPP method. Exact replication of columns in the figure substantiates the VOF-LPP predictability despite the conversion of a droplet to a particle.

Zoom In Zoom Out Reset image size

A clear identification of the similarity between the models can be also visualized from the comparison of number of droplets present in different azimuthal sectors, as predicted from VOF and VOF-LPP based schemes. Figure 7(c) depicts a good match between the predictions by VOF and VOF-LPP schemes. A comparison of total number of droplets in the domain, as visualized from cross sectional view at t* = 0.1767, is also shown in figure 7(d). For different size ranges number of droplets present in the domain is shown in the Figure. Comparison of schemes in this perspective shows a very good match which establishes the predictability of VOF-LPP based scheme as equivalent to VOF based scheme.

Furthermore, to establish the equivalence between the schemes, comparisons of axial velocities at two distinct locations have been made. The same is depicted in figure 8. Here, symbols represent VOF based results, and corresponding lines are for the VOF-LPP scheme. AA and BB are the axial planes chosen as representative with respect to the inlet at 0.8 and 2.5 times of tube diameter (D), respectively. As the momentum flux of liquid film is more than the gas phase, one can observe two peaks in the vicinity of the film interface and, after that, gradually fall and eventually become almost steady axial velocity inside the core. Both liquid and gas velocities are shown equivalence when compared between schemes. To show the similarity between the VOF and VOF-LPP based simulations, droplet count when observed from AA and BB towards the domain are shown in figure 8(b). A good match in predictions of droplet count shows the similarity between the schemes. Moreover, a total count of droplet, present in the domain when observed from axial view has been also compared between the predictions obtained from VOF and VOF-LPP schemes (figure 8(c)). One can observe that the both the schemes are having similar predictability of droplet count.

Zoom In Zoom Out Reset image size

A close observation of the droplets cluster also reveals reattachment of isolated liquid mass in wavy film after a certain period of sustainment in the gaseous core. Such phenomenon helps in increasing the film thickness and delays critical stages like dry-out involving heat transfer. In figure 9, a representative droplet shows the occurrence of deposition in the phenomenon. For clear visualization of the phenomenon, the dashed rectangular block has been used in different time frames around the event. The presence of strong vortex motion of gas trapped inside the bag like structures leads to the aberration of ligament crest (figure 9(a)). Further thinning of ligament reaches its critical limit and provides a favourable environment for 'pinch off' a liquid element (figure 9(b)). Detached liquid element, due to surface tension, molds into a spherical drop from elongated lamella shape. Isolated drop, due to neighbouring gaseous phase dynamics, starts its journey towards the wavy-annular film (figures 9(c) and (d)). After coming close to the film, droplet reunites (figures 9(e) and (f)) with film and retards film depletion rate.

Zoom In Zoom Out Reset image size

Depending on the entrainment rate (figure 3) and deposition (figure 9), film thickness varies, and the transition of annular to droplet flow occurs. Efforts have also been made to understand the reason behind the movement of the droplets by analysing the velocity vectors around it. Though there is a strong upward moving wind around the droplet, a horizontal outward velocity vector is observed in figures 9(g) and (h) to push the droplet towards the core. Such horizontal velocity vectors are generated due to the bag like lamella from which droplet has generated. The entrapment of air below the lamella orients the velocity vector from vertical to horizontal, targeting bypass of bag entrapment. In this process, the air pushes droplets towards film and results in the deposition. Depositions are only possible for the droplets originating very near to the unperturbed film region.

Close surveillance has been done to demarcate the infusion trajectory by a detached liquid element from the finite length scale ligament in the gaseous core. Deposition sequences of a droplet to the liquid film are observed for establishing similarity between VOF and coupled VOF-LPP model in figure 10. Here, the particle represented the detached small length scale liquid element and subjected to the same external forces that influence the isolated liquid element dynamics, separated from the ligament. In the cross-sectional view, (figure 10(a)), route of the drop using both VOF and VOF-LPP is almost linear and shows very close agreement. The path traced in the longitudinal view of the same drop is non-linear (figure 10(b)), but particle traced follows the same trajectory.

Zoom In Zoom Out Reset image size

Exercises reported in figures 5 and 10 successfully establish the equivalence of popular VOF and proposed VOF-LPP method in annular to droplet transformation. Shearing off of liquid element from the film generates a wide range in size and shapes of droplets. In the process of individual dynamics of these droplets, the further breakup of liquid mass is inevitable. Fast-moving air deforms larger size lamellas and eventually may fragment into multiple daughter droplets, called secondary atomization. Point by point, theoretical depiction about the bending of fluid-structure to the pinch-off of droplets from ligament has been depicted in figure 11 for a representative case.

Zoom In Zoom Out Reset image size

The entrapped air below a larger size lamella (figure 11(a)) makes efforts to bypass it, and in this process, it deforms the liquid mass (figure 11(b)). In locations, where droplets aberration becomes critical, a possibility of pinch-off has been observed. Depending on the size of the lamella, a number of droplets generated vary and even successive pinch-offs are observed. In the case of present effort, secondary atomization falls into the regime of vibrational breakup (Guildenbecher et al 2009) (0 < We < 11). As a representation, in figure 11(c), three droplets have taken birth from the lamella shown in figure 11(a). These droplets are a smaller entity and follow their dynamics (figure 11(d)). In the VOF-LPP method, provisions are kept to check the droplet volume and compare with respect to minimum cell criteria, reducing below which in terms of volume, a droplet will be converted into a particle. A proper selection of this criterion avoids converting lamella to the particle for significantly larger size liquid mass, which can further disintegrate (figure 11(c)). In figure 11(e), a large-sized lamella has not converted into the particle and supports the efficacy of the selected criterion in the present case study. It continues in grid-based structure and breaks into daughter droplets. However, as these daughters are the smaller than criterion set for a droplet to particle conversion, insertion of particles are observed in figure 11(f). Later on, these particles follow equations (8) and (9) for traversal in the domain (figure 11(g)).

Velocity vectors in the gaseous medium are observed to understand the reasons behind the disintegration of a lamella. The dynamics of the air around the liquid lamella, reported in figure 11(a), is shown in figure 12 from the top view of the conduit. It can be observed that the typical shape of the lamella has accumulated air in its bag-like shape, which thins down some section and further progress towards pinch off. The velocity vector is aligned with the lamella periphery, as shown in figures 12(a) and (b). Figure 12(c) shows the diverging nature of the velocity vector near the neck formation. On the other hand, air tries to push the daughter droplets away from each other, as shown in figure 12(d).

Zoom In Zoom Out Reset image size

Apart from the droplet disintegration, one may also observe occasional merging between two droplets to form larger sized ones in the cluster. Such merging of droplets depends on local velocity characteristics of the gas, and successful collisions are rare. As representation in figures 13(a)–(c) longitudinal view of stages of a droplet, merging are shown inside the gaseous core. Drops are represented by two distinct colors of red and blue for differentiation. Both droplets possess the unidirectional momentum towards each other with a different magnitude that allows them to move closer (figure 13(a)). As the process continued, drops reduce air film thickness in-between figure 13(b) leading towards the integration of mass figure 13(c). The unified droplet then rounds off under surface tension and follows own dynamics, even though it can break further. It is to be noted that lubrication forces between coalescing drops are not considered in present VOF formulation.

Zoom In Zoom Out Reset image size

It may deviate VOF simulation from the actual merging process. But as merging is rare in operating velocity ranges of consideration and incorporation of long-range forces is computationally expensive, the same has not targeted in present scope. The present particle-based methodology does not have provision for merging; hence, a similar phenomenon has not been observed in VOF-LPP method. But, as commented earlier, occurrences of merging situations are rare which imparts almost no error in droplet swarm prediction using VOF-LPP. Entrained mass of droplets is one crucial parameter for the study of droplet annular flow. Accumulation of mass in core signifies the reduction of film thickness and progress towards dry out. It is interesting to note the rate at which droplet is being generated from a three-dimensional wavy annular rim. Results obtained from numerical simulation are plotted in figure 14, which show the monotonous increase of mass accumulation as time progresses for a wide range of liquid to gas velocity ratio. Here, the entrained mass has been plotted as a percentage of total liquid film at the beginning of the simulation.

Zoom In Zoom Out Reset image size

A comparative look at different curves in figure 14 explains higher shearing rate for a larger ratio of water to air velocity. As shown in figure 3(a), much more protrusion formation can be observed at higher velocity ratios that filling up the cross-sectional area faster. Besides, it can also establish an increase of entrainment rate as relative velocity increases. But an increase in entrainment rate is not directly proportional to relative velocity due to the deposition mechanism of droplets whenever population becomes crowded. Though, present simulations have not been extended to dry out, a continuous increase in the mass accumulation with time hints towards depletion of overall film thickness. At later stages, in figure 14, the curve slope becomes constant; the ideas can be extrapolated to predict complete dry out time. It can be commented that dry out appears faster in a case of higher liquid to gas velocity. At the same time with$\,{v_l}/{v_g}$ = 20, a tendency of saturation in entrainment rate has been noticed which hints towards sustainment of annular regime in place of complete conversion to droplets.

The entrainment rate can also be seen faster from the cross-sectional viewpoint at the higher liquid to gas velocity. After 0.8 ms from the first droplet pinch-off, cross-sectional views are compared for three different velocity ratios and depicted in figure 15. It needs to be noted that the peak point of the disturbance zone in the annular film shift upward as time progresses due to liquid inflow (figure 6). Hence, comparisons of different axial planes are made keeping the location of this peak point (y_bend) as a reference. At three different axial location comparison are made. A critical look columnar figure establishes that droplet population density increases as the cross-sectional plane becomes nearer to the peak point.

Zoom In Zoom Out Reset image size

A careful comparison between horizontal rows establishes higher entrainment rate at a larger value of ${v_l}/{v_g}$.Throughout the present study, during the formation and breakage of a droplet, the dispersion of multiple-size droplets is observed. Bigger size droplets fragmented into smaller daughter drop while traversing in the cluster. A successful merging of droplets (noticed in VOF based simulation) has also been found out to produce larger sized droplet in a few occasions. Hence, a population balance based on a mass of droplet will be one interesting feature to be noted. We traced the droplet frequencies as a function of its characteristics diameter at the initiation and fully developed instances of the swarm for ${v_l}/{v_g}$ = 20. Same has been plotted in figures 16(a) and (b), which upon comparison shows a definite increase in the droplet population strength in the annular core. The droplet population shows the bi-modal nature of distribution in the fully developed stage, which corroborates shedding of both pinched off spherical droplets and dislodged lamellas. The first populated peak is appeared at the lower range of diameters, whereas the presence of a large number of lamellas in the annular core generates the second peak.

Zoom In Zoom Out Reset image size

Continuing in this effort, we have shown the droplet population as four equispaced intervals of diameter. The number count of droplets in these diameter intervals are shown with the progress of time in figures 17(a)–(c). As smallest size liquid droplets sustain its life cycle, avoiding further breakage, therefore, the number distribution is higher for 0.5 to1.5 mm diameter droplets. Though some higher-sized lamellas are being formed, the number density of larger size drops never shoot up due to subsequent breakage with time. It can be also seen from an increasing number of pillar heights in the sub-figures that mass of entrained droplets continuously increases and the majority of the droplets are formed at a lower size range. In figures 17(a)–(c), we have shown population density for three different velocity ratios which successfully establishes a higher entrainment rate with higher liquid inertia.

Zoom In Zoom Out Reset image size

At ${v_l}/{v_g}$ = 40 due to high liquid flux, some larger-sized lamellas are formed increasing the population density of higher-level intervals. The comparative nature of figures 17(a)–(c) depicts that entrainment explicitly depends on velocity ratio (or momentum flux ratio) of the fluids. In the domain, both symmetric and asymmetric breakage of the droplet, have been noticed. During the calculation of population balance, groups are decided based on the nominal diameter of the droplet. During breakage, if some intermediate size is generated, its frequency is distributed among neighbors with proper weight by conserving number and volume constant. Same has been followed for the merging process. In the domain with time, death and birth of a particular size group in the population is noticed. In case of fragmentation of a droplet at maximum tertiary daughters have been observed for the described velocity ranges.

Equivalence of Eulerian VOF simulation and coupled VOF-LPP prediction has been depicted in the context of wavy annular flow along with droplet swarm for air-water two-phase flow. Conversion of a smaller sized detached liquid entity into particle has shown tremendous potential in predictability and computational efficiency. The multitude of entrainment stages has given rise formation of droplets and lamellas, which contribute volume and area filling nature of dispersed liquid swarm. The identical trajectory followed by drops in VOF based simulation and particles in VOF-LPP scheme gives the confidence to predict the chaotic scenario of droplet population using a later method with higher computational efficiency. Droplets near the film inside the tube have found to rejoin back due to radial outward air flux. The breakup of bag-shaped lamella into daughter droplets is investigated from both VOF and VOF-LPP simulation. Occasional merging between droplets, which showed a very feeble influence on cluster dynamics, has also found in grid-based VOF simulation. Simulation of different gas and liquid flow rates showed faster entrainment rate in case of higher shear acting at the annular interface. A faster rate of volumetric filling of the gaseous core by droplets is observed at a high value of a velocity ratio between the liquid to gas. Study of droplet population dynamics showed a temporal increase in droplet swarm span where mainly smaller size spherical liquid entities are dominated. Lamellas fragmented into smaller droplets to keep the population density of larger entity limited. The growth of population has found to be faster in a higher ratio of liquid to gas velocity. A constant slope in mass accumulation due to entrainment will allow one to extrapolate the idea for obtaining critical stages as dry out. Finally, it can be concluded that during the prediction of multi-scale features like droplet annular flow VOF-LPP scheme is equally efficient compared to the pure VOF approach. But coupled Eulerian-Lagrangian scheme is for more efficient in terms of computational effort.

This work has been supported by the project №. 36(1)/14/08/2015-BRNS/36052 funded by the Board of Research in Nuclear Sciences, Department of Atomic Energy, India.