Published: 31 August 2016

Modeling of flow of two-phase mixture in curved channel pipeline

V. E. Lyalin1
A. N. Krasnov2
1Kalashnikov Izhevsk State Technical University, Izhevsk, Russia
2Ufa State Petroleum Technological University, Ufa, Russia
Corresponding Author:
V. E. Lyalin
Views 27
Reads 13
Downloads 1134


This article describes the three-dimensional modeling of steady motion of gas and condensed phase particles of varying severity in a curvilinear orthogonal coordinate system gas channels. A numerical method for solving the equations that describe the motion of the particles, based on the difference scheme and the application of the two-level iterative process. Built multivariate mathematical model of viscous steady flow of gas and the equilibrium of hydrate particles in axisymmetric taking into account transport and energy dissipation. Obtained pressure field in the pipe by solving the model equations using SIMPLE method. Numerical calculations of the velocity field of the gas and the dispersed phase, the trajectory of the mass flow deposited on particles of different sizes with the wall during the wet gas in the curved pipe. The calculation of particle trajectories in different kinds of process equipment show at the option of either abrasive wear parts located in a stream, or an intense build-up and accumulation of hydrates. The results can be applied in the analysis of impact of solid particles of various sizes on the elements of gas pipelines.

1. Introduction

According to studies hydrocarbon hydrates are white crystalline solids, like snow, and at seal similar to ice [1]. The condensed phase in the flow of natural gas hydrates in addition may contain other solid impurities (slag, sand, etc.) Mechanical impact at high velocity a mixture of gas and particles can be exposed structural elements and stop measuring gas equipment. Assessment of the impact can be made on the basis of solving the equations of motion of a two-phase mixture in the pipeline elements.

2. The equations of two-phase flow of gas in the curved pipes

The system of equations describing the steady flow of viscous heat-conducting gas in an arbitrary coordinate system is as follows:

ρV=0, ρVV+P=0, ρ(eV)+P(V)+q=0,

where ρ – density of the gas; V – The velocity vector; e – Specific internal energy; P – Pressure tensor related to pressure p and viscous stress tensor S relation P=pI+S. The viscous stress tensor is defined as:


where I – the unit tensor, μ – dynamic viscosity coefficient. Tensor components S are shown in many literatures, for example in [2]. Heat flow q is determined by the Fourier law: q=-λT, where λ – coefficient of thermal conductivity, T – temperature. For the system of Eq. (2), you must add the equation of state of the gas in the form p=ρRT, where R – the gas constant. In the calculation of flows with low speeds (up to Mach numbers M<0.3), it is advisable to consider for incompressible. In this case, in the Eq. (1) ρ=const and the equation for the energy can not be seen due to the small change in the flow temperature.

Together with the Eq. (1) it is necessary to solve the equation for the turbulent viscosity. It was used model of turbulent viscosity [3].

The equation for the turbulent viscosity is as follows:


where νT – turbulent viscosity; μ=ρ(νT+ν0), ΓT – A term describing the generation and dissipation of turbulent viscosity [3]; ν0 – Molecular viscosity.

To calculate the flow in curved pipes need to solve the problem of flow in a three-dimensional setting. Solving such problems is difficult enough. To evaluate the mechanical interaction of a condensed phase with the walls will be considered two-dimensional setting. The system of Eqs. (1), (2) can be written in a flat or cylindrical coordinates with axial symmetry (x – longitudinal, r – transverse coordinates). For the curved regions of the flow plane x, r applied coordinate transformation taking curvilinear computational domain into a rectangle. The difference mesh is constructed using a complex boundary element method [4].

To write equations used curvilinear orthogonal coordinate ξ, η. System variables ξ, η associated with Cartesian coordinates x, r as follows: ξ=ξ(x,r), η=η(x,r). The components of the metric tensor are equal g11=g22=D. In the case of orthogonal coordinates xξ=rη, xη=-rξ, D=xξ2+rξ2.

The equations of continuity, momentum and energy in the coordinate system ξ, η are of the form [5]:


where J=D; u, v – Projections of the velocity vector v on the axis ξ, η; E=ρe+0,5V2 – the total energy; T – temperature. Axisymmetric flow corresponds to γ=1, a flat – γ=0. The coefficients in the right sides of the equations in the coordinate system ξ, η is written as follows:

P11=μ2J-1(uξ+J-1Jηv)-23divV, P22=μ2J-1(vη+J-1Jξu)-23divV,
P33=μ2r-γ(J-1rξu+J-1rηv)-23divV, P12=μD-1Juη+Juξ-2uJη+vJξ,

Odds Q1, Q2, Q3 describe the exchange of momentum and energy between phases:

Q2=nBuu-us, Q3=nBuv-vs,

where, the index s belongs to the dispersed phase; G – The mass evaporation rate drops; cs, ρl – heat capacity and density of droplets of substance; n – Condensed number of particles per unit volume; Bu, BT – Heat transfer coefficients and the resistance between the phases.

Together with the Eq. (3) it is necessary to solve the equation for the turbulent viscosity. In the plane x, r orthogonal difference grid built complex boundary element method.

The basis of the solution of a system of Eq. (3) laid SIMPLE method [6], which is kind of focused on the use of curvilinear coordinates, implemented in [5, 7].

The equations for the dispersed phase is also recorded in a curvilinear coordinate system:

(JrγρsusEs)ξ+(JrγρsvsEs)η=J2rγQ1, (Jrγnsus)ξ+(Jrγnsvs)η=0.

Resistance and heat transfer coefficients are determined by the formulas [6]:

Bu=9μ2ρlrs2f, f=1+0.25Re+0.0117Re, Re=v-vsrsρμ,
BT=BuNucp3Prcsf, Nu=2+0,459Re 0,55Pr0,33.

3. The boundary conditions

On the permeable boundaries of the fixed costs of the gas and condensed phase, the enthalpy of the phase and the initial particle size. On impermeable boundaries puts a condition of sticking. At the outlet of the pipe for subsonic and where the pressure was set. On the axis of symmetry – the symmetry condition.

Changing the number n of particles per unit volume is determined by the speed lag of gas (including the interaction of the dispersed phase and gas). We write the equations of momentum and energy transfer:


where ςs – coordinate in the direction from the sth path; Ws=(Us)2+(Vs)2; Us=usξx+vsξr; Vs=usηx+vsηr; fs=(us,vs,Ts); Cs=A1sA1sA1T; Bs=A1su+axgA1sv+aygA1TsT .

Driving integration of the system Eq. (5) is shown in Fig. 1.

The integration is performed over the interval of the particle trajectory. The segment begins with the diagonal cell comprising an anchor point i,j, at 0 particles with known parameters, and ends at the node (i,j), where the parameters to be calculated. To integrate the use of implicit Runge-Kutt second order of approximation:


An index of 0 corresponds to the beginning of the segment of the trajectory. Required parameters φs0 at the point 0 determined by the formula:


using known values φs at the ends of the diagonal. The values αs, Δςs, iu, jv are given by:

αs,i+iu,j=1ΔξiUWs,i,j+UWs,0, αs,i,j+jv=1ΔηjVWs,i,j+VWs,0,
Δςs=2αs,i+iu,j+αs,i,j+jv, Δξi=ξi-ξi+iu, Δηj=ηj-ηj+jv,
iu=sign(-Us,i,j), jv=sign(-Vs,i,j).

The equation for ns(ξ,η) Eq. (4) is solved by countercurrent scheme [6]. The difference analog of the Eq. (4) has the form:


where aWi,j=max[δηjyi-1/2,jDi-1/2,jUi-1/2,j,0], aEi,j=max[-δηjyi+1/2,jDi+1/2,jUi+1/2,j,0], aSi,j=max[δξiyi,j-1/2Di,j-1/2Vi,j-1/2,0], aNi,j=max[-δξiyi,j+1/2Di,j+1/2Vi,j+1/2,0],


The boundary conditions for the system Eq. (6).

Since the Eq. (6) recorded along the paths, which are characteristics that only the boundary conditions are placed on parts of the border, where the particles start their movement. Here you can set the initial parameters of the particles: u0s, v0s – the initial velocity of the particles; T0s – The initial temperature of the particles.

For the Eq. (4) applies, counterflow circuit, so the boundary conditions are set similarly.

Thus, to solve the equations that describe the motion of particles, obtained difference scheme that monitors the direction of flow. To solve differential equations has the following two-level iterative process.

As an initial approximation for the particle velocities are taken gas velocity. Then, depending on the course of nature, there are two variants of the algorithm for solving this problem.

1st option. Within this field, there is no areas with a return (circulation) over. In this case, we realized marching along the longitudinal coordinate algorithm. Iterations on points held on the transverse coordinate. As the first value Bsi,j, Csi,j is taken Bs0, Cs0. Then iterate to refine the position of the point 0 on the diagonal of each cell and solved the Eq. (2-6). After that solved the Eq. (5) under the scheme Eq. (6) at each point i, j and specified values Bsi,j, Csi,j. In this embodiment, the parameters in the computer memory is stored on the particle layer on the layer i and i- 1.

2nd option. The computational domain contains zones with reverse flow. In this case, point by point iteration are conducted throughout the computational domain. computer memory requirements in this embodiment is considerably higher than in the first embodiment. In order to achieve the residual from the equation of continuity ~10-4 takes only about 15 global iterations.

The presented numerical method allows the calculation of the velocity field of the gas and the dispersed phase, the trajectories of particle motion (by solving the equations dη/dξ=Vs/Us when found Us, Vs), the mass flow of particles deposited on the walls G.

Fig. 1A scheme for integrating the equations of motion of particles

A scheme for integrating the equations of motion of particles

4. The current results of the calculations in complex configuration channels

For the numerical solution of the problem of motion of the gas and particles at the site of bending the pipe [4], a difference mesh, as shown in Fig. 2. Lines ξ=const are longitudinal, a line η=const orthogonal thereto.

The movement of the gas phase was calculated at a Reynolds number 50000. The vector velocity field is shown in Fig. 3.

As follows from the results of calculations, the flow field is not symmetric. First, during urged against the wall of a small radius, and a loan in the opposite direction. As a result, the flow deflection from the center of the eddy currents are formed, as shown in Fig. 3. The velocity profile is characteristic of a turbulent profile.

Particle motion is calculated for an equivalent particle diameter from 10 microns to 500 microns. The trajectories of particles of different sizes are shown in Figs. 4-7.

Particles of small size (10 microns) monitor the gaseous phase line current and the walls of virtually interacting. Larger particles 100 microns in diameter collide with the pipe wall after bending.

Large particles (300 and 500mkm) have a more direct path to the bending portion of the tube and collide with the wall with a high intensity. Low speed 300 micron particles can be trapped vortex gas flow and to be in it for a while. Large particles with a diameter of 500 mm, almost all come through the vortex. Only particles from the wall region with a very low speed to make a whirlwind return movement.

Fig. 2Curved orthogonal net curved pipe

Curved orthogonal net curved pipe

Fig. 3The velocity field in the curved section of the pipe

The velocity field in the curved section of the pipe

Fig. 4The particle trajectories diameter 10 mm

The particle trajectories diameter 10 mm

Fig. 5The particle trajectories a diameter of 100 microns

The particle trajectories  a diameter of 100 microns

Fig. 6The particle trajectories diameter of 300 microns

The particle trajectories diameter of 300 microns

Fig. 7The particle trajectories diameter of 500 microns

The particle trajectories diameter of 500 microns

As a result of solving the problem of the motion of a two-phase mixture of particles and gas through the curved portion of the pipe calculated mass flows falling on the wall of particles of different sizes. The results are shown in Fig. 8.

Fig. 8 shows the distribution of massive flows of precipitating particles G in several sizes. Coordinate S measured at the tube wall. Mass flow rate G is related to the mass flow rate of particles on a straight section of pipe (ρsWs). These data can be used to assess mechanical impact on the particles pipeline wall at different degrees of drying of natural gas.

Narrowing the flow cross section can be various kinds of control and measurement devices. The computational domain and the difference curvilinear grid shown in Fig. 9.

Fig. 8Distribution of mass flows deposited on particles of the wall of the pipe

Distribution of mass flows deposited on particles of the wall of the pipe

Fig. 9The difference in net contraction

The difference in net contraction

Fig. 10The flow field in the constriction

The flow field in the constriction

Fig. 11The particle trajectories diameter 10 mm

The particle trajectories  diameter 10 mm

Fig. 12The particle trajectories diameter of 500 microns

The particle trajectories diameter  of 500 microns

Fig. 13Mass flow deposited particulate

Mass flow deposited particulate

Calculation of the gas phase of motion parameters. The results are illustrated in Fig. 10.

Before narrowing of the pipe cross-section is observed stagnation zone with poor circulation flow in the wall region. After opening for jet nature is the presence of extensive areas (up to several pipe diameters) of the reverse flow. The movement of the condensed particles is shown in Fig. 11 and 12.

Small particles in this device follow the gas and do not fall into an obstacle.

Coarse particles have an intense effect on the structural member. precipitation intensity distribution for frequent diameter of 100-500 microns is shown in Fig. 13.

Here, the coordinate of S corresponds to a distance from the pipe wall. Particles larger than 100 microns fall almost entire surface with obstacles fairly high intensity.

5. Conclusions

Figs. 6-8 demonstrate the possibility of increased erosion ash pipelines in areas of bending and rotation. The calculation of particle trajectories in the technological equipment of various kinds of shows at the option of either abrasive wear parts located in a stream, or an intense build-up and accumulation of hydrates.

The calculated results can be applied in the analysis of impact of solid particles of various sizes on the elements of gas pipelines.


  • Strizhov I. N., Hodanovich I. E. Gas Production. Institute of Computer Science, Moscow-Izhevsk, 2003, p. 376, (in Russian).
  • Oran Je., Boris Dzh. Numerical Simulation of Reacting Flows. Second Edition Cambridge University Press, 2001, p. 640.
  • Sekundov A. N. The use of differential equations for the turbulent viscosity and analysis of non-self plane currents. Bulletin of the Academy of Sciences USSR, Vol. 5, 1971, p. 114-127, (in Russian).
  • Gromadka II T., Lej Ch The Complex Boundary Element Method in Applied Sciences. Mir, Moscow, 1990, p. 303, (in Russian).
  • Benderskij B. Ja, Tenenev V. A. Experimental and numerical study of flows in axisymmetric channels of complex shape with blowing. Bulletin of the Russian Academy of Sciences, Vol. 2, 2001, p. 24-28, (in Russian).
  • Patankar S. Numerical Methods for Solving Problems of Heat Transfer and Fluid Dynamics. Jenergoizdat, Moscow, 1984, p. 150, (in Russian).
  • Benderskij B. Ja, Tenenev V. A. Spatial subsonic flows in areas with complex geometry. Math Modeling, Vol. 13, Issue 8, 2001, p. 35-39, (in Russian).
  • Amimul Ahsan Two Phase Flow, Phase Change and Numerical Modeling. InTech, 2011, p. 596.

About this article

05 August 2016
08 August 2016
31 August 2016
Flow induced structural vibrations
model of gas flow
curved pipe channels
dispersed gas particles