A New Analytical Procedure to solve Two phase Flow in Tubes

A new formulation for a proposed solution to the 3D Navier-Stokes Equations in cylindrical 1 co-ordinates coupled to the continuity and level set convection equation is presented in terms of an 2 additive solution of the three principle directions in the radial, azimuthal and z directions of flow and 3 a connection between the level set function and composite velocity vector for the additive solution is 4 shown. For the case of a vertical tube configuration with small inclination angle, results are obtained 5 for the level set function defining the interface in both the radial and azimuthal directions. It is found 6 that the curvature dependent part of the problem alone induces sinusoidal azimuthal interfacial 7 waves wheras when the curvature is small oscillating radial interfacial waves occur. The implications 8 of two extremes indicate the importance of looking at the full equations including curvature. 9


Introduction
Annular two phase flow is frequently experienced in industrial applications, for example, in power generation plants and heat removal devices.Applications include transfer lines in pipes (gas-liquid oil wells), evaporators and condensers.This flow regime involves a liquid phase flow in the form of a thin film along the pipe wall and a core region of a gas phase where there are droplets entrained within the gas.An important feature of annular pipe flows is the formation of waves at the core/film interface.Thus the Navier Stokes equations which describe the flow problem must be addressed inorder to determine if the 3-D velocity components exhibit wave-like characteristics.The Navier Stokes equations have been dealt with extensively in the literature for both analytical [1,2] and numerical solutions [3,4].The level set method, has been used originally as a numerical technique for tracking interfaces and shapes [5,6] and has been increasingly applied to various areas of engineering and applied mathematics.In the level set method, contours or surfaces are represented as the zero level set of a higher dimensional function called a level set function.This can be the distance from the particular phase of material to the interface.For example in fracture mechanics level set methods have been used to track the shape around a crack in two and three dimensions that is propagating with a sharp kink [7].In addition, various applications in image segmentation have been used with corresponding active curve evolution algorithms [8,9].Reachability analysis is frequently used to study the safety of control systems.Using exact reachability operators for nonlinear hybrid systems is presented in [10].An algorithm for determining reachable sets and synthesizing control laws is implemented using level set methods in [10].Various models used to compute the interaction of 3D incompressible fluids with elastic membranes or bodies, rely on the use of level set functions [11], to capture the fluid-solid interfaces and to measure elastic stresses that have been used.In [11] the computation of equilibrium shapes of biological vesicles is presented and numerical simulations of spontaneous cardiomoyocyte contractions is presented.A conservative method of level set type for moving interfaces in divergence free velocity fields is presented in [12,13].The method in [13] was coupled to a Navier-Stokes solver for incompressible two phase flow with surface tension.Wave phenomena is known to exist at the interface of two phase immiscible flows [14].In the present paper, we present a level set method for moving interfaces for such velocity fields which are coupled to Navier-Stokes equations for two phase flows in tubes.The need for a cylindrical co-ordinate system is apparent as most of the industrial applications involve assemblies of round pipes and therefore a simple model using a straight tube is necessary to obtain insight into the more complex problem found there .The novelty of the present work is to reveal an analytical approach in solving the 3D cylindrical Navier-Stokes equations where the three principle directions of flow, in radial, azimuthal and longitudinal directions are summed to form a new composite vector velocity expression.In this light, we propose to solve a curvature only formulation of the governing equation for the level set function and one in which curvature is added to the governing equation in the large time evolution for the level set function and hence composite velocity formulation.The importance of these two extremes is to get some insight into the nature of and existence of interfacial waves and level set function for the interface and velocity profiles in annular flow.Also it is hoped that the method outlined in this work will lead to eventually solve the full problem analytically for the composite velocity with full curvature expression in cylindrical coordinates included.To the best of our findings this remains an open problem on the subject of two phase flows.Finally, the method is new in the sense that having an analytical solution allows us to solve for velocity and coupled temperature fields easier than employing numerical methods for certain cases.Our results for the limiting cases demonstrate how the velocity and level set function exhibits wave like phenomea in both the radial and azimuthal directions of flow.The problem is solved by assuming the level set function is a product of a decreasing exponential term with a function of radial and azimuthal coordinates.It is shown that when the rate of decrease of time dependent term is very high the level set function exhibits oscillations with higher frequency and lower amplitude and eventually approaches a hyperbolic tan function which is a steady state level set function.Further results showing composite velocity are also shown.

Level Sets in Cylindrical Co-Ordinates
Let φ be a level set function [5,15].The gradient of the level set function in cylindrical co-ordinates is defined as: The mean curvature, κ, of the interface defined by the zero isocontour of the level set function φ [5,15], is the divergence of the normal to the interface given by Thus it can be expressed as: The mean curvature κ of a dynamic surface φ(r, θ, z, t) = ψ(r, θ, t) + z in cylindrical co-ordinates is, where The geometric trace of a dynamic closed surface that bounds an open set can be represented implicitly as and for two phase flow has two separate regions where φ > 0 and φ < 0 respectively [15].The surface evolution is determined by: In this work the function φ is what is called a filter which has an effective width and will be defined on a small region near the interface.The form of the Navier-Stokes Equations that are not coupled to the interface of the two phases are: where for cylindrical (r, θ, z) coordinate system, Laplace operator has the form the gradient is given by Equation ( 1) and u is the three dimensional velocity of the flow field and is described further below.

A New Composite Velocity Formulation
The 3D cylindrical incompressible unsteady Navier-Stokes equations coupled to interface convection equation are written in expanded form, for each component, u r , u θ and u z , where the remaining non linear terms appearing in full Navier-Stokes equations are suppressed for the time being: and where u r is the radial component of velocity, u θ is the azimuthal component and u z is the component along the pipe in the direction of fully developed upward flow, ρ is density, µ is dynamic viscosity, Fg r , Fg θ , Fg z are body forces on fluid and σ is surface interfacial tension.Introducing the following definition, Multiplying Equations ( 8)-( 10) by unit vectors e r , e θ and k respectively and adding Equations ( 8)- (10) gives the following equation, for L = u r e r + u θ e θ + u z k and L, The level set function φ is governed by Equation ( 7), which when expanded becomes in cylindrical co-ordinates: The continuity equation in cylindrical co-ordinates is Multiply Equation ( 11) by Multiply Level set function convection Equation ( 12) by ρ µ L: Adding Equations ( 14) and (15) gives By product rule we rewrite the previous equation as: It is noted that since L is a composite velocity term it must have units of length per time, and since φ is usually taken to be a distance function from the interface we can consider the following expression in terms of φ : where µ ρ has SI units of m 2 /s and L is defined previously and is of dimension 1/cm and φ is dimensionless, (See Appendix A for decomposition of µ ρ ) where ρ 1 , µ 1 are density and dynamic viscosity of fluid 1 and ρ 2 , µ 2 are density and dynamic viscosity of fluid 2 separated by interface and we write It is assumed that the density and viscosity of fluid 1 is negligible compared to fluid 2. From Level set Equation (7), From the first line of Equation ( 17) using Equation ( 18) we have a derivative term in t, first we write: and for ρ 1 relatively small in comparison to ρ 2 and for µ 1 relatively small to µ 2 we obtain and for derivative terms in r, θ, z , for ρ 1 small relative to ρ 2 and for µ 1 small relative to µ 2 from the first line of Equation ( 17) it can be proven that,

Special Case Solution
In this section, it is assumed that the tube is in a vertical configuration with a small inclination angle, with assumption that ρ 2 and µ 2 are nonzero terms for fluid 2 whereas fluid 1 has approximately zero density and viscosity.This is an idealization of the full problem which is amenable only to numerical treatment.It is worthy to notice that part of Equation ( 22) can be rewritten using the time dependent continuity equation Equation (13), Use of Equation ( 17), Equations ( 18)-( 24) gives a non-linear PDE Using Equation ( 4) and incorporating the curvature κ, and using Equation ( 18) and dot product in Equation ( 19), with the assumption of small inclination angle, Equation (25) for fully developed flow in the vertical z direction becomes, Fz is force of gravity in inclined tube.Incorporating Equation (18) and partial fraction decomposition for µ ρ found in Appendix A, Equation(26) becomes for ρ 1 and µ 1 negligibly small in comparison to ρ 2 and µ 2 respectively of fluid 2, where we have solved Equation (26) using the assumption of small inclination of tube, and have subsitituted L = ∇φ and then set φ = e αt F(r, θ) for α < 0 in Equation(26).This leads to Equation (27).
It can be proven that F(r, θ) is multiplicatively separable when term in t is dropped or t gets large since exponential term is in the denominator for α < 0 , thus, F(r, θ) = f (r)g(θ) and we obtain the following, It can be shown that g(θ) is an arbitrary function of θ, for which one chooses g(θ) to be constant and large inorder to be consistent with the omission of the term in e αt in Equation ( 27).

Consideration of Curvature Alone
Secondly we solve the curvature pde given in Equation ( 4).For small inclination angle of tube Equation (4) reduces to, where c 3 is sufficiently large.For constant c 2 negative there is a sinusoidal component of the azimuthal part of φ(r, θ, t) = f 1(r) f 2(θ) f 3(t).

Results and Discussion
An analytical solution is found to exist for Equation (28) and is given by the following where

Conclusions
It is worth noting that the interfacial oscillations occurring as extremes of two problems one for curvature coupled to composite velocity and the other for curvature alone presents the daunting problem of solving the full equation of Equation (25) without the assumption of small inclination angle.There is a plethora of results for computational multiphase flow using level set methods.The advantage of the present work lies in that analytical results are possible for the two extreme cases presented.It is conjectured at this point that the combination or full Equation (25) without the small inclination angle, i.e., approaching a horizontal tube configuration of flow, is non separable due to the inherent complexity of Equations ( 4) and (25) combined.We can expect that there be a very complex relationship between azimuthal and radial components of L. Work on the complete problem for this complex relationship is in progress for future studies.

Figure 2 .
Figure 2. f versus r for α = −1000 in level set function φ, R i = 0.85 cm is the radial value at the gas/liquid interface.

Figure 3 .
Figure 3. f versus r for α = −5000 in level set function φ, R i = 0.85 cm is the radial value at the gas/liquid interface.

Figure 4 .
Figure 4. f versus r for α = −10,000 in level set function φ, R i = 0.85 cm is the radial value at the gas/liquid interface.

Figure 5 .
Figure 5. f versus r for α = −100,000 in level set function φ, R i = 0.85 cm is the radial value at the gas/liquid interface.

Figure 6 .
Figure 6.f versus r for α = −1,000,000 in level set function φ, R i = 0.85 cm is the radial value at the gas/liquid interface.