- Email: [email protected]
Correcting interface turbulence viscosity using CFD modeling for predicting stratified gas–liquid flow shear stress in horizontal pipes
European Journal of Mechanics / B Fluids 79 (2020) 202–211
Contents lists available at ScienceDirect
European Journal of Mechanics / B Fluids journal homepage: www.elsevier.com/locate/ejmflu
Correcting interface turbulence viscosity using CFD modeling for predicting stratified gas–liquid flow shear stress in horizontal pipes ∗
Ping Zheng a , , Liang Zhao b a b
College of Petroleum Engineering, Liaoning Shihua University, Fushun, Liaoning 113001, China College of Civil Aviation Safety Engineering, Civil Aviation Flight University of China, Guanghan, Sichuan 618307, China
article
info
Article history: Received 30 June 2017 Received in revised form 6 September 2019 Accepted 17 September 2019 Available online 19 September 2019 Keywords: Stratified gas–liquid flow Turbulence viscosity Liquid holdup Pressure drop Shear stress CFD modeling
a b s t r a c t The prediction of different shear stresses is one of the great challenges of turbulent–turbulent stratified two-phase flow in horizontal pipes. In this work, VOF method, near-wall differential viscosity and local turbulence viscosity distribution coefficient function are introduced and offer an efficient tool to correct the interface turbulence viscosity. The results show that the new method can better predict the shear stresses, liquid holdup and pressure drop of stratified two-phase flow. The fitting relationship between interfacial and wall friction factor (fi /fW ) is in good agreement with the experimental data. It is found that fi /fW is predicted with a relative error of 12.62% by the new method, which is much less than that by any other method when the gas and liquid superficial Reynolds numbers are 8000≤ReSG ≤90000 and 5000≤ReSL ≤170000. It provides a reliable method for achieving the closure of stratified flow to predict the shear stresses. © 2019 Published by Elsevier Masson SAS.
1. Introduction Stratified gas–liquid flow in horizontal pipes has gained considerable importance because of its widespread applications in petroleum, chemical and nuclear engineering. The pressure drop, liquid holdup and velocity are the important factors of stratified gas–liquid flow, which are required to select the right pipe size and design appropriate facilities in chemical industry, steam generation, refrigeration, and nuclear reactor [1]. It should be stressed that stratified flow phenomena is strongly dominated not only by gas-wall and the liquid-wall shear stresses, but also by interfacial shear stress [2]. No unified conclusion regarding the closure relations of shear stresses was made due to complex momentum and energy transmission phenomena at gas–liquid interface. Directly numerical simulations for gas-wall, liquid-wall and interfacial shear stresses, which are in accord with the experimental observations cannot be obtained accurately in literature. There are two kinds of theoretical models for predicting stratified gas–liquid flow: global model and complementary local model. The objective of the global models is to calculate holdup and pressure drop in the entire stratified flow pattern [3]. However, the predictions of global models were unsatisfactory in some particular regions of flow regime [4]. Taitel and Dukler’s [5] one-dimensional two-fluid model, which was one of the earliest global models of stratified flow, gave good prediction in laminar ∗ Corresponding author. E-mail address: [email protected] (P. Zheng). https://doi.org/10.1016/j.euromechflu.2019.09.012 0997-7546/© 2019 Published by Elsevier Masson SAS.
and low turbulent flows, but relatively low accuracy was attained in larger gas–liquid turbulent flow due to its basic assumptions: the friction factor of gas–liquid interface was equal to that of the gas-wall, and liquid-wall friction factor was a single-phase-based result. Some investigators considered liquid-wall friction factor was failed to be determined by single-phase flow [6,7]. In this case, local two-fluid models accounting for the distribution of two phases and the shape of the gas–liquid interface were necessary. Tzotzi and Andritsos [7] developed a new local model based on Taitel and Dukler’s [5] model for wavy flow, in which stratified flow was divided into two-dimensional (2-D) wave and Kelvin– Helmholtz (K-H) wave sub-regimes. The friction factor equations for the two kinds of waves were presented. But there still existed great errors in high Reynolds number turbulent flow. Currently, with the rapid advance of computer technology, Computational Fluid Dynamics (CFD) approaches have become more and more affordable. Scholars have conducted much research on predicting the gas–liquid interfacial shear stress by CFD calculations. Meknassi et al. [8] expressed the gas–liquid interaction in terms of interfacial equivalent roughness and well predicted the experiments of Strand [9] and Lopez [10]. Newton and Behnia [11] assumed gas–liquid interface was smooth, and allowed the predictions of pressure gradient and liquid holdup by empirical wall damping functions. In 2001, they proposed an equation of empirical friction factor to predict the pressure drop and liquid holdup in wavy stratified flow [1]. However, their studies were mainly focused on gas–liquid interfacial shear stress, few works have considered the closure relations of stratified flow
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
to simulate all the different shear stresses including walls and gas–liquid interface. The surface tension calculation in the Volume of Fluid (VOF) method is improved from the standard continuum-surface-force (CSF) model in three different ways: the balanced force and pressure implementation [12–15], the improved interface representation [12,16–18] and height function modeling [19,20]. Renardy et al. [12] developed PROST algorithm which represented accurately the body force due to surface tension. The simulation showed convergence because the algorithm could reduce spurious currents. Francois et al. [13] proposed a VOF-based balanced-force algorithm to model interfacial flow with surface tension. They found the algorithm could achieve an exact balance of surface tension and pressure gradient forces. An equilibrium (static) drop in two and three dimensions was used to validate the method. But better methods for the curvature estimation were required before surface tension force was calculated. Popinet [14] presented a new combination of a VOF and adaptive quad/octree spatial discretization to accurately predict the surface-tensiondriven flows. A new balanced-force CSF surface tension scheme and height-function curvature estimation were used to recover an exact equilibrium solution between surface tension and pressure gradient by providing the shape of interface. Abu-Al-Saud et al. [15] presented a well-balanced surface tension model to predict interfacial flow. In the model, surface tension was considered and fluid momentum was conserved. The results showed that accurate solutions can be obtained for thermocapillary flows. Rider et al. [16] described the piecewise linear interface reconstruction and presented a new second-order-accurate algorithm for the volume tracking of interfaces. Pilliod et al. [17] presented a comprehensive framework for the design and implementation of VOF interface tracking algorithms. A second-order, unsplit, VOF advection algorithm was introduced and compared with the standard operator-split advection algorithm. The results were shown that the unsplit algorithm exhibited the better solution in some special regions. Liovic et al. [18] presented interface tracking-based solution algorithms to simulate surface tensiondominated interfacial flow. VOF interface tracking were used to develop curvature estimation and stability of interfacial flow simulation. Gerrits [19] firstly proposed the height function model to change the curvature calculation. This model could also reduce the spurious currents and the smoothed phase fraction, therefore improve the accuracy of surface tension calculation. Baltussen et al. [20] compared three different surface tension models for VOF method. The results showed that the height function method was the best for small and oscillating bubbles, the tensile force method was the best for large and stationary bubbles. And the height function and tensile force methods are both superior to the continuum-surface-force (CSF) method. Only a very limited number of works have dealt so far with stratified two-phase flow in pipelines by using VOF method. As we know, VOF method is volume-conserved and computes the volume fraction of a single phase instead of the interface. The main drawback of this method is spatial derivatives because the VOF function is discontinuous across the interface [21]. Some scholars improved the momentum source term and turbulence model to enhance the calculation accuracy of the stratified two-phase flow. Ghorai and Nigam [22] obtained more satisfactory solutions by employing interfacial equivalent roughness to revise the two-phase interface with VOF method, and a correlation which could predict the interfacial friction factor was presented. A three-dimensional gas-phase moving-wall model was presented with CFD software by Sidi-Ali and Gatignol [23], in which the gas–liquid interface was defined as non-slip moving wall with a steady velocity. The moving velocity, the interfacial friction factor and shear stress were simulated by using roughness velocity. However, the simulation error of gas-wall shear
203
stress was greater than 90% because it is assumed that the velocity of moving wall was equal to the average liquid velocity, which was not in accordance with the fact. In addition, no further information about gas-phase closure relation was given. In order to remedy the problem of conventional VOF method, Algebraic Interfacial Area Density (AIAD) model was adopted to revise SSW k − ω model in the simulation of free surface flow by Höhne and Mehlhoop [24], but the calculation accuracy was not of concern in the studies. Ningegowda et al. [25] concluded that the interfacial surface tension force calculation based on the LS function generates lesser parasitic currents as compared to the VOF function approach. The goal of this paper is to calculate the shear stresses, liquid holdup and pressure drop of stratified gas–liquid flow in pipelines, where the gas and liquid superficial Reynolds numbers are 8000 ≤ ReSG ≤ 90000 and 5000 ≤ ReSL ≤ 170000. RNG k − ε model and SST k − ω turbulent models are used during the simulation. Near-wall differential viscosity is taken into account. What is the most important is that a local turbulence viscosity, which involved a new distribution coefficient function is introduced to revise interfacial calculation. By VOF method, the volume fraction of a particular phase is calculated and the gas–liquid interface is tracked. All the predictions are conducted by ANSYS FLUENT 16.0 software. The gas velocity, wall and interfacial shear stresses and friction factors in stratified gas–liquid flow can be predicted accordingly. The numerical simulations are then compared with available experimental data and different methods. 2. Mathematical formulation 2.1. Stratified two-phase flow model 2.1.1. Momentum balance equations The stratified gas–liquid flow region in a near horizontal pipe is divided into two separate ‘conduits’ (see Fig. 1). The fully developed flow is expressed by two momentum balance equations for each phase [5]:
( −AG
dP
)
dz
( −AL
dP
)
dz
− τWG SG − τi Si + g ρG AG sin δ = 0
(1)
− τWL SL + τi Si + g ρL AL sin δ = 0
(2)
where δ is the angle of inclination of the pipe (equal to zero for horizontal pipe), g is the gravity acceleration, and ρG,L are the densities of gas (G) and liquid (L). If liquid height hL or liquid holdup α = AL /(AL + AG ) is decided, the cross-sectional areas of gas and liquid AG,L and the wetted perimeter of gas, liquid and interface SG,L,i can be obtained:
( ) 2hL θ = cos−1 1 −
(3)
D
AG =
1 4
π D2 − AL =
1 4
D2 (π − θ + cos θ sin θ)
(4.1)
SG = π D − SL = D (π − θ)
(4.2)
Si = D sin θ
(4.3)
where π is the circumference ratio, and D is the diameter of pipe. Taitel & Dukler transformed Eqs. (1) and (2) into a dimensionless form using eleven dimensionless parameters.
(dP /dz )L χ 2 == = (dP /dz )G
4 C D WL 4 C D WG
ρL USL D µL
)−mWL
ρG USG D µG
)−mWG
(
(
·
ρL |USL |USL
·
ρG |USG |USG
2
2
(5)
204
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
Fig. 1. Stratified gas–liquid flow in a near horizontal pipe.
where (dP/dz)L , G are pressure drop per unit length of single phase liquid flow and single phase gas flow, respectively. Lockhart– Martinelli parameter χ can be calculated unambiguously with the knowledge of flow rates, fluid properties, and pipe diameter [26]. Another parameter ˜ hL is the function of dimensionless liquid height, which can be written as follows:
and all cells completely filled with the gas phase are represented by F = 0. When the cell contains both gas and liquid phased, there will be a value between zero and one (0 < F <1). The advection equation for indicator function F is given as:
∇ (F U ) = F ∇ (U )
(15)
hL ˜ hL =
ρ = ρL F + ρG (1 − F )
(16)
µ = µL F + µG (1 − F )
(17)
(6)
D
Once gas-wall (WG), liquid-wall (WL) and interfacial shear stresses τWG,WL,i are determined, liquid holdup α and pressure drop −dP/dz can be calculated. Blasius-type equations based on single-phase flow to determine shear stresses were used:
τWG = fWG τWL = fWL τi = fi
ρG |UG | UG 2
ρL |UL | UL 2
,
U · ∇F = −
,
−mWG
fWG = CWG ReG
−mWL
fWL = CWL ReL
ρG |UG − UL | (UG − UL ) 2
,
(7) (8)
−mi
fi = Ci ReG
(9)
where UG,L are the average velocities, and fWG , WL , i are the friction factors. The Reynolds numbers of gas and liquid are, respectively, given as ReG = ρG DIG UG /µG and ReL = ρL DIL UL /µL , and hydraulic diameters of gas and liquid are, respectively, given as DIG = 4AG /(SG + Si ) and DIL = 4AL /SL . 2.1.2. VOF method For steady flow, the continuity and momentum equations of VOF method are given by Eqs. (10) and (11), respectively [21]:
∇ · (ρ U ) = 0
(10)
{ [ ]} ∇ · (ρ UU ) = −∇ P + ρ g + ∇ · µ ∇ U + (∇ U )T + Fσ
(11)
where U is the mean velocity vector, ρ is the density, µ is dynamic viscosity, P is the pressure, g is the gravity vector, the superscript T is the transpose operation of diffusion terms. Fσ is the volumetric surface tension force vector which can be determined through the continuum-surface-force (CSF) model [27]. Fσ = σ κδS n
κ=
1
|n|
[(
n
|n|
n = ∇ F /|∇ F |
)
As for the interface of two phases,
ρL
, U · ∇ (1 − F ) = −
S aG
ρG
(18)
where S aL , S aG are the external mass source into liquid and gas phase, which can be set as a constant or user-defined mass source term. In this work, the default value is zero [21]. 2.2. Turbulence modeling 2.2.1. RNG k − ε model RNG k − ε turbulence model is derived from the instantaneous Navier–Stokes equations, using a mathematical technique called renormalization group (RNG) methods [28]. RNG k − ε model is well suited for flows with large strain rates, large streamline curvature, low-Reynolds-number and near-wall flows [29], which has a similar form as the standard k − ε model. The transport equations for RNG k − ε model [30] are shown as:
∂ ∂ (ρ kUi ) = ∂ xj ∂ xj
[
∂ ∂ (ρε Ui ) = ∂ xj ∂ xj
[
] µeff ∂ k + Gk − ρε σk ∂ xj
(19)
] µeff ∂ε ε ε2 + C1ε Gk − C2ε ρ σε ∂ xj k k
(20)
where k is turbulence kinetic energy, ε is dissipation rate, Gk is calculated through: Gk = µt
(
∂ Uj ∂ Ui + ∂ xi ∂ xj
)
∂ Ui ∂ xj
(21)
(12)
the turbulent viscosity (eddy viscosity) µt is computed by combining k and ε as follows:
(13)
µt = ρ Cµ
]
∇ |n| − ∇ n
S aL
(14)
where σ is the surface tension, κ the interfacial curvature, δS is Dirac surface delta function given by |∇ F | and n is the normal vector. An indicator function F is defined for the volume fraction. All cells containing only the liquid phase will have a value of F = 1,
k2
(22)
ε
the effective viscosity µeff is:
( µeff = µ 1 +
√
Cµ k
µ
√
ε
)2 (23)
For the model constants the following values are chosen [30]: C1ε = 1.42, C2ε = 1.68, Cµ = 0.0845, σk = 1.0, σε = 1.3.
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
205
2.2.2. Shear-stress transport (SST) k − ω model SST k − ω model has been shown to be more accurate than standard k − ω model for a wider class of flows, including those with adverse pressure gradients [31]. The transport equations for k and ω for SST k − ω model [32] are shown as:
∂ ∂ (ρ kUi ) = ∂ xi ∂ xj
[(
µt µ+ σk
)
] ∂k + Gk − Yk ∂ xj
(24)
∂ ∂ (ρωUi ) = ∂ xj ∂ xj
[( ) ] µt ∂ω µ+ + Gω − Yω σω ∂ xj ρ ∂ k ∂ω + 2(1 − F1 ) ωσω,2 ∂ xj ∂ xj
(25)
where σk and σω are the turbulent Prandtl numbers for k and ω, which are shown as:
σk = σω =
1
(26)
F1 /σk,1 + (1 − F1 ) /σk,2 1
(27)
F1 /σω,1 + (1 − F1 ) /σω,2
F1 = tanh
⎧[ ⎨
( √
(
min max
(
;
500µ
0.09ωd ρ d2 ω
⎩
D+ ω = max
k
2ρ
∂ k ∂ω ; 10−10 σω,2 ω ∂ xj ∂ xj
) ;
4ρ k 2 σω,2 D+ ωd
)]4 ⎫ ⎬ ⎭
(
1
α∗
, aSFω2
F2 = tanh
⎩ α = ∗
∗ α∞
(
ρk Ret = µω
( max 2
500µ
)2 ⎫ ⎬
0.09ωd ρ d2 ω
⎭
k
α0∗ + Ret /Rk 1 + Ret /Rk
;
0 < γL < 1 γL = 1, γL = 0
(34)
FD = 1 − 2 2.0 + c1 + ec0 (|γL −0.5|−γ0 )
1
√
µt ,
(29)
where S is the strain rate magnitude, F2 is defined as:
⎧ ⎨
FD µt ,
where γL = 1 represents the fluid are all liquid phase, γL = 0 represents the fluid are all gas phase, and 0 < γL < 1 represents that the fluid are gas and liquid phase; FD is distribution coefficient function as follows:
(30)
)
{
(28)
)
1
k
ω max
no-slip wall. That is to say, the turbulence viscosity of moving interface tended to zero. Similar to the above methods, the interface turbulence viscosity was corrected by defining the smaller turbulence viscosity µ∗t in flow region to decrease the interface eddy diffusion and modify interface turbulence calculation. The smaller turbulence viscosity µ∗t is defined as:
µ∗t =
where d is the distance from the next surface. The turbulent viscosity (eddy viscosity) µt is defined as:
µt = ρ
Fig. 2. Distribution coefficient function FD vs. liquid volume fraction γL .
(31)
[
(33)
For the model constants the following values are chosen [32]: ∗ a1 = 0.31, Rk = 6, α0∗ = 0.024, α∞ = 1, σk,1 = 1.176, σk,2 = 1.0, σω,1 = 2.0, σω,2 = 1.168.
+ c2
(35)
where γ0 is the range of the minimum values, here γ0 = 0.4, the changing speed of FD depends on c0 , the smaller c0 is, the more intense FD changes, here c0 = 70; c1 and c2 are revised to achieve FD = 1 when γL = 0 or γL = 1; when c1 = 1.65×10−2 and c2 = 1.82×10−3 , µt * reaches the minimum value of 0.01, FD is defined with the help of UDF (user defined function) for turbulence viscosity and shown in Fig. 2. 2.3.2. Near-wall differential viscosity modification Given near-wall flow is low-Reynolds-number flow, the improved differential viscosity in near wall region was used as follows:
) (32)
]−1
ρ2k d √ εµ
(
1.72
) = √(
µeff µ
(
)3
µeff µ
)
− 1 + Cν
( d
µeff µ
) (36)
where Cν ≈ 100. The modification of turbulence viscosity in Eqs. (34) ∼(36) can enhance the VOF calculation accuracy of mixing flows in near wall region.
2.3. Correction equations 2.3.1. Local turbulence viscosity correction Not considering the momentum transfer at the gas–liquid interface in the VOF method, only a single-phase momentum equation is calculated. As a result, the interface turbulence viscosity is too large to predict the interfacial shear stress in continuousphase simulation. Issa [33] increased dissipation rate ε by assuming the gradient of turbulence kinetic energy k was equal to zero, and thereby the turbulence viscosity at the interface was extremely small. Accordingly, his simulations are in good agreement with the experimental results. As pointed out by SidiAli and Gatignol [23], the gas–liquid interface was regarded as
2.3.3. Turbulence damping correction Since SST k − ω model has little effect on the simulation of the near-interface region, which has the similar results as k −ε model, appreciable errors may occur. A new turbulence damping Sq is introduced to enhance the interface calculation accuracy [24], which is added to the right-hand side of Eq. (30), reads
Sq = Aq ∆nβρq
6Bµq
(
βρq ∆n2
Aq = 2.0γq ⏐∇γq ⏐
⏐
⏐
)2 (37)
(38)
206
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
where Aq is the area density of phase q, ∆n is the cell height normal to interface, β = 0.075, B is damping coefficient, whose default value is equal to 10, γq is volume fraction of phase q.
Table 1 Meshes for 2D model.
∆z/m
∆y/m
∆z/∆y
Iteration step number
2.4. Boundary conditions
1 2 3
80 104 40 083 30 070
0.012 0.024 0.032
0.0045 0.0045 0.0045
2.67 5.33 7.11
5000 3000 2400
At the inlet, the gas superficial velocity and the liquid superficial velocity are given by:
4 5 6
105 944 54 306 35 648
0.012 0.024 0.032
0.0030 0.0030 0.0030
4.00 8.00 10.67
5600 4200 2800
USG = USG0 , USL = USL0
(40)
3. Numerical method In the present work, the commercial software ANSYS Fluent 16.0 is used as the CFD solver. Single-precision, pressurebased, steady and explicit solver are used for the calculation due to steady stratified interface. To resolve the velocity–pressure coupling, PISO segregated algorithm is employed. The secondorder discretization scheme is employed for momentum, because it can be used to obtain better results for tetrahedral and quad/hex meshes, especially for complex flows [21]. The Modified HRIC discretization is used for calculation of volume fraction. The under-relaxation factors of momentum, density and body force are lowered to 0.8 for better convergence [21]. Residuals of continuity, momentum, energy of the solution, velocity and turbulent kinetic energy and turbulence dissipation equations are performed within 1.0 × 10−5 relative error as the criterion for convergence of the solution. 3.1. Two-dimensional pipe In our paper, we firstly establish a two-dimensional pipe model according to Sidi-Ali and Gatignol’s [23] experimental data. The pipe diameter is 0.1m and the length is 31 m. The gas superficial velocity USG0 is 2.18 m/s, the liquid superficial velocity USL0 is 0.16 m/s. The surface tension coefficient of gas and liquid is 0.0736 N m−1 . Meshes and different methods were calculated and compared in two-dimensional pipe, which is the evidence of choosing the appropriate three-dimensional pipe model. 3.1.1. Meshes generation Different meshes of two-dimensional pipe are shown in Table 1, where ∆z is length of mesh and ∆y is width of mesh. The ratio of mesh length and width ∆z /∆y has a significant effect on the simulation. If ∆y is set to a constant value, when the velocity is low, ∆z /∆y should be decreased; when the velocity is high, ∆z /∆y should be increased. For Meshes 1 ∼ 3, ∆y maintains 0.0045; for Meshes 4 ∼ 6, ∆y maintains 0.0030, and ∆z is changing from 0.012, 0.024 to 0.036. Sidi-Ali and Gatignol’s [23] model well predicted the interfacial shear stress, but Taitel and Dukler’s model [5] has higher accuracy in predicting gas-wall and liquid-wall shear stress. For Sidi-Ali and Gatignol’s [23] model, the correlations for interfacial shear stress and friction factor can be written as: 0.522 τi = 0.5fi ρG (UG − UL )2 , fi = 2.82Re− G
(41)
For Taitel and Dukler’s model [5], in Eqs. (7) and (8), CWG,WL = 0.046, mWG,WL = −0.2. Thus, the correlations for wall shear stress and friction factor can be written as:
τq = 0.5fq ρq Uq 2 , fq = 0.046Req −0.2
Nodes number
(39)
An atmospheric pressure is the imposed pressure at the exit boundary, and the diffusion flux for all variables in exit direction is equal to zero:
∂ (U , k, ε, ω) = 0 ∂n
Meshes
(42)
where i represents interface of gas and liquid, q represents gas and liquid phase. By Eqs. (1) ∼ (6) and Eqs. (41) and (42), the simulation value of liquid holdup α is 0.353, and the pressure drop dP/dz is 1.98 Pa/m. Relative errors of liquid holdup α and pressure drop -(dP/dz) calculated by Meshes 1, 2 and 4 are shown in Fig. 3(a). It is found that the liquid holdup α calculated by Meshes 3 is 0.283, by Meshes 5 is 0.241 and by Meshes 6 is 0.215, they are rather low and dissatisfactory, therefore they are all neglected in Fig. 3(a). Relative errors of gas-wall and liquid-wall shear stress obtained by different meshes are shown in Fig. 3(b). It can be seen that extremely low ∆z /∆y has a great effect on the prediction of wall shear stress and pressure drop. For Meshes 1, the relative error of the pressure drop -(dP/dz) is extremely high, which is not reasonable. The best predictions of pressure drop -(dP/dz) occur in Meshes 2 and Meshes 4. For Meshes 2 and Meshes 4, the liquid holdup α calculated by Meshes 2 is closed to that of Meshes 4; however, the wall shear stress calculated by Meshes 2 is much reasonable. As a result, Meshes 2 is the most suitable meshes for calculation. 3.1.2. Confirmation of the best simulation method Seven different simulation methods of two-dimensional pipe are shown in Table 2. (1) Verification test Assuming the z axis is the direction of length of the pipe, the liquid volume fraction of the cross section at z = 26 m and outlet cross section at z = 31 m changing with iteration steps are observed and showed as Fig. 4. At first, liquid volume fraction γL of the cross section increases while step number increases by M 2 ∼ M 4 methods. In Fig. 4(a), as steps number arrives at about 3000, liquid volume fraction γL of the cross section at z = 26 m and outlet cross section at z = 31 m by M 2 ∼ M 3 methods are stable in a horizontal line. In Fig. 4(b), as steps number arrives at about 15 000, liquid volume fraction γL of the cross section at z = 26 m and outlet cross section at z = 31 m by M 4a ∼ M 4c methods are also stable. The stratified flow can be considered steady and the convergence is achieved. (2) Comparisons of different simulation methods Fig. 5(a) shows relative errors of liquid holdup α and pressure drop -(dP/dz) calculated by different methods. Fig. 5(b) shows relative errors of gas-wall and liquid-wall shear stresses in different methods. The relative errors of pressure drop and gas-wall shear stress predicted by M 4a ∼ M 4c are larger than that predicted using M 1 ∼ M 3 method. That means SST k − ω model cannot accurately meet Taitel and Dukler’s [5] experimental results. From Fig. 5(a), it can be seen that M 2 has the lowest relative errors of liquid holdup and pressure drop in M 1 ∼ M 4 methods. By means of RNG k − ε model containing the local turbulence viscosity and near-wall differential viscosity, the relative errors of liquid holdup and pressure drop are −1.41% and −4.46%, the relative errors of gas-wall and liquid-wall shear stress are 22.7% and −14.1%. As a result, RNG k − ε model is more conducive to discussing closure relationship of stratified flow than SST k − ω model. The simulations predicted by M 2 method provide the closest results to Taitel and Dukler’s [5] model.
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
207
Table 2 Cases of different turbulence methods. Methods M M M M M M
1 2 3 4a 4b 4c
RNG k-ε
SST k-ω
√ √
Local turbulence viscosity and near-wall differential viscosity
Turbulence damping
Iteration step number
√ (B = 0.6) √ (B = 0.75) √ (B = 1.0)
2 800 3 000 3 000 15 000 17 600 19 000
√ √ √ √ √
Fig. 3. Relative errors calculated by different meshes.
Fig. 4. Liquid volume fraction of the cross section at z = 26 m and z = 31 m changing with iteration step number.
As we know, there occurs a problem for calculation of interface in stratified flow using VOF method. If the surface tension force plays an important role, an artificial velocity field at the interface called ‘‘parasitic currents’’ may arise. Parasitic currents can destabilize the interface significantly. In this way can we conclude that the same method as M 2 method is used for the simulation of three-dimensional pipe.
was done with the center line liquid volume fraction γL−C vs. the y axis. It can be observed that the difference is very small indicating the mesh independence of the solution. In view of computational accuracy and cost, the 1.3 million mesh (refinement) for Run 1 and the 1.55 million mesh (refinement) for Run 2 are adopted. The meshes for Run 1 is shown as Fig. 7. 4. Results and discussion
3.2. Three-dimensional pipe 4.1. Gas velocity Based on the previous conclusions of two-dimensional model and the existing experimental results [9], three-dimensional pipe model is established. The gas superficial velocity are 1.7 m/s (ReSG = 11,800) and 5.5 m/s (ReSG = 38,200), the liquid superficial velocity is 0.1 m/s (ReSL = 13,000). The observation region is limited by the length of pipe z ∈ [25,26]. The two flows are defined as Run 1 and Run 2. The structured meshes are used and local mesh refinement near the gas–liquid wall is taken into account. Four types of mesh cells number are selected for Run 1 and Run 2, and the mesh independence study, as shown in Fig. 6,
Fig. 8 represents the stream wise velocity Uz , which is calculated by M 1 method, M 2 method and Strand’s [9] experimental results. We can see that the maximal gas velocity calculated by M 1 method is obviously greater than the experimental results and that calculated with M 2 method. M 2 method improves the prediction of gas velocity gradient at the interface which is in good agreement with Strand’s [9] experimental results. Ayati et al. [34] measured the mean axial velocity composed by the gas velocity and liquid profile by advanced PIV method,
208
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
Fig. 5. Relative errors calculated by different methods.
Fig. 6. Mesh independence study in 3D models.
Fig. 9. Mean velocity profiles measured by Ayati et al. [34].
as shown in Fig. 9. We can see the results show that the mean axial velocity has the same profile as the stream wise velocity Uz in Fig. 8 and agrees better with our simulations. Fig. 10 represents Uy , the center line gas velocity in the direction of the y axis (z = 26 m), which is calculated by M 1 method, M 2 method and Meknassi et al.’s [8] simulation results. It is assumed that Y is non-dimensional height as follows [24]: Y = Fig. 7. 3D horizontal pipe model and its meshes.
y − hL D 2
− hL
(43)
where D is diameter of pipes, hL is liquid height. It can be observed that the velocity of secondary flow is rather low with M 1 method, which is one order of magnitude smaller than that calculated by Meknassi et al. [8]. This is due to the fact that M 2 method takes account of the secondary flow and the gas velocity simulated by M 2 method agrees well with Meknassi et al.’s [8] simulation. 4.2. Wall shear stress As we know, the pressure drop of stratified flow is directly affected by the wall friction. The interface shear stress hinders the gas flow, whereas it improves the liquid flow [5,34]. Accordingly, the key of predicting pressure drop is to enhance the prediction accuracy of liquid holdup and wall shear stress. The average value of total shear stress τW can be given by: Fig. 8. Comparison of the stream wise velocity Uz as a function of the vertical coordinate y.
τW =
SG τG + SL τL SG + SL
(44)
where SG,L are the gas and liquid wetted perimeters, τG and τL can be obtained after interfacial height is calculated. The equivalent interfacial shear stress τi∗ is determined by the following
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
209
Fig. 10. The center line component velocity in the direction of the y axis of gas flow.
relations: 2τi∗ Si = (AG − AL )
( −
dP
)
dz
− τG SG + τL SL
(45)
where Si is the length of cross section interface, AG,L are gas and liquid cross section area. Table 3 shows the simulations of liquid holdup, pressure drop, wall shear stress and relative errors calculated by M 1 method, M 2 method and Strand’s [9] experimental results. We can see that the liquid holdup calculated by M 1 method is lower than Strand’s [9] experimental results, and the pressure drop errors are rather high. For Run1 and Run2, the relative errors of shear stress, liquid holdup and pressure drop calculated by M 2 method are much lower than those using M 1 method. M 2 method can successfully predict the shear stress, liquid holdup and pressure drop of three-dimensional pipes.
Fig. 11. Comparisons of fi /fW with M 2 method, Ghorai and Nigam’s model and Strand’s experiments.
4.3. Prediction of friction factor relations Assuming the superficial gas velocity USG is from 1.2 m/s to 8.6 m/s, the superficial liquid velocity USL is 0.1 m/s, the relation between wall friction factor fW and equivalent interfacial friction factor fi can be derived by Ghorai and Nigam’s [22] model as follows: fi fW
= 1.204
(
USG
)0.245
USL
(46)
fi /fW can be derived by M 2 method as follows: fi fW
= 0.669USG 1.094 USL 0.169
(47)
Fig. 11 shows fi /fW calculated by M 2 method, Ghorai and Nigam’s [22] model and Strand’s [9] experiments. The superficial gas velocity USG is from 4.19 m/s to 12.90 m/s, the superficial liquid velocity USL is from 0.05 m/s to 0.13 m/s. The Reynoldsnumber range of Ghorai and Nigam’s [28] model is 8000 ≤ ReSG ≤ 90,000 and 5000 ≤ ReSL ≤ 17,000. The average root error of fi /fW calculated by Ghorai and Nigam’s [22] model is 34.52%, while the average root error of fi /fW calculated by M 2 method is 12.62%. M 2 method predicts the better results of friction factor than those with Ghorai and Nigam’s [22] model. Fig. 12 shows fi /fW vs. the superficial liquid velocity USL under the same or similar superficial gas velocity USG condition, which is in good agreement with the results of Meknassi et al. [8]. As USL was increased, fi /fW is found to increase when USG > 5.0 m/s; as USL was decreased, fi /fW is found to decrease when USG < 5.0 m/s. When USG = 5.0 m/s, fi /fW has no obvious trend. Since 83% of Strand’s experiments were performed when USG > 5.0 m/s, the positive exponent of USL in Eq. (47) presents more reasonable results when USG > 5.0 m/s.
Fig. 12. Comparison of fi /fW values with different liquid superficial velocities.
5. Conclusions Numerical simulation of RNG k − ε turbulence model and VOF method containing the local turbulence viscosity and near-wall differential viscosity was used for the analysis of two-dimensional and three-dimensional stratified gas–liquid flow in horizontal pipes. Mesh independence study and validation of numerical simulation were both carried out. Calculated results of gas velocity, wall shear stress were compared with Strand’s experimental values and with Meknassi et al.’s [8] theoretic prediction. The prediction of friction factor relations are also obtained by this new method. The calculated results showed that RNG k − ε model can be more beneficial to achieve the closure of stratified flow. Not depending on experimental results, a new distribution
210
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211 Table 3 Results of M 1 method, M 2 method and experimental data. Run
1
2
Method
Result
τG (Pa)
τL (Pa)
τi or τi∗ (Pa)
α (–)
Strand
Value
0.035
0.166
0.083
0.386
3.00
M 1
Value Error/%
0.039 +11.4
0.181 +9.04
0.143 +72.3
0.324 −16.1
3.95 +31.7
M 2
Value Error/%
0.037
0.153
+5.71
−8.00
0.086 +3.61
0.371 −3.88
2.82 −6.00
Strand
Value
0.189
0.582
0.470
0.252
13.8
M 1
Value Error/%
0.112 −40.7
0.408 −29.9
0.391 −16.8
0.218 −13.6
10.8 −21.7
M 2
Value Error/%
0.179
0.575
−5.29
−1.13
0.437 −6.91
0.264 +4.53
13.0 −5.80
coefficient function FD was defined. Calculated velocity of stratified flow is in good agreement with experimental values. The new method successfully predicted the shear stress (including gas-wall, liquid-wall, equivalent interfacial shear stress), liquid holdup and pressure drop, the simulation error of fi /fW is 12.62% with 8000 ≤ ReSG ≤ 90,000 and 5000 ≤ ReSL ≤ 17,000. Nomenclature Aq AG AL B C D d -dP/dz F FD Fσ
[kg/m3 ] [m2 ] [m2 ] [–] [–] [m] [m] [Pa/m] [–] [–] [N/m3 ]
fi fW g G Gb
[–] [–] [m/s2 ] [–] [–]
Gk Gω h hL ˜ hL
[–] [–] [m] [m] [–]
k L n
[m2 /s2 ] [–] [–]
P q Re ReSG ReSL S SG SL S aG
[Pa] [–] [–] [–] [–] [–] [m] [m] [–]
S aL
[–]
Sq Si
[–] [m]
The area density of phase q Cross-sectional area of gas Cross-sectional area of liquid Damping coefficient Coefficient Diameter of pipes The distance from the nearest wall Pressure drop per unit length Indicator function Distribution coefficient function The volumetric surface tension force vector Interfacial friction factor Wall facial friction factor Gravity vector Gas phase The generation of k due to the mean velocity The generation of k due to buoyancy The generation of ω Grid spacing Liquid height Function of dimensionless liquid height Turbulence kinetic energy Liquid phase The divergence of the interface unit normal Pressure Phase (gas or liquid) Reynolds number Gas superficial Reynolds number Liquid superficial Reynolds number Strain rate magnitude Wetted perimeter of gas Wetted perimeter of liquid The external mass source into gas phase The external mass source into liquid phase Turbulence damping The length of cross section interface
U USG USL Uy
[m/s] [m/s] [m/s] [m/s]
Uz x, y, z Y YM
[m/s] [m] [–] [–]
Yk [–] Yω [–] Greek symbols α [–] αk [–]
αε
[–]
β δ δS γL−C ε ρ ρq ρG ρL µ µG µL µt µt ∗ µeff τG τL τi ∗ τij τW ω ∆z ∆y χ
[–] [rad] [–] [–] [m2 /s3 ] [kg/m3 ] [kg/m3 ] [kg/m3 ] [kg/m3 ] [Pa s] [Pa s] [Pa s] [Pa s] [Pa s] [Pa s] [Pa] [Pa] [Pa] [Pa] [Pa] [s−1 ] [m] [m] [–]
−(dP/dz) (Pa m−1 )
Velocity gradient vector Gas superficial velocity Liquid superficial velocity Center line gas velocity in the direction of y Stream wise velocity Cartesian coordinate Non-dimensional height The contribution of the fluctuating dilatation to the dissipation rate The dissipation of k The dissipation of ω Liquid holdup The inverse effective Prandtl numbers for k The inverse effective Prandtl numbers for ε Constant Angle of inclination of the pipe The delta function Center line of liquid volume fraction Turbulent dissipation rate Density Density of phase q Density of gas Density of liquid Dynamic viscosity Dynamic viscosity of gas Dynamic viscosity of liquid Turbulent (or eddy) viscosity Smaller turbulence viscosity Effective viscosity Shear stress of gas Shear stress of liquid Equivalent interfacial shear stress Turbulent stress tensor Average value of total shear stress Specific dissipation rate Length of mesh Width of mesh Lockhart–Martinell parameter
Acknowledgments The authors are profoundly grateful to the financial supports of the Scientific Technology Research Project of Education Department of Liaoning Province, China (Grant No. L2016017, No. L2017LZD004) and Natural Science Foundation of Liaoning Province, China (Grant No. 20170540591)
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
References [1] C.H. Newton, M. Behnia, A numerical model of stratified-wavy gas-liquid pipe flow, Chem. Eng. Sci. 56 (2001) 6851–6861. [2] A. Ullmann, N. Brauner, Closure relations for two-fluid models for twophase stratified smooth and stratified wavy flows, Int. J. Multiphas. Flow 32 (2006) 82–105. [3] A. Liné, D. Lopez, Two-fluid model of wavy separated two-phase flow, Int. J. Multiphas. Flow 23 (1997) 1131–1146. [4] P.L. Spedding, E. Benard, G.F. Donnelly, Prediction of pressure drop in Multiphas. horizontal pipe flow, Int. Commun. Heat Mass 33 (2006) 1053–1062. [5] Y. Taitel, A.E. Dukler, A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow, AIChE J. 22 (1976) 47–55. [6] J.E. Kowalski, Wall and interfacial shear stress in stratified flow in a horizontal pipe, AIChE J. 32 (1987) 274–281. [7] C. Tzotzi, N. Andritsos, Interfacial shear stress in wavy stratified gas-liquid flow in horizontal Pipes, Int. J. Multiphas. Flow 54 (2013) 43–54. [8] F. Meknassi, R. Benkirane, A. Liné, L. Masbernat, Numerical modelling of wavy stratified two-phase flow in pipes, Chem. Eng. Sci. 55 (2000) 4681–4697. [9] O. Strand, An Experimental Investigation of Stratified Two-Phase Flow in Horizontal Pipes (Ph.D. Thesis), University of Oslo, Norway, 1993. [10] D. Lopez, Ecoulements diphasiques a phases separees a faible contenu de liquide, I.N.P. Toulouse, France, 1994. [11] C.H. Newton, M. Behnia, Numerical calculation of turbulent stratified-gasliquid pipe flows, Int. J. Multiphas. Flow 26 (2000) 327–337. [12] Y. Renardy, M. Renardy, Prost: a parabolic reconstruction of surface tension for the volume-of-fluid method, J. Comput. Phys. 183 (2002) 400–421. [13] M.M. Francois, S.J. Cummins, E.D. Dendy, D.B. Kothe, J.M. Sicilian, M.W. Williams, A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comput. Phys. 213 (2006) 141–173. [14] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Phys. 228 (2009) 5838–5866. [15] M.O. Abu-Al-Saud, S. Popinet, H.A. Tchelepi, A conservative and well-balanced surface tension model, J. Comput. Phys. 371 (2018) 896–913. [16] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (1998) 112–152. [17] J.E. Pilliod, E.G. Puckett, Second-order accurate volume-of-fluid algorithms for tracking material interfaces, J. Comput. Phys. 199 (2004) 465–502.
211
[18] P. Liovic, M. Francois, M. Rudman, R. Manasseh, Efficient simulation of surface tension-dominated flows through enhanced interface geometry interrogation, J. Comput. Phys. 229 (2010) 7520–7544. [19] J. Gerrits, Dynamics of Liquid-Filled Spacecraft: Numerical Simulation of Coupled Solid–Liquid Dynamics, University of Groningen, 2001. [20] M.W. Baltussen, J.A.M. Kuipers, N.G. Deen, A critical comparison of surface tension models for the volume of fluid method, Chem. Eng. Sci. 109 (2014) 65–74. [21] ANSYS FLUENT, ANSYS FLUENT Theory Guide Release 16.0. ANSYS Inc., 2014. [22] S. Ghorai, K.D.P. Nigam, CFD modelling of flow profiles and interfacial phenomena in two-phase flow in pipes, Chem. Eng. Process. 45 (2006) 55–65. [23] K. Sidi-Ali, R. Gatignol, Interfacial friction factor determination using CFD simulations in a horizontal stratified two-phase flow, Chem. Eng. Sci. 65 (2010) 5160–5169. [24] T. Höhne, J.-P. Mehlhoop, Validation of closure models for interfacial drag and turbulence in numerical simulations of horizontal stratified gas-liquid flows, Int. J. Multiphas. Flow 62 (2014) 1–16. [25] B.M. Ningegowda, B. Premachandran, A Coupled Level Set and Volume of Fluid method with multi-directional advection algorithms for two-phase flows with and without phase change, Int. J. Heat Mass Transfer 79 (2014) 532–550. [26] R.W. Lockhart, R.C. Martinelli, Proposed correlation of data for isothermal two-phase, two-component flow in pipes, Chem. Eng. Prog. 45 (1949) 39–48. [27] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [28] S.B. Pope, Turbulent Flows, Cambridge University Press, 2000. [29] C.L. William, R.S. James, An analysis of practical RANS simulations for spinnaker aerodynamics, J. Wind. Eng. Ind. Aerodyn. 96 (2008) 143–165. [30] Q.F. Hou, Z.S. Zou, Comparison between standard and renormalization group k-ε models in numerical simulation of swirling flow tundish, ISIJ Int. 45 (2005) 325–330. [31] C. Laurent, I. Mary, V. Gleize, A. Lerat, D. Arnal, DNS database of a transitional separation bubble on a flat plate and application to RANS modeling validation, Comput. Fluids 61 (2012) 21–30. [32] D.C. Wilcox, Turbulence Modeling for CFD, DCW Industries, Inc., La Canada, California, 2006. [33] R.I. Issa, Prediction of turbulent, stratified, two-phase flow in inclined pipes and channel, Int. J. Multiphas. Flow 14 (1988) 141–154. [34] A.A. Ayati, J. Kolaas, A. Jensen, G.W. Johnson, A PIV investigation of stratified gas-liquid flow in a horizontal pipe, Int. J. Multiphas. Flow 61 (2014) 129–143.
Contents lists available at ScienceDirect
European Journal of Mechanics / B Fluids journal homepage: www.elsevier.com/locate/ejmflu
Correcting interface turbulence viscosity using CFD modeling for predicting stratified gas–liquid flow shear stress in horizontal pipes ∗
Ping Zheng a , , Liang Zhao b a b
College of Petroleum Engineering, Liaoning Shihua University, Fushun, Liaoning 113001, China College of Civil Aviation Safety Engineering, Civil Aviation Flight University of China, Guanghan, Sichuan 618307, China
article
info
Article history: Received 30 June 2017 Received in revised form 6 September 2019 Accepted 17 September 2019 Available online 19 September 2019 Keywords: Stratified gas–liquid flow Turbulence viscosity Liquid holdup Pressure drop Shear stress CFD modeling
a b s t r a c t The prediction of different shear stresses is one of the great challenges of turbulent–turbulent stratified two-phase flow in horizontal pipes. In this work, VOF method, near-wall differential viscosity and local turbulence viscosity distribution coefficient function are introduced and offer an efficient tool to correct the interface turbulence viscosity. The results show that the new method can better predict the shear stresses, liquid holdup and pressure drop of stratified two-phase flow. The fitting relationship between interfacial and wall friction factor (fi /fW ) is in good agreement with the experimental data. It is found that fi /fW is predicted with a relative error of 12.62% by the new method, which is much less than that by any other method when the gas and liquid superficial Reynolds numbers are 8000≤ReSG ≤90000 and 5000≤ReSL ≤170000. It provides a reliable method for achieving the closure of stratified flow to predict the shear stresses. © 2019 Published by Elsevier Masson SAS.
1. Introduction Stratified gas–liquid flow in horizontal pipes has gained considerable importance because of its widespread applications in petroleum, chemical and nuclear engineering. The pressure drop, liquid holdup and velocity are the important factors of stratified gas–liquid flow, which are required to select the right pipe size and design appropriate facilities in chemical industry, steam generation, refrigeration, and nuclear reactor [1]. It should be stressed that stratified flow phenomena is strongly dominated not only by gas-wall and the liquid-wall shear stresses, but also by interfacial shear stress [2]. No unified conclusion regarding the closure relations of shear stresses was made due to complex momentum and energy transmission phenomena at gas–liquid interface. Directly numerical simulations for gas-wall, liquid-wall and interfacial shear stresses, which are in accord with the experimental observations cannot be obtained accurately in literature. There are two kinds of theoretical models for predicting stratified gas–liquid flow: global model and complementary local model. The objective of the global models is to calculate holdup and pressure drop in the entire stratified flow pattern [3]. However, the predictions of global models were unsatisfactory in some particular regions of flow regime [4]. Taitel and Dukler’s [5] one-dimensional two-fluid model, which was one of the earliest global models of stratified flow, gave good prediction in laminar ∗ Corresponding author. E-mail address: [email protected] (P. Zheng). https://doi.org/10.1016/j.euromechflu.2019.09.012 0997-7546/© 2019 Published by Elsevier Masson SAS.
and low turbulent flows, but relatively low accuracy was attained in larger gas–liquid turbulent flow due to its basic assumptions: the friction factor of gas–liquid interface was equal to that of the gas-wall, and liquid-wall friction factor was a single-phase-based result. Some investigators considered liquid-wall friction factor was failed to be determined by single-phase flow [6,7]. In this case, local two-fluid models accounting for the distribution of two phases and the shape of the gas–liquid interface were necessary. Tzotzi and Andritsos [7] developed a new local model based on Taitel and Dukler’s [5] model for wavy flow, in which stratified flow was divided into two-dimensional (2-D) wave and Kelvin– Helmholtz (K-H) wave sub-regimes. The friction factor equations for the two kinds of waves were presented. But there still existed great errors in high Reynolds number turbulent flow. Currently, with the rapid advance of computer technology, Computational Fluid Dynamics (CFD) approaches have become more and more affordable. Scholars have conducted much research on predicting the gas–liquid interfacial shear stress by CFD calculations. Meknassi et al. [8] expressed the gas–liquid interaction in terms of interfacial equivalent roughness and well predicted the experiments of Strand [9] and Lopez [10]. Newton and Behnia [11] assumed gas–liquid interface was smooth, and allowed the predictions of pressure gradient and liquid holdup by empirical wall damping functions. In 2001, they proposed an equation of empirical friction factor to predict the pressure drop and liquid holdup in wavy stratified flow [1]. However, their studies were mainly focused on gas–liquid interfacial shear stress, few works have considered the closure relations of stratified flow
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
to simulate all the different shear stresses including walls and gas–liquid interface. The surface tension calculation in the Volume of Fluid (VOF) method is improved from the standard continuum-surface-force (CSF) model in three different ways: the balanced force and pressure implementation [12–15], the improved interface representation [12,16–18] and height function modeling [19,20]. Renardy et al. [12] developed PROST algorithm which represented accurately the body force due to surface tension. The simulation showed convergence because the algorithm could reduce spurious currents. Francois et al. [13] proposed a VOF-based balanced-force algorithm to model interfacial flow with surface tension. They found the algorithm could achieve an exact balance of surface tension and pressure gradient forces. An equilibrium (static) drop in two and three dimensions was used to validate the method. But better methods for the curvature estimation were required before surface tension force was calculated. Popinet [14] presented a new combination of a VOF and adaptive quad/octree spatial discretization to accurately predict the surface-tensiondriven flows. A new balanced-force CSF surface tension scheme and height-function curvature estimation were used to recover an exact equilibrium solution between surface tension and pressure gradient by providing the shape of interface. Abu-Al-Saud et al. [15] presented a well-balanced surface tension model to predict interfacial flow. In the model, surface tension was considered and fluid momentum was conserved. The results showed that accurate solutions can be obtained for thermocapillary flows. Rider et al. [16] described the piecewise linear interface reconstruction and presented a new second-order-accurate algorithm for the volume tracking of interfaces. Pilliod et al. [17] presented a comprehensive framework for the design and implementation of VOF interface tracking algorithms. A second-order, unsplit, VOF advection algorithm was introduced and compared with the standard operator-split advection algorithm. The results were shown that the unsplit algorithm exhibited the better solution in some special regions. Liovic et al. [18] presented interface tracking-based solution algorithms to simulate surface tensiondominated interfacial flow. VOF interface tracking were used to develop curvature estimation and stability of interfacial flow simulation. Gerrits [19] firstly proposed the height function model to change the curvature calculation. This model could also reduce the spurious currents and the smoothed phase fraction, therefore improve the accuracy of surface tension calculation. Baltussen et al. [20] compared three different surface tension models for VOF method. The results showed that the height function method was the best for small and oscillating bubbles, the tensile force method was the best for large and stationary bubbles. And the height function and tensile force methods are both superior to the continuum-surface-force (CSF) method. Only a very limited number of works have dealt so far with stratified two-phase flow in pipelines by using VOF method. As we know, VOF method is volume-conserved and computes the volume fraction of a single phase instead of the interface. The main drawback of this method is spatial derivatives because the VOF function is discontinuous across the interface [21]. Some scholars improved the momentum source term and turbulence model to enhance the calculation accuracy of the stratified two-phase flow. Ghorai and Nigam [22] obtained more satisfactory solutions by employing interfacial equivalent roughness to revise the two-phase interface with VOF method, and a correlation which could predict the interfacial friction factor was presented. A three-dimensional gas-phase moving-wall model was presented with CFD software by Sidi-Ali and Gatignol [23], in which the gas–liquid interface was defined as non-slip moving wall with a steady velocity. The moving velocity, the interfacial friction factor and shear stress were simulated by using roughness velocity. However, the simulation error of gas-wall shear
203
stress was greater than 90% because it is assumed that the velocity of moving wall was equal to the average liquid velocity, which was not in accordance with the fact. In addition, no further information about gas-phase closure relation was given. In order to remedy the problem of conventional VOF method, Algebraic Interfacial Area Density (AIAD) model was adopted to revise SSW k − ω model in the simulation of free surface flow by Höhne and Mehlhoop [24], but the calculation accuracy was not of concern in the studies. Ningegowda et al. [25] concluded that the interfacial surface tension force calculation based on the LS function generates lesser parasitic currents as compared to the VOF function approach. The goal of this paper is to calculate the shear stresses, liquid holdup and pressure drop of stratified gas–liquid flow in pipelines, where the gas and liquid superficial Reynolds numbers are 8000 ≤ ReSG ≤ 90000 and 5000 ≤ ReSL ≤ 170000. RNG k − ε model and SST k − ω turbulent models are used during the simulation. Near-wall differential viscosity is taken into account. What is the most important is that a local turbulence viscosity, which involved a new distribution coefficient function is introduced to revise interfacial calculation. By VOF method, the volume fraction of a particular phase is calculated and the gas–liquid interface is tracked. All the predictions are conducted by ANSYS FLUENT 16.0 software. The gas velocity, wall and interfacial shear stresses and friction factors in stratified gas–liquid flow can be predicted accordingly. The numerical simulations are then compared with available experimental data and different methods. 2. Mathematical formulation 2.1. Stratified two-phase flow model 2.1.1. Momentum balance equations The stratified gas–liquid flow region in a near horizontal pipe is divided into two separate ‘conduits’ (see Fig. 1). The fully developed flow is expressed by two momentum balance equations for each phase [5]:
( −AG
dP
)
dz
( −AL
dP
)
dz
− τWG SG − τi Si + g ρG AG sin δ = 0
(1)
− τWL SL + τi Si + g ρL AL sin δ = 0
(2)
where δ is the angle of inclination of the pipe (equal to zero for horizontal pipe), g is the gravity acceleration, and ρG,L are the densities of gas (G) and liquid (L). If liquid height hL or liquid holdup α = AL /(AL + AG ) is decided, the cross-sectional areas of gas and liquid AG,L and the wetted perimeter of gas, liquid and interface SG,L,i can be obtained:
( ) 2hL θ = cos−1 1 −
(3)
D
AG =
1 4
π D2 − AL =
1 4
D2 (π − θ + cos θ sin θ)
(4.1)
SG = π D − SL = D (π − θ)
(4.2)
Si = D sin θ
(4.3)
where π is the circumference ratio, and D is the diameter of pipe. Taitel & Dukler transformed Eqs. (1) and (2) into a dimensionless form using eleven dimensionless parameters.
(dP /dz )L χ 2 == = (dP /dz )G
4 C D WL 4 C D WG
ρL USL D µL
)−mWL
ρG USG D µG
)−mWG
(
(
·
ρL |USL |USL
·
ρG |USG |USG
2
2
(5)
204
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
Fig. 1. Stratified gas–liquid flow in a near horizontal pipe.
where (dP/dz)L , G are pressure drop per unit length of single phase liquid flow and single phase gas flow, respectively. Lockhart– Martinelli parameter χ can be calculated unambiguously with the knowledge of flow rates, fluid properties, and pipe diameter [26]. Another parameter ˜ hL is the function of dimensionless liquid height, which can be written as follows:
and all cells completely filled with the gas phase are represented by F = 0. When the cell contains both gas and liquid phased, there will be a value between zero and one (0 < F <1). The advection equation for indicator function F is given as:
∇ (F U ) = F ∇ (U )
(15)
hL ˜ hL =
ρ = ρL F + ρG (1 − F )
(16)
µ = µL F + µG (1 − F )
(17)
(6)
D
Once gas-wall (WG), liquid-wall (WL) and interfacial shear stresses τWG,WL,i are determined, liquid holdup α and pressure drop −dP/dz can be calculated. Blasius-type equations based on single-phase flow to determine shear stresses were used:
τWG = fWG τWL = fWL τi = fi
ρG |UG | UG 2
ρL |UL | UL 2
,
U · ∇F = −
,
−mWG
fWG = CWG ReG
−mWL
fWL = CWL ReL
ρG |UG − UL | (UG − UL ) 2
,
(7) (8)
−mi
fi = Ci ReG
(9)
where UG,L are the average velocities, and fWG , WL , i are the friction factors. The Reynolds numbers of gas and liquid are, respectively, given as ReG = ρG DIG UG /µG and ReL = ρL DIL UL /µL , and hydraulic diameters of gas and liquid are, respectively, given as DIG = 4AG /(SG + Si ) and DIL = 4AL /SL . 2.1.2. VOF method For steady flow, the continuity and momentum equations of VOF method are given by Eqs. (10) and (11), respectively [21]:
∇ · (ρ U ) = 0
(10)
{ [ ]} ∇ · (ρ UU ) = −∇ P + ρ g + ∇ · µ ∇ U + (∇ U )T + Fσ
(11)
where U is the mean velocity vector, ρ is the density, µ is dynamic viscosity, P is the pressure, g is the gravity vector, the superscript T is the transpose operation of diffusion terms. Fσ is the volumetric surface tension force vector which can be determined through the continuum-surface-force (CSF) model [27]. Fσ = σ κδS n
κ=
1
|n|
[(
n
|n|
n = ∇ F /|∇ F |
)
As for the interface of two phases,
ρL
, U · ∇ (1 − F ) = −
S aG
ρG
(18)
where S aL , S aG are the external mass source into liquid and gas phase, which can be set as a constant or user-defined mass source term. In this work, the default value is zero [21]. 2.2. Turbulence modeling 2.2.1. RNG k − ε model RNG k − ε turbulence model is derived from the instantaneous Navier–Stokes equations, using a mathematical technique called renormalization group (RNG) methods [28]. RNG k − ε model is well suited for flows with large strain rates, large streamline curvature, low-Reynolds-number and near-wall flows [29], which has a similar form as the standard k − ε model. The transport equations for RNG k − ε model [30] are shown as:
∂ ∂ (ρ kUi ) = ∂ xj ∂ xj
[
∂ ∂ (ρε Ui ) = ∂ xj ∂ xj
[
] µeff ∂ k + Gk − ρε σk ∂ xj
(19)
] µeff ∂ε ε ε2 + C1ε Gk − C2ε ρ σε ∂ xj k k
(20)
where k is turbulence kinetic energy, ε is dissipation rate, Gk is calculated through: Gk = µt
(
∂ Uj ∂ Ui + ∂ xi ∂ xj
)
∂ Ui ∂ xj
(21)
(12)
the turbulent viscosity (eddy viscosity) µt is computed by combining k and ε as follows:
(13)
µt = ρ Cµ
]
∇ |n| − ∇ n
S aL
(14)
where σ is the surface tension, κ the interfacial curvature, δS is Dirac surface delta function given by |∇ F | and n is the normal vector. An indicator function F is defined for the volume fraction. All cells containing only the liquid phase will have a value of F = 1,
k2
(22)
ε
the effective viscosity µeff is:
( µeff = µ 1 +
√
Cµ k
µ
√
ε
)2 (23)
For the model constants the following values are chosen [30]: C1ε = 1.42, C2ε = 1.68, Cµ = 0.0845, σk = 1.0, σε = 1.3.
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
205
2.2.2. Shear-stress transport (SST) k − ω model SST k − ω model has been shown to be more accurate than standard k − ω model for a wider class of flows, including those with adverse pressure gradients [31]. The transport equations for k and ω for SST k − ω model [32] are shown as:
∂ ∂ (ρ kUi ) = ∂ xi ∂ xj
[(
µt µ+ σk
)
] ∂k + Gk − Yk ∂ xj
(24)
∂ ∂ (ρωUi ) = ∂ xj ∂ xj
[( ) ] µt ∂ω µ+ + Gω − Yω σω ∂ xj ρ ∂ k ∂ω + 2(1 − F1 ) ωσω,2 ∂ xj ∂ xj
(25)
where σk and σω are the turbulent Prandtl numbers for k and ω, which are shown as:
σk = σω =
1
(26)
F1 /σk,1 + (1 − F1 ) /σk,2 1
(27)
F1 /σω,1 + (1 − F1 ) /σω,2
F1 = tanh
⎧[ ⎨
( √
(
min max
(
;
500µ
0.09ωd ρ d2 ω
⎩
D+ ω = max
k
2ρ
∂ k ∂ω ; 10−10 σω,2 ω ∂ xj ∂ xj
) ;
4ρ k 2 σω,2 D+ ωd
)]4 ⎫ ⎬ ⎭
(
1
α∗
, aSFω2
F2 = tanh
⎩ α = ∗
∗ α∞
(
ρk Ret = µω
( max 2
500µ
)2 ⎫ ⎬
0.09ωd ρ d2 ω
⎭
k
α0∗ + Ret /Rk 1 + Ret /Rk
;
0 < γL < 1 γL = 1, γL = 0
(34)
FD = 1 − 2 2.0 + c1 + ec0 (|γL −0.5|−γ0 )
1
√
µt ,
(29)
where S is the strain rate magnitude, F2 is defined as:
⎧ ⎨
FD µt ,
where γL = 1 represents the fluid are all liquid phase, γL = 0 represents the fluid are all gas phase, and 0 < γL < 1 represents that the fluid are gas and liquid phase; FD is distribution coefficient function as follows:
(30)
)
{
(28)
)
1
k
ω max
no-slip wall. That is to say, the turbulence viscosity of moving interface tended to zero. Similar to the above methods, the interface turbulence viscosity was corrected by defining the smaller turbulence viscosity µ∗t in flow region to decrease the interface eddy diffusion and modify interface turbulence calculation. The smaller turbulence viscosity µ∗t is defined as:
µ∗t =
where d is the distance from the next surface. The turbulent viscosity (eddy viscosity) µt is defined as:
µt = ρ
Fig. 2. Distribution coefficient function FD vs. liquid volume fraction γL .
(31)
[
(33)
For the model constants the following values are chosen [32]: ∗ a1 = 0.31, Rk = 6, α0∗ = 0.024, α∞ = 1, σk,1 = 1.176, σk,2 = 1.0, σω,1 = 2.0, σω,2 = 1.168.
+ c2
(35)
where γ0 is the range of the minimum values, here γ0 = 0.4, the changing speed of FD depends on c0 , the smaller c0 is, the more intense FD changes, here c0 = 70; c1 and c2 are revised to achieve FD = 1 when γL = 0 or γL = 1; when c1 = 1.65×10−2 and c2 = 1.82×10−3 , µt * reaches the minimum value of 0.01, FD is defined with the help of UDF (user defined function) for turbulence viscosity and shown in Fig. 2. 2.3.2. Near-wall differential viscosity modification Given near-wall flow is low-Reynolds-number flow, the improved differential viscosity in near wall region was used as follows:
) (32)
]−1
ρ2k d √ εµ
(
1.72
) = √(
µeff µ
(
)3
µeff µ
)
− 1 + Cν
( d
µeff µ
) (36)
where Cν ≈ 100. The modification of turbulence viscosity in Eqs. (34) ∼(36) can enhance the VOF calculation accuracy of mixing flows in near wall region.
2.3. Correction equations 2.3.1. Local turbulence viscosity correction Not considering the momentum transfer at the gas–liquid interface in the VOF method, only a single-phase momentum equation is calculated. As a result, the interface turbulence viscosity is too large to predict the interfacial shear stress in continuousphase simulation. Issa [33] increased dissipation rate ε by assuming the gradient of turbulence kinetic energy k was equal to zero, and thereby the turbulence viscosity at the interface was extremely small. Accordingly, his simulations are in good agreement with the experimental results. As pointed out by SidiAli and Gatignol [23], the gas–liquid interface was regarded as
2.3.3. Turbulence damping correction Since SST k − ω model has little effect on the simulation of the near-interface region, which has the similar results as k −ε model, appreciable errors may occur. A new turbulence damping Sq is introduced to enhance the interface calculation accuracy [24], which is added to the right-hand side of Eq. (30), reads
Sq = Aq ∆nβρq
6Bµq
(
βρq ∆n2
Aq = 2.0γq ⏐∇γq ⏐
⏐
⏐
)2 (37)
(38)
206
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
where Aq is the area density of phase q, ∆n is the cell height normal to interface, β = 0.075, B is damping coefficient, whose default value is equal to 10, γq is volume fraction of phase q.
Table 1 Meshes for 2D model.
∆z/m
∆y/m
∆z/∆y
Iteration step number
2.4. Boundary conditions
1 2 3
80 104 40 083 30 070
0.012 0.024 0.032
0.0045 0.0045 0.0045
2.67 5.33 7.11
5000 3000 2400
At the inlet, the gas superficial velocity and the liquid superficial velocity are given by:
4 5 6
105 944 54 306 35 648
0.012 0.024 0.032
0.0030 0.0030 0.0030
4.00 8.00 10.67
5600 4200 2800
USG = USG0 , USL = USL0
(40)
3. Numerical method In the present work, the commercial software ANSYS Fluent 16.0 is used as the CFD solver. Single-precision, pressurebased, steady and explicit solver are used for the calculation due to steady stratified interface. To resolve the velocity–pressure coupling, PISO segregated algorithm is employed. The secondorder discretization scheme is employed for momentum, because it can be used to obtain better results for tetrahedral and quad/hex meshes, especially for complex flows [21]. The Modified HRIC discretization is used for calculation of volume fraction. The under-relaxation factors of momentum, density and body force are lowered to 0.8 for better convergence [21]. Residuals of continuity, momentum, energy of the solution, velocity and turbulent kinetic energy and turbulence dissipation equations are performed within 1.0 × 10−5 relative error as the criterion for convergence of the solution. 3.1. Two-dimensional pipe In our paper, we firstly establish a two-dimensional pipe model according to Sidi-Ali and Gatignol’s [23] experimental data. The pipe diameter is 0.1m and the length is 31 m. The gas superficial velocity USG0 is 2.18 m/s, the liquid superficial velocity USL0 is 0.16 m/s. The surface tension coefficient of gas and liquid is 0.0736 N m−1 . Meshes and different methods were calculated and compared in two-dimensional pipe, which is the evidence of choosing the appropriate three-dimensional pipe model. 3.1.1. Meshes generation Different meshes of two-dimensional pipe are shown in Table 1, where ∆z is length of mesh and ∆y is width of mesh. The ratio of mesh length and width ∆z /∆y has a significant effect on the simulation. If ∆y is set to a constant value, when the velocity is low, ∆z /∆y should be decreased; when the velocity is high, ∆z /∆y should be increased. For Meshes 1 ∼ 3, ∆y maintains 0.0045; for Meshes 4 ∼ 6, ∆y maintains 0.0030, and ∆z is changing from 0.012, 0.024 to 0.036. Sidi-Ali and Gatignol’s [23] model well predicted the interfacial shear stress, but Taitel and Dukler’s model [5] has higher accuracy in predicting gas-wall and liquid-wall shear stress. For Sidi-Ali and Gatignol’s [23] model, the correlations for interfacial shear stress and friction factor can be written as: 0.522 τi = 0.5fi ρG (UG − UL )2 , fi = 2.82Re− G
(41)
For Taitel and Dukler’s model [5], in Eqs. (7) and (8), CWG,WL = 0.046, mWG,WL = −0.2. Thus, the correlations for wall shear stress and friction factor can be written as:
τq = 0.5fq ρq Uq 2 , fq = 0.046Req −0.2
Nodes number
(39)
An atmospheric pressure is the imposed pressure at the exit boundary, and the diffusion flux for all variables in exit direction is equal to zero:
∂ (U , k, ε, ω) = 0 ∂n
Meshes
(42)
where i represents interface of gas and liquid, q represents gas and liquid phase. By Eqs. (1) ∼ (6) and Eqs. (41) and (42), the simulation value of liquid holdup α is 0.353, and the pressure drop dP/dz is 1.98 Pa/m. Relative errors of liquid holdup α and pressure drop -(dP/dz) calculated by Meshes 1, 2 and 4 are shown in Fig. 3(a). It is found that the liquid holdup α calculated by Meshes 3 is 0.283, by Meshes 5 is 0.241 and by Meshes 6 is 0.215, they are rather low and dissatisfactory, therefore they are all neglected in Fig. 3(a). Relative errors of gas-wall and liquid-wall shear stress obtained by different meshes are shown in Fig. 3(b). It can be seen that extremely low ∆z /∆y has a great effect on the prediction of wall shear stress and pressure drop. For Meshes 1, the relative error of the pressure drop -(dP/dz) is extremely high, which is not reasonable. The best predictions of pressure drop -(dP/dz) occur in Meshes 2 and Meshes 4. For Meshes 2 and Meshes 4, the liquid holdup α calculated by Meshes 2 is closed to that of Meshes 4; however, the wall shear stress calculated by Meshes 2 is much reasonable. As a result, Meshes 2 is the most suitable meshes for calculation. 3.1.2. Confirmation of the best simulation method Seven different simulation methods of two-dimensional pipe are shown in Table 2. (1) Verification test Assuming the z axis is the direction of length of the pipe, the liquid volume fraction of the cross section at z = 26 m and outlet cross section at z = 31 m changing with iteration steps are observed and showed as Fig. 4. At first, liquid volume fraction γL of the cross section increases while step number increases by M 2 ∼ M 4 methods. In Fig. 4(a), as steps number arrives at about 3000, liquid volume fraction γL of the cross section at z = 26 m and outlet cross section at z = 31 m by M 2 ∼ M 3 methods are stable in a horizontal line. In Fig. 4(b), as steps number arrives at about 15 000, liquid volume fraction γL of the cross section at z = 26 m and outlet cross section at z = 31 m by M 4a ∼ M 4c methods are also stable. The stratified flow can be considered steady and the convergence is achieved. (2) Comparisons of different simulation methods Fig. 5(a) shows relative errors of liquid holdup α and pressure drop -(dP/dz) calculated by different methods. Fig. 5(b) shows relative errors of gas-wall and liquid-wall shear stresses in different methods. The relative errors of pressure drop and gas-wall shear stress predicted by M 4a ∼ M 4c are larger than that predicted using M 1 ∼ M 3 method. That means SST k − ω model cannot accurately meet Taitel and Dukler’s [5] experimental results. From Fig. 5(a), it can be seen that M 2 has the lowest relative errors of liquid holdup and pressure drop in M 1 ∼ M 4 methods. By means of RNG k − ε model containing the local turbulence viscosity and near-wall differential viscosity, the relative errors of liquid holdup and pressure drop are −1.41% and −4.46%, the relative errors of gas-wall and liquid-wall shear stress are 22.7% and −14.1%. As a result, RNG k − ε model is more conducive to discussing closure relationship of stratified flow than SST k − ω model. The simulations predicted by M 2 method provide the closest results to Taitel and Dukler’s [5] model.
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
207
Table 2 Cases of different turbulence methods. Methods M M M M M M
1 2 3 4a 4b 4c
RNG k-ε
SST k-ω
√ √
Local turbulence viscosity and near-wall differential viscosity
Turbulence damping
Iteration step number
√ (B = 0.6) √ (B = 0.75) √ (B = 1.0)
2 800 3 000 3 000 15 000 17 600 19 000
√ √ √ √ √
Fig. 3. Relative errors calculated by different meshes.
Fig. 4. Liquid volume fraction of the cross section at z = 26 m and z = 31 m changing with iteration step number.
As we know, there occurs a problem for calculation of interface in stratified flow using VOF method. If the surface tension force plays an important role, an artificial velocity field at the interface called ‘‘parasitic currents’’ may arise. Parasitic currents can destabilize the interface significantly. In this way can we conclude that the same method as M 2 method is used for the simulation of three-dimensional pipe.
was done with the center line liquid volume fraction γL−C vs. the y axis. It can be observed that the difference is very small indicating the mesh independence of the solution. In view of computational accuracy and cost, the 1.3 million mesh (refinement) for Run 1 and the 1.55 million mesh (refinement) for Run 2 are adopted. The meshes for Run 1 is shown as Fig. 7. 4. Results and discussion
3.2. Three-dimensional pipe 4.1. Gas velocity Based on the previous conclusions of two-dimensional model and the existing experimental results [9], three-dimensional pipe model is established. The gas superficial velocity are 1.7 m/s (ReSG = 11,800) and 5.5 m/s (ReSG = 38,200), the liquid superficial velocity is 0.1 m/s (ReSL = 13,000). The observation region is limited by the length of pipe z ∈ [25,26]. The two flows are defined as Run 1 and Run 2. The structured meshes are used and local mesh refinement near the gas–liquid wall is taken into account. Four types of mesh cells number are selected for Run 1 and Run 2, and the mesh independence study, as shown in Fig. 6,
Fig. 8 represents the stream wise velocity Uz , which is calculated by M 1 method, M 2 method and Strand’s [9] experimental results. We can see that the maximal gas velocity calculated by M 1 method is obviously greater than the experimental results and that calculated with M 2 method. M 2 method improves the prediction of gas velocity gradient at the interface which is in good agreement with Strand’s [9] experimental results. Ayati et al. [34] measured the mean axial velocity composed by the gas velocity and liquid profile by advanced PIV method,
208
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
Fig. 5. Relative errors calculated by different methods.
Fig. 6. Mesh independence study in 3D models.
Fig. 9. Mean velocity profiles measured by Ayati et al. [34].
as shown in Fig. 9. We can see the results show that the mean axial velocity has the same profile as the stream wise velocity Uz in Fig. 8 and agrees better with our simulations. Fig. 10 represents Uy , the center line gas velocity in the direction of the y axis (z = 26 m), which is calculated by M 1 method, M 2 method and Meknassi et al.’s [8] simulation results. It is assumed that Y is non-dimensional height as follows [24]: Y = Fig. 7. 3D horizontal pipe model and its meshes.
y − hL D 2
− hL
(43)
where D is diameter of pipes, hL is liquid height. It can be observed that the velocity of secondary flow is rather low with M 1 method, which is one order of magnitude smaller than that calculated by Meknassi et al. [8]. This is due to the fact that M 2 method takes account of the secondary flow and the gas velocity simulated by M 2 method agrees well with Meknassi et al.’s [8] simulation. 4.2. Wall shear stress As we know, the pressure drop of stratified flow is directly affected by the wall friction. The interface shear stress hinders the gas flow, whereas it improves the liquid flow [5,34]. Accordingly, the key of predicting pressure drop is to enhance the prediction accuracy of liquid holdup and wall shear stress. The average value of total shear stress τW can be given by: Fig. 8. Comparison of the stream wise velocity Uz as a function of the vertical coordinate y.
τW =
SG τG + SL τL SG + SL
(44)
where SG,L are the gas and liquid wetted perimeters, τG and τL can be obtained after interfacial height is calculated. The equivalent interfacial shear stress τi∗ is determined by the following
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
209
Fig. 10. The center line component velocity in the direction of the y axis of gas flow.
relations: 2τi∗ Si = (AG − AL )
( −
dP
)
dz
− τG SG + τL SL
(45)
where Si is the length of cross section interface, AG,L are gas and liquid cross section area. Table 3 shows the simulations of liquid holdup, pressure drop, wall shear stress and relative errors calculated by M 1 method, M 2 method and Strand’s [9] experimental results. We can see that the liquid holdup calculated by M 1 method is lower than Strand’s [9] experimental results, and the pressure drop errors are rather high. For Run1 and Run2, the relative errors of shear stress, liquid holdup and pressure drop calculated by M 2 method are much lower than those using M 1 method. M 2 method can successfully predict the shear stress, liquid holdup and pressure drop of three-dimensional pipes.
Fig. 11. Comparisons of fi /fW with M 2 method, Ghorai and Nigam’s model and Strand’s experiments.
4.3. Prediction of friction factor relations Assuming the superficial gas velocity USG is from 1.2 m/s to 8.6 m/s, the superficial liquid velocity USL is 0.1 m/s, the relation between wall friction factor fW and equivalent interfacial friction factor fi can be derived by Ghorai and Nigam’s [22] model as follows: fi fW
= 1.204
(
USG
)0.245
USL
(46)
fi /fW can be derived by M 2 method as follows: fi fW
= 0.669USG 1.094 USL 0.169
(47)
Fig. 11 shows fi /fW calculated by M 2 method, Ghorai and Nigam’s [22] model and Strand’s [9] experiments. The superficial gas velocity USG is from 4.19 m/s to 12.90 m/s, the superficial liquid velocity USL is from 0.05 m/s to 0.13 m/s. The Reynoldsnumber range of Ghorai and Nigam’s [28] model is 8000 ≤ ReSG ≤ 90,000 and 5000 ≤ ReSL ≤ 17,000. The average root error of fi /fW calculated by Ghorai and Nigam’s [22] model is 34.52%, while the average root error of fi /fW calculated by M 2 method is 12.62%. M 2 method predicts the better results of friction factor than those with Ghorai and Nigam’s [22] model. Fig. 12 shows fi /fW vs. the superficial liquid velocity USL under the same or similar superficial gas velocity USG condition, which is in good agreement with the results of Meknassi et al. [8]. As USL was increased, fi /fW is found to increase when USG > 5.0 m/s; as USL was decreased, fi /fW is found to decrease when USG < 5.0 m/s. When USG = 5.0 m/s, fi /fW has no obvious trend. Since 83% of Strand’s experiments were performed when USG > 5.0 m/s, the positive exponent of USL in Eq. (47) presents more reasonable results when USG > 5.0 m/s.
Fig. 12. Comparison of fi /fW values with different liquid superficial velocities.
5. Conclusions Numerical simulation of RNG k − ε turbulence model and VOF method containing the local turbulence viscosity and near-wall differential viscosity was used for the analysis of two-dimensional and three-dimensional stratified gas–liquid flow in horizontal pipes. Mesh independence study and validation of numerical simulation were both carried out. Calculated results of gas velocity, wall shear stress were compared with Strand’s experimental values and with Meknassi et al.’s [8] theoretic prediction. The prediction of friction factor relations are also obtained by this new method. The calculated results showed that RNG k − ε model can be more beneficial to achieve the closure of stratified flow. Not depending on experimental results, a new distribution
210
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211 Table 3 Results of M 1 method, M 2 method and experimental data. Run
1
2
Method
Result
τG (Pa)
τL (Pa)
τi or τi∗ (Pa)
α (–)
Strand
Value
0.035
0.166
0.083
0.386
3.00
M 1
Value Error/%
0.039 +11.4
0.181 +9.04
0.143 +72.3
0.324 −16.1
3.95 +31.7
M 2
Value Error/%
0.037
0.153
+5.71
−8.00
0.086 +3.61
0.371 −3.88
2.82 −6.00
Strand
Value
0.189
0.582
0.470
0.252
13.8
M 1
Value Error/%
0.112 −40.7
0.408 −29.9
0.391 −16.8
0.218 −13.6
10.8 −21.7
M 2
Value Error/%
0.179
0.575
−5.29
−1.13
0.437 −6.91
0.264 +4.53
13.0 −5.80
coefficient function FD was defined. Calculated velocity of stratified flow is in good agreement with experimental values. The new method successfully predicted the shear stress (including gas-wall, liquid-wall, equivalent interfacial shear stress), liquid holdup and pressure drop, the simulation error of fi /fW is 12.62% with 8000 ≤ ReSG ≤ 90,000 and 5000 ≤ ReSL ≤ 17,000. Nomenclature Aq AG AL B C D d -dP/dz F FD Fσ
[kg/m3 ] [m2 ] [m2 ] [–] [–] [m] [m] [Pa/m] [–] [–] [N/m3 ]
fi fW g G Gb
[–] [–] [m/s2 ] [–] [–]
Gk Gω h hL ˜ hL
[–] [–] [m] [m] [–]
k L n
[m2 /s2 ] [–] [–]
P q Re ReSG ReSL S SG SL S aG
[Pa] [–] [–] [–] [–] [–] [m] [m] [–]
S aL
[–]
Sq Si
[–] [m]
The area density of phase q Cross-sectional area of gas Cross-sectional area of liquid Damping coefficient Coefficient Diameter of pipes The distance from the nearest wall Pressure drop per unit length Indicator function Distribution coefficient function The volumetric surface tension force vector Interfacial friction factor Wall facial friction factor Gravity vector Gas phase The generation of k due to the mean velocity The generation of k due to buoyancy The generation of ω Grid spacing Liquid height Function of dimensionless liquid height Turbulence kinetic energy Liquid phase The divergence of the interface unit normal Pressure Phase (gas or liquid) Reynolds number Gas superficial Reynolds number Liquid superficial Reynolds number Strain rate magnitude Wetted perimeter of gas Wetted perimeter of liquid The external mass source into gas phase The external mass source into liquid phase Turbulence damping The length of cross section interface
U USG USL Uy
[m/s] [m/s] [m/s] [m/s]
Uz x, y, z Y YM
[m/s] [m] [–] [–]
Yk [–] Yω [–] Greek symbols α [–] αk [–]
αε
[–]
β δ δS γL−C ε ρ ρq ρG ρL µ µG µL µt µt ∗ µeff τG τL τi ∗ τij τW ω ∆z ∆y χ
[–] [rad] [–] [–] [m2 /s3 ] [kg/m3 ] [kg/m3 ] [kg/m3 ] [kg/m3 ] [Pa s] [Pa s] [Pa s] [Pa s] [Pa s] [Pa s] [Pa] [Pa] [Pa] [Pa] [Pa] [s−1 ] [m] [m] [–]
−(dP/dz) (Pa m−1 )
Velocity gradient vector Gas superficial velocity Liquid superficial velocity Center line gas velocity in the direction of y Stream wise velocity Cartesian coordinate Non-dimensional height The contribution of the fluctuating dilatation to the dissipation rate The dissipation of k The dissipation of ω Liquid holdup The inverse effective Prandtl numbers for k The inverse effective Prandtl numbers for ε Constant Angle of inclination of the pipe The delta function Center line of liquid volume fraction Turbulent dissipation rate Density Density of phase q Density of gas Density of liquid Dynamic viscosity Dynamic viscosity of gas Dynamic viscosity of liquid Turbulent (or eddy) viscosity Smaller turbulence viscosity Effective viscosity Shear stress of gas Shear stress of liquid Equivalent interfacial shear stress Turbulent stress tensor Average value of total shear stress Specific dissipation rate Length of mesh Width of mesh Lockhart–Martinell parameter
Acknowledgments The authors are profoundly grateful to the financial supports of the Scientific Technology Research Project of Education Department of Liaoning Province, China (Grant No. L2016017, No. L2017LZD004) and Natural Science Foundation of Liaoning Province, China (Grant No. 20170540591)
P. Zheng and L. Zhao / European Journal of Mechanics / B Fluids 79 (2020) 202–211
References [1] C.H. Newton, M. Behnia, A numerical model of stratified-wavy gas-liquid pipe flow, Chem. Eng. Sci. 56 (2001) 6851–6861. [2] A. Ullmann, N. Brauner, Closure relations for two-fluid models for twophase stratified smooth and stratified wavy flows, Int. J. Multiphas. Flow 32 (2006) 82–105. [3] A. Liné, D. Lopez, Two-fluid model of wavy separated two-phase flow, Int. J. Multiphas. Flow 23 (1997) 1131–1146. [4] P.L. Spedding, E. Benard, G.F. Donnelly, Prediction of pressure drop in Multiphas. horizontal pipe flow, Int. Commun. Heat Mass 33 (2006) 1053–1062. [5] Y. Taitel, A.E. Dukler, A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow, AIChE J. 22 (1976) 47–55. [6] J.E. Kowalski, Wall and interfacial shear stress in stratified flow in a horizontal pipe, AIChE J. 32 (1987) 274–281. [7] C. Tzotzi, N. Andritsos, Interfacial shear stress in wavy stratified gas-liquid flow in horizontal Pipes, Int. J. Multiphas. Flow 54 (2013) 43–54. [8] F. Meknassi, R. Benkirane, A. Liné, L. Masbernat, Numerical modelling of wavy stratified two-phase flow in pipes, Chem. Eng. Sci. 55 (2000) 4681–4697. [9] O. Strand, An Experimental Investigation of Stratified Two-Phase Flow in Horizontal Pipes (Ph.D. Thesis), University of Oslo, Norway, 1993. [10] D. Lopez, Ecoulements diphasiques a phases separees a faible contenu de liquide, I.N.P. Toulouse, France, 1994. [11] C.H. Newton, M. Behnia, Numerical calculation of turbulent stratified-gasliquid pipe flows, Int. J. Multiphas. Flow 26 (2000) 327–337. [12] Y. Renardy, M. Renardy, Prost: a parabolic reconstruction of surface tension for the volume-of-fluid method, J. Comput. Phys. 183 (2002) 400–421. [13] M.M. Francois, S.J. Cummins, E.D. Dendy, D.B. Kothe, J.M. Sicilian, M.W. Williams, A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comput. Phys. 213 (2006) 141–173. [14] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Phys. 228 (2009) 5838–5866. [15] M.O. Abu-Al-Saud, S. Popinet, H.A. Tchelepi, A conservative and well-balanced surface tension model, J. Comput. Phys. 371 (2018) 896–913. [16] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (1998) 112–152. [17] J.E. Pilliod, E.G. Puckett, Second-order accurate volume-of-fluid algorithms for tracking material interfaces, J. Comput. Phys. 199 (2004) 465–502.
211
[18] P. Liovic, M. Francois, M. Rudman, R. Manasseh, Efficient simulation of surface tension-dominated flows through enhanced interface geometry interrogation, J. Comput. Phys. 229 (2010) 7520–7544. [19] J. Gerrits, Dynamics of Liquid-Filled Spacecraft: Numerical Simulation of Coupled Solid–Liquid Dynamics, University of Groningen, 2001. [20] M.W. Baltussen, J.A.M. Kuipers, N.G. Deen, A critical comparison of surface tension models for the volume of fluid method, Chem. Eng. Sci. 109 (2014) 65–74. [21] ANSYS FLUENT, ANSYS FLUENT Theory Guide Release 16.0. ANSYS Inc., 2014. [22] S. Ghorai, K.D.P. Nigam, CFD modelling of flow profiles and interfacial phenomena in two-phase flow in pipes, Chem. Eng. Process. 45 (2006) 55–65. [23] K. Sidi-Ali, R. Gatignol, Interfacial friction factor determination using CFD simulations in a horizontal stratified two-phase flow, Chem. Eng. Sci. 65 (2010) 5160–5169. [24] T. Höhne, J.-P. Mehlhoop, Validation of closure models for interfacial drag and turbulence in numerical simulations of horizontal stratified gas-liquid flows, Int. J. Multiphas. Flow 62 (2014) 1–16. [25] B.M. Ningegowda, B. Premachandran, A Coupled Level Set and Volume of Fluid method with multi-directional advection algorithms for two-phase flows with and without phase change, Int. J. Heat Mass Transfer 79 (2014) 532–550. [26] R.W. Lockhart, R.C. Martinelli, Proposed correlation of data for isothermal two-phase, two-component flow in pipes, Chem. Eng. Prog. 45 (1949) 39–48. [27] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [28] S.B. Pope, Turbulent Flows, Cambridge University Press, 2000. [29] C.L. William, R.S. James, An analysis of practical RANS simulations for spinnaker aerodynamics, J. Wind. Eng. Ind. Aerodyn. 96 (2008) 143–165. [30] Q.F. Hou, Z.S. Zou, Comparison between standard and renormalization group k-ε models in numerical simulation of swirling flow tundish, ISIJ Int. 45 (2005) 325–330. [31] C. Laurent, I. Mary, V. Gleize, A. Lerat, D. Arnal, DNS database of a transitional separation bubble on a flat plate and application to RANS modeling validation, Comput. Fluids 61 (2012) 21–30. [32] D.C. Wilcox, Turbulence Modeling for CFD, DCW Industries, Inc., La Canada, California, 2006. [33] R.I. Issa, Prediction of turbulent, stratified, two-phase flow in inclined pipes and channel, Int. J. Multiphas. Flow 14 (1988) 141–154. [34] A.A. Ayati, J. Kolaas, A. Jensen, G.W. Johnson, A PIV investigation of stratified gas-liquid flow in a horizontal pipe, Int. J. Multiphas. Flow 61 (2014) 129–143.