The meshless numerical simulation of Kelvin–Helmholtz instability during the wave growth of liquid–liquid slug flow

https://doi.org/10.1016/j.camwa.2020.08.006Get rights and content

Highlights

  • Kelvin–Helmholtz instability (KHI) was simulated using meshless radial basis functions.

  • Fractional step method was used to solve the modified Navier–Stokes equations.

  • The numerical procedure can handle KHI problems under the various of fundamental parameters.

  • The density ratio and initial velocity difference play the important roles on the growth of KHI.

Abstract

Kelvin–Helmholtz instability (KHI) occurs at the interface of two fluids, in which the heavier fluid flows at the bottom. In the present numerical work, the KHI was analyzed by solving the modification of two-dimensional (2-D) incompressible Navier–Stokes equations. To capture the interface, the Cahn–Hilliard equation was implemented into the Navier–Stokes equations. The phenomena around the KHI of two and three-component fluids were investigated numerically by using radial basis function (RBF) combined with the domain decomposition method (DDM) in a primitive variable formulation. Here DDM is able to solve the large scale problem. On the other hand the calculation accuracy decreases with the increase of the number of the subdomains. For the above reason, in the present works, the domain was partitioned into 15 x 15 subdomains in order to reduce the decrease in accuracy due to the domain division. Next, fractional step method was used to solve the modification of the Navier–Stokes equations.

The numerical results indicate that the procedure used in the present work can easily handle the KHI problem under the variations of the interface thickness, density ratios, and initial velocity differences. Moreover, the interface evolutions obtained from the present method agree well with those of the finite difference method. The effects of the interface thickness, density ratio and magnitude of velocity difference on the KHI were also investigated. The decrease of the interface thickness produces a non-smooth concentration profile, and the increase of the interface thickness produces too much surface diffusion. It was found also that the increase of the density ratio reduces the growth of KHI. The interface rolls up are strongly affected by the initial horizontal velocity difference. Finally, the present study also shows that the RBF method is a reliable method to solve the KHI on the complex domains.

Introduction

Slug flows are often encountered in the practical applications during the hydrocarbons transportation pipeline system [1]. In general, it is intermittent and is formed by sequences of large gas bubbles followed by the liquid slug. It flows under the high frequencies in which large gas bubbles flow alternately with liquid slug [2], [3]. In the horizontal pipes, it is formed from a stratified flow due to the growth of hydrodynamic instabilities. Here small perturbations in the form of small waves naturally emerge. Those waves continue to grow, capturing the liquid that flows in front of them until the pipe cross-sectional area is fully blocked by the liquid, thus forming the slugs [4]. This wave growth mechanism is known as the Kelvin–Helmholtz instability (KHI) [5], [6].

Slug flow is one of the dangerous flow patterns in the petroleum industry. In order to ensure the safe operation of the pipeline, the knowledge of the slug flow mechanism is essential. The appearance of slug flow in the pipeline will increase the structural disadvantage that leads to crack and corrosion of the pipe [7], whereas high-frequency slug flows increase the rate of the pipeline corrosion [8]. Consequently, it will decrease the hydrocarbon production and increase the equipment maintenance cost [9]. From the viewpoint of multiphase-flow, the understanding of KHI is very important in pipeline design.

The researchers have developed the numerical methods to simulate the KHI in horizontal pipes or horizontal channel. Simmons et al. [10] explored the influence of KHI during the transition to slug flow in horizontal pipes. The basic mechanism of KHI is explained as the existence of the uniform shear stress between two different fluids. To simplify the calculation, the surface tension and viscosity were neglected, and, consequently, KHI is considered only as the fluid instability of a lighter fluid is placed on the top of a heavier one. San & Kara [11] calculated 2-D KHI in order to examine the performance of six Riemann flux formulas. Those formulas are the combination of the Rusanov scheme, Reo scheme, Harten–Lax–Van Leer scheme, first-order centered scheme, advection upstream splitting method, and Marquina scheme. Here, a single mode perturbation was used as a tested case. They concluded that the advection upstream splitting scheme results in the smallest numerical dissipation.

Patnaik et al. [12] developed a numerical simulation for KHI between two stratified fluids by using the finite difference method. They concluded that the minimum Richardson number required for the onset of KHI was at Reynolds number of 100. They claimed that their result is in a good agreement with that of the linear stability theory. On the other hand, Balabel et al. [13] introduced the level set method to solve the hydrodynamic instability. The numerical procedure was developed on the basis of the solution of Navier–Stokes equations and structural 2-D computational grid by using the control volume approach. The results showed that the onset of roll-up mechanism at the interface is largely affected by the considered flow regimes. Next, Lee et al. [14] examined the finite difference method to solve the KHI of multi-component fluids. It was modeled by both of modification of Navier–Stokes and the multi-component Cahn–Hilliard equations. They concluded that the linear growth rate for the KHI of two-component fluids was in good agreement with the analytical results. These finite difference and finite volume methods used the mesh to solve problems, and it can produce lot of numerical problems during the discretizing of complex domains with irregular boundaries [15]. Moreover, those methods were developed on the basis of the polynomial interpolations, whereas they are limited by the algebraic convergence rates [16].

Shadloo et al. [17] used the smoothed particle hydrodynamics (SPH) scheme to simulate the 2-D KHI problem of an incompressible two-phase immiscible fluid. They concluded that the numerical algorithm is able to handle the two-phase fluid flow under the variation of the density ratios. Yue et al. [18] studied also the growth of the KHI at the interface between two ideal gases using the SPH scheme. The linear growth rate in the numerical solution was found to be in a satisfactory agreement with that of the analytical solution, in which the average relative error was 13%. As the SPH is an approximate method, it is noticed that there are some disadvantages if it is not applied correctly in a specific problem. First, it is very sensitive to the number of particles used in the spatial approximation. Second, the precision of the SPH depends on the order of accuracy of the window function [19].

The above facts indicate that there is a strong disagreement from the previous works regarding the best solution to describe the physical phenomena around KHI. For this reason, a systematic numerical study to explain the involved phenomena around the transition to slug flow of liquid–liquid two-phase flow is needed. The aim of the present study is to carry out a numerical investigation of the KHI during the transition to slug flow. In the present numerical study, the radial basis function (RBF) method is used to simulate the phenomena using modification of the Navier–Stokes equations in terms of primitive variables and the Cahn–Hilliard equations in order to capture the evolution of the interface between two fluids.

An attractive and simple approach that is suitable for a wide range of the problems employs the capability of radial basis function (RBF) [20]. It has a good capability to solve the governing equations in the strong form, without any integration is required. In addition, it is also effective to solve the numerical terms with a low numerical diffusion [21]. The main challenge to implement the RBF method is the applicability to handle the large-scale domain. RBF method should be combined with the DDM, when handling large-scale domain. The basic idea of DDM is to divide one large domain into a number of smaller subdomains [22]. The DDM is able to handle large-scale problems, meanwhile it will decrease the calculation accuracy if more subdomains are employed [23].

Several numerical experiments were carried-out in order to examine the accuracy of RBF method combined with DDM. Mai-Duy et al. [22] used it to solve one and two-dimensional elliptic problems. The domain was divided into sixteen subdomains, and the results show a high level of accuracy and a fast rate of convergence. Therefore it presents more attractive than the one-domain method to solve the elliptic problems under a large computational domain. Zhou et al. [23] developed RBF method combine with DDM to solve the elliptic partial differential equations. They investigated the accuracy of the numerical solution during the increase of the numbers of subdomains. In case of 16 × 16 subdomains, they found that the root-mean square error (RMSE) was 0.005972, which represents a high accuracy of the numerical method.

In general, the RBF and SPH methods have several similarities. Those are as follows: (1) they are simple to implement, (2) they are a truly meshless method, and (3) they were developed on the basis of strong form of partial differential equations. Moreover, the RBF method has a better accuracy than the SPH method at the expense of more computational effort [24]. The need to invert the coefficient matrix makes the RBF more expensive than SPH method. However, the efficiency of the RBF method can be improved by using the domain decomposition method.

In the present numerical works, the combination of RBF and DDM was implemented, because it is a very efficient method and high accuracy with low numerical diffusion as suggested by [15] and [25]. The DDM is used to solve the multi small subdomain problems instead of single large domain problem. Next, the Cahn–Hilliard equation is used to capture the evolution of the interface, because it is a thermodynamically consistent equation in order to follow the interface evolution as reported by [26] and [27].

Section snippets

Radial Basis Function (RBF) method

The main advantage of RBF approximation is the ability to work for high dimensions of the arbitrary geometries. RBF is a function whose value depends only on the distance from some center point. By using the distance functions, it can be easily implemented to reconstruct a plane or surface using scattered data in two- and three-dimensional (2-D and 3-D) or higher dimensional spaces. It was initially developed for multivariate data and function interpolation. Hence, the meshfree nature has

Problem definition

For the simulation of KHI, it is considered that the 2-D motion of two immiscible fluids of different densities and velocities separated by an interface, in which the heavier fluid is at the bottom. In the numerical analysis, the Cartesian coordinates were selected. For simplicity, the x-dimension of the computational domain L (0 < x < L) is chosen to be equal to the domain height H (L = H). The fluids move counter-currently at the uniform velocity U as shown in Fig. 2. In the calculation,

Results and discussions

In the next section, the meshless method on the basis of the RBF in the framework of the fractional step method is used to solve the KHI of two incompressible and immiscible fluids. The instability is described by the modification of 2-D incompressible Navier–Stokes equations as shown in Eqs. (31)–(32). Here the Cahn–Hilliard equation is used for capturing the interface of working fluids.

In the present numerical study, the point independence tests were carried-out firstly in order to determine

Conclusions

The interface evolution of the KHI in a rectangular cavity and complex domains was investigated numerically by using RBF and domain decomposition method. The phenomena around KHI were modeled by the modification of Navier–Stokes equations. Here the convective Cahn–Hilliard​ equations were applied to capture the interface between the fluids. The numerical method on the basis of the implicit Euler scheme for temporal discretization and RBF method for spatial discretization were also utilized. The

CRediT authorship contribution statement

Eko Prasetya Budiana: Conceptualization, Methodology, Investigation, Programming, Writing - original draft, Formal analysis. Pranowo: Investigation, Data checking, Formal analysis. Indarto: Resources, Writing - review & editing, Funding acquisition, Supervision. Deendarlianto: Analysis validation, Writing - review & editing, Resources, Funding acquisition, Supervision.

Acknowledgments

The work was carried out and partially funded by the Indonesian Ministry of Research and Technology/National Agency for Research and Innovation , and Indonesian Ministry of Education and Culture, under World Class University Program managed by Institut Teknologi Bandung. Next it was funded partially also within research programs “World Class Researcher 2019” and the Ph.D. scholarship for Eko Prasetya Budiana (first author) from the Directorate General of Higher Education, Ministry of Research,

References (36)

  • LingL. et al.

    Preconditioning for radial basis functions with domain decomposition methods

    Math. Comput. Modelling

    (2004)
  • DuklerA.E. et al.

    Gas-liquid slug flow-knots and loose ends

  • FabreJ. et al.

    Modeling of two-phase slug flow

    Annu. Rev. Fluid Mech.

    (1992)
  • AnsariM.R.

    Dynamical behavior of slug initiation generated by short waves in two-phase air-water stratified flow

    ASME Heat Transfer Div.

    (1998)
  • FanZ. et al.

    Initiation of slugs in horizontal gas-liquid flows

    AIChE J.

    (1993)
  • SanO. et al.

    Evaluation of Riemann flux solvers for WENO reconstruction schemes: Kelvin–Helmholtz instability

    Comput. & Fluids

    (2014)
  • PatnaikP.C. et al.

    A numerical solution of Kelvin-Helmholtz waves of finite amplitude

    J. Fluid Mech.

    (1976)
  • BalabelA.

    Application of level set method for simulating Kelvin-Helmholtz instability including viscous and surface tension effects

    Int. J. Math. Comput. Sci.

    (2015)
  • View full text