Article pubs.acs.org/IECR
Heat Transfer and Phase Change Characterization of Multicomponent Vapor/Condensation Stratified Flow in a Hydrocarbon Pipeline Guoxi He* National Engineering Laboratory for Pipeline Safety/MOE Key Laboratory of Petroleum Engineering/Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum-Beijing, Beijing 102249, People’s Republic of China ABSTRACT: Heat transfer together with phase change in pipelines frequently occurs in petroleum and chemical industries. Numerical methods have been extensively used in modeling two-phase flow and become more complicated when the phase change process is involved. One-dimensional two-fluid models are unreliable for modeling phase behavior because they do only yield uniform velocity/temperature distributions. Thus, mathematical models based on bipolar coordinate systems have been developed, where nonisothermal flows and heat/mass transfer have been implemented to predict the varying vapor−liquid interface and component mass fraction. The PR equation is chosen for phase equilibrium calculation and the P−T flash is used to evaluate the physical properties, gasification rate, and enthalpy departure. The liquid holdup, velocity/temperature profile, and phase fraction could be predicted along the pipeline, which are important for transportation optimization and flow assurance. The numerical predictions fit well with available experimental and simulated data and could present more detailed and accurate results. vapor−liquid two-phase flow.14,15 The fundamental engineering parameters are the pressure drop, liquid holdup, phase fraction, phase flow rate, temperature, thermo-physical properties of the fluids, and pipeline geometry.16 Despite several decades of research into vapor−liquid pipe flow phenomena, the numerical calculation of these parameters in stratified fluctuating vapor− liquid pipeline flow is also difficult, because there is not enough information to describe the varying characteristics of the interface and to reveal the interaction of heat and mass transfer process between two phases through the varying interface and the unsteady phase change phenomenon.1 Not limited to hydraulic parameters, more recent attention has turned to nonisothermal flow in a pipe or plane channel where some numerical studies are also found.6,17−23 The detailed characteristic of heat transfer is taken into consideration in these mentioned models instead of an average value being represented for the temperature profiles in a pipe.24,25 According to their studies, the wall temperature distribution is different from the assumption of a fully developed isothermal state.26 The energy transfer model has been taken into the flow progress for the optimization of transportation, estimation of corrosion, or prediction of wax deposition.27,28 Concretely, the two-dimensional momentum and three-dimensional energy equations for both phases have been established for dynamic and thermal
1. INTRODUCTION The pipe flow together with phase change is commonly encountered in various heat and mass transfer processes over the past four decades, for instance, in petroleum and chemical processing, steam generating equipment, nuclear reactors, geothermal fields, heat exchangers, cooling systems, and solar energy systems.1−7 In petroleum transportation, two-phase flow characterization is a very common and economic technique, where vapor−liquid two-phase stratified flow is often observed in horizontal or slightly inclined systems.8,9 During the vapor− liquid pipeline transportation process, the liquid holdup and pressure drop in situ are required for determining the pigging period and pigging frequency, and the prediction of the liquid volume is necessary for optimizing the capacity of the liquid receiving facilities. What’s more, a good prediction of temperature contour in both phases and inner wall is essential to master the flow assurance issues, for instance, phase equilibrium, wax deposition, hydrate formation, and pipeline corrosion.10,11 Those predictions are related to the modeling of momentum, energy, and mass transfer for each phase at the wall and the interface, which are significantly different from the single phase flow.12 Therefore, to optimize the design and operation of the above-mentioned systems, the multiphase momentum, heat transfer, and phase equilibrium are in need of being modeled properly.13 It has been observed experimentally that when phase change occurs during the saturated vapor pipeline transportation process, the thermal gradients are created in the wall of the pipeline that lead to severe liquid condensation and stratified © XXXX American Chemical Society
Received: Revised: Accepted: Published: A
May 7, 2017 June 23, 2017 June 30, 2017 June 30, 2017 DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research numerical simulation.8,19,22,29 The smooth or wavy interface between phases was obtained in a different range of flow rates.29 For such a two-phase nonisothermal stratified flow, analytical and numerical heat transfer solutions limited to laminar flow and without phase change have been obtained for fully developed stratified flow under different thermal boundary conditions.19 Then, solutions that are more applicable to fully developed turbulent gas−liquid smooth stratified flow have been obtained through the use of a high Reynolds k − ε model.22 Recently, the steady-state axial momentum and energy equations coupled with a low Reynolds k − ε model were established and solved.28 The pressure drop, liquid height, and temperature field included in the solutions could match well with the experimental data. However, although the EOS (equation of state) was utilized in previous one-dimensional (1D) models to calculate the phase fraction,15,30−37 the flow rate, temperature, and pressure were not coupled with the varying liquid level. From the above numerical simulations, it can be concluded that most of the work is restricted to no phase-change flow regime. Although some studies on hydrodynamics are done for fully developed laminar flow with phase change or done for fully developed turbulent flow with heat transfer, no work is done for the turbulent flow with both heat transfer and phase change, which could absolutely lead to a varying interface. Furthermore, no numerical study has been found for the multicomponent compressible stratified fluid flowing nonisothermally in a longdistance pipeline with phase change. The primary intention of this present study is to develop and solve a model that allows the predictions of corresponding pressure gradient, liquid holdup, mass flow rate, velocity distribution, temperature fields, mole fraction, and physical properties. The model could systematically describe the complex pipe flow and heat/mass transfer process, which includes the two-phase mass conservation equation, three-dimensional steady axial momentum equation, three-dimensional steady heat transfer equation, and phase equilibrium equation. The results are compared to the experimental data and the predictions by other mechanistic models, which appeared in the literatures. In conclusion, this model could have important application for the design of downstream equipment, optimization of transportation rates, determination of pigging frequency, prediction of wax deposition, estimation of corrosion, and guarantee of flow assurance.
Figure 1. Vapor−liquid two phase flow coupled with heat and mass transfer.
as pseudocomponent C7+. (6) The effect of gravity on the P−T flash process is not considered. 2.2. Model Establishment. Several kinds of flow patterns are likely to form for vapor/condensation flowing in pipeline: stratified flow, slug flow, annular flow, and stratified−dispersed flow.39 It is hard to calculate the liquid level, the exact position of liquid film, and the migration patterns of condensed liquid (as shown in Figure 1) via one-dimensional model which could not present the detailed distribution of hydraulic and thermal parameters at pipe cross section. Moreover, the 1D model does not consider the turbulent flow which would lead to different numerical results.4,5,7 Hence, the control volume in three-dimensional coordinates is adopted to discretize the calculating area, as shown in Figure 2. In this section, the governing equations based on physical conservation are chosen and established in three-dimensional coordinates, which include mass conservation equation, momentum conservation equation, energy conservation equation, turbulent flow model and phase change model. 2.2.1. Mass Conservation Equation. The total mass flow rate G is the sum of vapor and liquid mass flow rates and remains constant during the flowing process: GG + G L ≡ G
(1)
The mass transfer rate of phase change within the pipe segment has been taken into consideration due to the vapor phase gradually condensing or the liquid evaporating during the flow process. The change of vapor mass flow rate ΔGG is equal to the opposite number of that of the liquid, ΔGL: ΔGG + ΔG L = ΔG ≡ 0
(2)
The liquid volume fraction HL, vapor volume fraction HG = 1 − HL, pressure, vapor, and liquid velocities as well as temperature are the principal unknowns in the cross-sectional area along the pipeline. The derivation of mass conservation equation is presented here on the basis of the two-fluid approach. The mass conservation equations for each phase within the control volume are given as
2. MATHEMATICAL MODEL OF STRATIFIED VAPOR-LIQUID PIPE FLOW WITH PHASE CHANGE 2.1. Problem Description and Model Assumptions. Because pressure and temperature drop along the pipeline, the saturated vapor of hydrocarbons would condenses gradually, as shown in Figure 1. The condensed liquid would attach to the pipe wall as liquid film or accumulate at the lower part of the pipeline, which could result in continuous change of liquid level.30,31,38 Thus, a three-dimensional model of stratified vapor−liquid two-phase flow coupled with phase change is proposed in this paper. Assumptions are made as follows: (1) Precipitation of condensed liquid is a flash evaporation equilibrium process that occurs in a short moment. (2) Regardless of the attachment to inner wall of the pipe, all the condensed liquid accumulates at the bottom of the pipe. (3) Two-phase flow in vapor−liquid pipeline has a stratified flow pattern as well as a stable developing flow area in every calculated pipeline segment. (4) The smooth vapor−liquid interface model is adopted to describe the interface shape. (5) The heavy components of hydrocarbon are simplified
HLρL
∂ρ ∂wL ∂H + ρL wL L + HLwL L ∂z ∂z ∂P
+ HLwL
HGρG
∂TL
P
∂ρG ∂TG
P
T
∂TL ΔG L = A pipeΔz ∂z
∂ρ ∂wG ∂H + ρG wG G + HGwG G ∂P ∂z ∂z
+ HGwG B
∂ρL
∂P ∂z
∂TG ΔGG = A pipeΔz ∂z
(3)
T
∂P ∂z
(4)
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
Figure 2. Schematic illustration of stratified vapor−liquid flow with phase change.
2.2.2. Momentum Conservation Equation. The law of momentum conservation is a universal principle for any flow system that the varying rate of momentum is equal to the sum of the forces imposed on the control volume. Considering the compressibility of vapor and liquid, the equation of momentum is as follows: ∂(ρG,L wG,L) ∂t =−
(
⎛ ∂uG,L ∂wG,L ⎞ ∂wG,L ( − ρG,L u′G,L w′G,L ) = μt ,G,L ⎜ + ⎟ = μt ,G,L ∂x ⎠ ∂x ⎝ ∂z ⎛ ∂vG,L ∂wG,L ⎞ ∂wG,L ( − ρG,L v′G,L w′G,L ) = μt ,G,L ⎜ + ⎟ = μt ,G,L z y ∂ ∂ ∂y ⎝ ⎠
→) + div(ρG,L wG,Lu⎯⎯⎯G,L
(10)
∂τyz ,G,L ∂τzz ,G,L ∂τxz ,G,L ∂p + + FG,L + + ∂y ∂z ∂z ∂x
τxz ,G,L = μm ,G,L
where −ρG,L u′G,L w′G,L and −ρG,L v′G,L w′G,L represent viscous stresses caused by turbulent flow:
∂uG,L ∂z
+
∂wG,L ∂x
)
a n d τyz ,G,L = μm ,G,L
(
∂vG,L ∂z
With the analysis above, the final momentum equation for vapor and liquid stratified flow is as follows:
(5)
+
∂wG,L ∂y
)
⎛ ∂w ⎞⎞ ⎛ ∂wG,L ⎞⎞ ∂ ⎛ ∂⎛ ⎜⎜Γw,G,L⎜ G,L ⎟⎟⎟ ⎜Γw,G,L⎜ ⎟⎟ + ∂x ⎝ ⎝ ∂x ⎠⎠ ∂y ⎝ ⎝ ∂y ⎠⎠ ∂ρG,L wG,LwG,L dp = − ρG,L g sin θ + dz ∂z
represent viscous stress introduced by molecular viscosity, μm denotes molecular dynamic viscosity, and FG,L = ρG,Lg sin θ. In flow region of development, there is no migration or blending between two-flow lamination: ∂uG,L ∂x
∂vG,L
=0
∂y
=0
where Γw,G,L means effective diffusion coefficient, which is the effective viscosity given as the sum of the molecular and eddy viscosity, Γw,G,L = μm,G,L + μt,G,L. 2.2.3. Turbulent Flow Model. Following a similar approach as Xiao et al.40 and Reboux et al.,41 for two-phase turbulent flow, the eddy viscosity is modeled using the LES (large eddy simulation) turbulence model on the basis of the assumption of nonisotropic turbulence. Meanwhile, changes are made to take account of the progressive attenuation of turbulence close to the wall. The LES model is applicable to the flow with different Reynold number. And the subgrid scale viscosity is written as
(6)
For a steady-state flow in the pipeline, the following formulas can be made: ∂(ρG,L wG,L) ∂t
=0
∂τzz ,G,L ∂z
=0
(7)
Thus, the momentum equation for vapor and liquid can be simplified as ⎛ ∂w ⎞⎞ ⎛ ∂wG,L ⎞⎞ ∂ ⎛ ∂⎛ ⎜⎜μm ,G,L ⎜ G,L ⎟⎟⎟ ⎜μm ,G,L ⎜ ⎟⎟ + ∂x ⎝ ⎝ ∂x ⎠⎠ ∂y ⎝ ⎝ ∂y ⎠⎠ ∂ρG,L wG,LwG,L dp = − ρG,L g sin θ + dz ∂z
μt ,G,L = ρG,L fint,W (CSDSΔ)2 |S|G,L
(12)
The Smagorinsky constant is CS = 0.1−0.2, and the resolved rateof-strain tensor is
(8)
According to the theory of CFD, the N−S equation is applicable to any kind of flow. A direct calculation of the N−S equation requires a high computer capacity, which is not practical in engineering. Hence, an item of Reynolds is introduced:
|S|G,L = (2SijSij)G,L
1/2
⎡⎛ ∂w ⎞2 ⎛ ∂w ⎞2 ⎤1/2 G,L G,L ⎥ = ⎢⎜ ⎟ ⎟ +⎜ ⎢⎣⎝ ∂x ⎠ ⎝ ∂y ⎠ ⎥⎦
(13)
The filter width Δ, is defined as Δ = (Vol)1/3, where Vol is the volume of the computational cell. Damping functions have been introduced into the expressions of turbulent viscosity. Near a wall, a wall-damping function is required to perish the eddy viscosity on the wall:
⎛ ∂w ⎞⎞ ⎛ ∂wG,L ⎞⎞ ∂ ⎛ ∂⎛ ⎜⎜μm ,G,L ⎜ G,L ⎟⎟⎟ ⎜μm ,G,L ⎜ ⎟⎟ + ∂x ⎝ ⎝ ∂x ⎠⎠ ∂y ⎝ ⎝ ∂y ⎠⎠ ∂ ∂ ( −ρG,L u′G,L w′G,L ) + ( −ρG,L v′G,L w′G,L ) ∂x ∂y ∂ρG,L wG,LwG,L dp = − ρG,L g sin θ + dz ∂z
(11)
+
⎛ d+ ⎞ DS = 1 − exp⎜ ⎟ ⎝ 25 ⎠
(9) C
(14) DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research and d+ is the dimensionless distance to pipe wall or vapor−liquid interface and calculated by d+ =
Prt,G,L
(ρτint,W )0.5 d μm
⎡ ⎛μ ⎞ ⎛μ ⎞2 t,G,L t,G,L ⎢ ⎜ ⎟ ⎜ ⎟ = ⎢0.5882 + 0.228⎜ ⎟ − 0.0441⎜ μ ⎟ μ ⎝ m,G,L ⎠ ⎝ m,G,L ⎠ ⎣ −1 ⎛ ⎛μ ⎞⎞⎤ m,G,L ⎥ ⎟⎟ × ⎜⎜1 − exp( −5.165)⎜⎜ ⎟⎟⎥ ⎝ μt,G,L ⎠⎠⎦ ⎝
(15)
Fulgosi et al.42 provides an exponential dependence of f int on the dimensionless distance to the interface or pipe wall:
(19)
2.2.5. Phase Change Model. The pressure and temperature drop along pipelines, which incurs the change of ratio of gas mass volume and liquid mass volume. It is prone to evoke P−T flash evaporation in each component. The mass percentage would change accordingly, influencing the parameters like molar mass, density, viscosity, heat capacity, and thermal conductivity. The vapor enthalpy would change during the process of phase change and dissipate into the vapor or liquid in the form of phase change heat, which results in variation of fluid temperature.33,48,49 The equation of state in P−T Flash involves the PR (Peng− Robinson) equation, which was proposed by Peng−Robinson in 1976.50 It can predict the molar volume more accurately than SRK equation and be applied for polar compounds. Apart from that, the equation is applicable to vapor and liquid at the same time and widely used in the calculation of phase equilibrium.
fint,W = 1 − exp[−1.3 × 10−4d + − 3.6 × 10−4(d +)2 − 1.08 × 10−5(d +)3 ]
(16)
Parameter d means the distance to pipe wall or vapor−liquid interface. The calculation of d is based on correlation:43,44 ⎛ D d = min⎜⎜ − ⎝2
⎞ ⎛ ⎞2 D x 2 + ⎜y − + hL⎟ , |d|⎟⎟ ⎝ ⎠ 2 ⎠
(17)
2.2.4. Heat Transfer Model. It is of great significance to study turbulent flow and the heat transfer mechanism due to their frequent occurrence in many industrial applications, such as heat exchangers, vapor turbines, cooling systems, and nuclear reactors.28 In this study, a quasi-steady state of temperature profile is calculated where axial thermal conduction is neglected. Because the operating temperature is influenced by ambient temperature, fluid properties and hydraulic parameters such as flow rate and liquid level, a heat transfer model for vapor−liquid flow is established. Energy equation means the increase of energy, which is equal to the result of heat flux entering the representative-element volume deducting the work from internal force. Ignoring axial heat conduction and viscous dispersive item, the equation of energy conservation is as follows:
p=
(20)
R2Tc 2 α(Tr) pc
(20a)
a(T ) = acα(Tr) = 0.45724
RTc pc
(20b)
α(Tr) = [1 + m(1 − Tr 0.5)]2
(20c)
m = (0.37464 + 1.54226ω − 0.26992ω 2)
(20d)
b = 0.07780
⎛ ∂T ⎞⎞ ⎛ ∂TG,L ⎞⎞ ∂ ⎛ ∂⎛ ⎜⎜ΓT ,G,L⎜ G,L ⎟⎟⎟ ⎜ΓT ,G,L⎜ ⎟⎟ + ∂x ⎝ ⎝ ∂x ⎠⎠ ∂y ⎝ ⎝ ∂y ⎠⎠ =
a(T ) RT − V−b V (V + b) + b(V − b)
During the P−T flash calculation process, the required parameters are the component of light hydrocarbon, gasification rate, density, molar mass, and enthalpy. The relation between input and output is explained in Figure 3:
∂ρG,L cP ,G,LwG,LTG,L (18)
∂z
ΓT,G,L means effective diffusion coefficient, which is the effective heat transfer coefficient defined as ΓT ,G,L = μm,G,L cP ,G,L Prm,G,L
μm,G,L cP ,G,L
and
Prm,G,L μt,G,L cP ,G,L Prt,G,L
+
μt,G,L cP ,G,L Prt,G,L Figure 3. Input and output of P−T flash calculation.
are respectively molecular thermal
The fluid in the pipe includes Nc sorts of hydrocarbon, i = 1, 2, ..., Nc. The total energy remains constant during the phase change process. When vaporization occurs, the excess heat of the vapor phase is absorbed from the liquid phase. When condensa tion occurs, excess heat is released from the vapor phase. LH represents the latent heat of vaporization or condensation, which is calculated by33
μ
diffusion and eddy diffusion of heat transfer. Prm = cP λm is the molecular Prandtl number. Prt is the turbulent Prandtl number, which is assumed by Jones and Launder.45 as 0.9. However, Hishida et al.46 obtained the conclusion that the value of Prt increases with the distance increase to pipe wall. Only if d+ > 30, Prt is 0.9. When d+ = 10, the value of Prt rises to 1.6 according to the measurement of circular tube of turbulent heat transfer experiments.28 Another calculating of turbulent Prandtl number is given by Kays47 and Mansoori22 on the basis of the direct numerical simulation:
⎧ EG − E L , when condensation occurs LH = ⎨ ⎩ E L − EG , when vaporization occurs ⎪
⎪
D
(21)
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
The relationship between liquid level hL and liquid circulation area AL is as follows:
Therefore, the latent heat may result in the temperature change of the vapor−liquid system. The enthalpy difference in the process of vaporization or condensation can be obtained by virtue of P−T flash calculation. The latent heat helps to retard the temperature change. The process can be calculated as follows:
∫A
A pipeHL = AL = −
ρG,L wG,LL Hlηlξ dη dξ = cP ,GGGΔTG + cP ,LG LΔTL
⎛ 2h ⎞ 1 2 D arccos⎜1 − L ⎟ ⎝ 4 D ⎠
2 ⎛ 2h ⎞ 2h ⎞ 1 2⎛ D ⎜1 − L ⎟ 1 − ⎜1 − L ⎟ ⎝ 4 ⎝ D ⎠ D ⎠
(27)
G,L
The value of HL as well as the calculating method differs in each cross section, as depicted in Figure 5. As shown in Figure 5, the evaporation fraction, density, and molar weight of the vapor and liquid phase can be obtained with the known mole fraction of each component because the temperature and pressure are the same as inlet data in the first cross-section. The liquid hold-up can be calculated as follows:
(22)
The parameters of flow state and related dimensionless parameters are given by GG =
∫A
ρG wGlηlξ dη dξ
GL =
G
∫A
ρL wLlηlξ dη dξ
L
(23)
The properties of fluids are calculated by
HL =
(e , X , Y , ρG , ρL , ρMIX , MG , ML , MMIX , EG , E L , EMIX ) = PR EOS(Z ,P ,T )
μm,G = μm,G (TG ,ρG )
λG = λG(TG)
μm,L = μm,L (PL ,TL ,ML ,ρL ,X ) cP ,L = cP ,L(ρL ,TL)
1+
e ρG ML (1 − e)
(28)
The temperature, pressure, and mole fraction vary in the second and other subsequent cross-sections due to the influence of thermal conduction, pressure drop, and phase change but can be derived from the former section. Hence, it requires the calculation of the evaporation fraction, density, and molar weight in each grid cell. The liquid holdup can be calculated as follows:
(24)
cP ,G = cP ,G(PG ,TG ,MG)
1 ρL MG
(25)
∑j
λL = λL(ρL ,TL)
HL = (26)
∑j
Both the gas and liquid phase are not ideal fluids. The vapor phase is a mixture of multicomponent light hydrocarbons. Thus, the property of the vapor phase is a combination of their quality weighting. The PR EOS equation is acknowledged as currently the most accurate formula to calculate the density of mixed vapor. If the pressure, temperature, and relative density of light hydrocarbon are known, the viscosity can be calculated using experimental formulas. Both the thermal conductivity and the heat capacity at constant pressure for the vapor mixture are related to temperature and pressure, which also can be calculated using the experimental formulas. For the liquid phase, the density is calculated by PR EOS. The other physical properties including viscosity, thermal conductivity and specific heat are respectively a function of temperature, density and other critical properties of components. The hL at cross section of vapor−liquid flow in the pipeline is depicted in Figure 4:
ρMIX, j
ρMIX, j ML, j MMIX, jρL, j
⎡ ML,j
MMIX, j ⎢ ⎣ ρL,j
Aj (1 − ej)
Aj (1 − ej) +
MG, j ρG, j
⎤ Aj ej ⎥ ⎦
(29)
Then, the liquid level hL can be obtained by eq 27. 2.3. Bipolar Cylindrical Coordinate System and Governing Equations. On the basis of the theory of flow and heat transfer, turbulent flow, and phase equilibrium, the model is solved by multiphysical field coupling numerical simulation. The noncircular liquid and vapor domains in stratified pipe flow, as shown schematically in Figures 1 and 2, can be simply modeled with the bipolar coordinate system, which is helpful in solving the problem caused by the inhomogeneity of boundaries.28 2.3.1. Bipolar Cylindrical Coordinate System. The bipolar cylindrical coordinate is composed of two orthogonal circles in rectangular coordinate,8 as shown in Figure 6a. As the flow field in both phases is bounded by a circular pipe wall and a plane interface, the calculation domain has been converted to rectangle form from the anomalous physical domain by adopting the bipolar coordinate system, seen in Figure 6b. Its analytical expression is as follows: (x − a coth η)2 + y 2 = (a csc hη)2 x 2 + (y − a cot ξ)2 = (a csc η)2
(30)
On the basis of eq 30, the conversion can be denoted as shown in a Figure 6a, where cosh η − cos ξ = lη = lξ and lz = 1 are Lame coefficient and defined as scale factors in bipolar coordinate system. According to the geometric relationship, the liquid level is hL and the vapor flows concurrently above it. The angle γ is half of the circumferential angle corresponding to the vapor−liquid interface and calculated by28
⎛ 2h ⎞ γ = arccos⎜1 − L ⎟ ⎝ D ⎠
Figure 4. Illustration of pipe cross-section with light hydrocarbon vapor. E
(31) DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
Figure 5. Illustration of the way to obtain the mole fraction by P−T flash in each grid cell.
Figure 6. Bipolar coordinate system and the conversion of calculation domain.
heat transfer and phase change involve vapor−liquid interface condition, wall boundary condition, symmetrical boundary condition, and inlet boundary condition. As for the vapor− liquid interface, to illustrate the mutual influences among flow, heat transfer, and phase change, the equal interfacial shear stress has been prescribed as the vapor−liquid interface condition, which is related to the fluid properties and velocity distribution and can be calculated as
And the computational domain with bounds is expressed respectively as Vapor: ξ ∈ [γ , π ]
η ∈ ( −∞ , +∞)
Liquid: ξ ∈ [π , π + γ ] η ∈ ( −∞ , +∞)
(32)
In the numerical calculation, the value of ηmax must be precise, Newton and Behnia16 propose ηmax = 6. The velocity in the pipeline is subjected to influences of pressure and liquid holdup. The temperature in turn affects the velocity through affecting on physical properties of the fluid. Momentum equations are set up in this paper to solve the velocity distribution from mutual influences of temperature and phase change. And the mass conservation equation is utilized to solve the pressure gradient and mass flow rates along the pipeline. 2.3.2. Governing Equations in Bipolar Cylindrical Coordinate. With the analysis above, the model in bipolar cylindrical coordinates can be expressed as the equations in Table 1. In fully developed compressible and stratified smooth two-phase pipe flow, the time-averaged steady-state conservation equations for each phase in an arbitrarily shaped domain mapped by an orthogonal η − ξ − z coordinate system can be written as shown in Table 1. In Table 1, dp/dz is the pressure gradient in the axial direction, Pa/m. In the condition of steady flow state, the axial pressure gradient dp/dz in liquid is balanced by the shear stress at the τWL and interface τint, respectively, in vapor by τWG and τint. 2.4. Boundary Condition. The boundary conditions included in vapor−liquid two-phase stratified pipe flow with
τint,G = (μm + μt )G
∂wG ∂ξ
= (μm + μt )L int
∂wL ∂ξ
= τint,L int
(33)
The temperature and heat flux of two phases are respectively equal at vapor liquid interface. qint,G = qint,L
Tint,G = Tint,L
(34)
At the pipe wall, the nonslip condition is applied for velocity in both phases.
wW,G,L = 0
(35)
The temperature boundary condition at the pipe wall of vapor and liquid phases is convective heat transfer, and its coefficient of the pipeline outer wall remains constant. αG,L(TW − TG,L) = λG,L F
dTG,L dn
W
(36)
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research Table 1. Summary of Thermal-Hydraulic Models no.
equation
model presented in this paper HL,GρL,G
1
∂wL,G ∂z
mass conversation
+ ρL,G wL,G
∂ρL,G
+ HL,GwL,G
∂TL,G
∂HL,G ∂z
∂TL,G p
+ HL,GwL,G
∂z
=
∂ρL,G ∂p
∂p ∂z
T
ΔG L,G A pipeΔz
ΔGG + ΔG L = ΔG ≡ 0
2
momentum conversation
⎛ ∂wG,L ⎞⎞ ⎛ ∂wG,L ⎞⎞ 1 ∂ ⎛ 1 ∂ ⎛ ⎜Γw,G,L⎜ ⎜Γw,G,L⎜ ⎟⎟ + ⎟⎟ lηlξ ∂ξ ⎝ ⎝ ∂ξ ⎠⎠ lηlξ ∂η ⎝ ⎝ ∂η ⎠⎠ ∂ρG,L wG,LwG,L dp = − ρG,L g sin θ + ∂z dz
Γw,G,L = μm,G,L + μt,G,L
3
turbulence model
μt,G,L = ρG,L fint,W (CSDSΔ)2 |S|G,L 2 ⎡⎛ ⎛ 1 ∂wG,L ⎞2 ⎤ 1 ∂wG,L ⎞ ⎟⎟ + ⎜⎜ ⎟⎟ ⎥ |S|G,L = (2SijSij)G,L1/2 = ⎢⎜⎜ ⎢⎝ l ∂η ⎠ ξ l ∂ ⎝ ⎠ ⎥⎦ η ξ ⎣
1/2
⎛ D d = min⎜⎜ − ⎝2
⎞ ⎛ ⎞2 D x 2 + ⎜y − + hL⎟ , |d|⎟⎟ ⎝ ⎠ 2 ⎠ ⎛ d+ ⎞ DS = 1 − exp⎜ ⎟ ⎝ 25 ⎠
CS = 0.1−0.2
d+ =
(ρτint,W )0.5 d μm
fint,W = 1 − exp[−1.3 × 10−4d + − 3.6 × 10−4(d +)2 − 1.08 × 10−5(d +)3 ]
4
energy conversation
⎛ ∂TG,L ⎞⎞ ⎛ ∂TG,L ⎞⎞ ∂ρG,L wG,LcP ,G,LTG,L 1 ∂ ⎛ 1 ∂ ⎛ ⎜ΓT ,G,L⎜ ⎜ΓT ,G,L⎜ ⎟⎟ + ⎟⎟ = lηlξ ∂ξ ⎝ ∂z ⎝ ∂ξ ⎠⎠ lηlξ ∂η ⎝ ⎝ ∂η ⎠⎠ μm,G,L cP ,G,L
ΓT ,G,L =
Prm,G,L
+
μt,G,L cP ,G,L
Prm,G,L =
Prt,G,L
cP ,G,L(TG,L) μm,G,L (TG,L) λG,L(TG,L)
⎡ ⎞2 ⎛μ ⎞ ⎛μ t,G,L ⎟ − 0.0441⎜ t,G,L ⎟ Prt,G,L = ⎢⎢0.5882 + 0.228⎜⎜ ⎟ ⎜μ ⎟ ⎝ m,G,L ⎠ ⎝ μm,G,L ⎠ ⎣ −1 ⎛ ⎞⎞⎤ ⎛μ m,G,L ⎟⎥ ⎜ ⎟ ⎜ × ⎜1 − exp( −5.165)⎜ ⎟⎟⎥ ⎝ μt,G,L ⎠⎠⎦ ⎝
5
phase change model
(e , X , Y , ρG , ρL , ρMIX , MG , ML , MMIX , EG , E L , EMIX ) = PR EOS(Z ,P ,T )
∫AG,L ρG,L wG,LLHlηlξ dη dξ = cP ,GGGΔTG + cP ,LGLΔTL ⎧ EG − E L , condensation LH = ⎨ ⎩ E L − EG , vaporization ⎪ ⎪
∑j HL = ∑j
6
auxiliary equations
MMIX, jρL, j
ρMIX, j ⎡ ML, j A (1 MMIX, j ⎢ ⎣ ρL, j j
A pipeHL = GG =
ρMIX, j ML, j
Aj (1 − ej) − ej) +
MG, j ρG, j
⎤ Aj ej ⎥ ⎦
2 ⎛ ⎛ 2h ⎞ 1 ⎛ 2h ⎞ 2h ⎞ 1 2 D arccos⎜1 − L ⎟ − D2⎜1 − L ⎟ 1 − ⎜1 − L ⎟ ⎝ ⎝ 4 D ⎠ 4 ⎝ D ⎠ D ⎠
∫AG ρG wGlηlξ dη dξ
μm,G = μm,G (TG ,ρG )
GL =
λG = λG(TG)
μm,L = μm,L (PL ,TL ,ML ,ρL ,X )
G
∫AL ρL wLlηlξ dη dξ cP ,G = cP ,G(PG ,TG ,MG)
λL = λL(ρL ,TL)
cP ,L = cP ,L(ρL ,TL)
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research The gradients of velocity and temperature are zero at the symmetrical boundary. ∂wG,L ∂η
=0 η=0
∂TG,L ∂η
=0 η=0
(37)
At the inlet of the pipeline, the velocity and temperature field of two phases are respectively equal to the pipe inlet values. wG,L = winlet,G,L
TG,L = Tinlet
(38)
3. NUMERICAL METHOD FOR MODEL SOLVING In bipolar cylindrical coordinates, the general governing equation is as follows: ∂ρwϕ 1 ∂ ⎛ ∂ϕ ⎞ 1 ∂ ⎛ ∂ϕ ⎞ + Sϕ ⎜Γϕ ⎟ = ⎜Γϕ ⎟ + ∂z lηlξ ∂η ⎝ ∂η ⎠ lξlη ∂ξ ⎝ ∂ξ ⎠
(39)
Second-order upwind difference scheme is used for discretization of convective term. As for the diffusion term, it is the central difference scheme. The general discrete format of the governing equation is as follows: aPϕP = aNϕN + aSϕS + aW ϕW + aEϕE + Sϕ
(40)
aP = aN + aS + aW + aE + aB
(41)
Because the model is coupled, the calculation of velocity, temperature, and pressure requires iteration until they converge. In the study, the mass conservation equation is rendered as a judging condition for calculation of the pressure gradient. The corresponding flowchart for programming in every pipe segment Δz is shown in Figure 7. After the convergence of velocity, turbulence, and temperature field, the total mass flow conservation is used to adjust the pressure drop and the rate of gasification calculated via P−T flash is applied to restrict the liquid level. The model solving process consists of four iteration cycles, resulting in four convergence criteria for this simulation. (1) The convergence condition of the iterative solution of the velocity field and the turbulence model is that the infinite norm of the difference between the two iterations of the velocity field is less than 1 × 10−3, that is, ∥wn − wn−1∥∞ ≤ 10−3 or the number of iterations Nwiteration ≥ 300. (2) The convergence condition of the temperature field is ∥Tn − Tn−1∥∞ ≤ 10−3 or the number of iterations NTiteration ≥ 30. (3) The convergence condition of the pressure gradient is that the mass flow is equal to the inlet mass flow, i.e., ∥Gn/Ginlet − 1∥∞ ≤ 10−3 or the number of iterations NGiteration ≥ 30. (4) The convergence condition of the liquid height is ∥hn − hn−1∥∞ ≤ 10−3 or the number of iterations Nhiteration ≥ 30. The condition that the total mass flow remains constant and equal to the inlet mass flow is the key parameter to check convergence, which means when ∥Gn/Ginlet − 1∥∞ ≤ 10−3, the whole model is considered to converge, and the result is credible.
Figure 7. Flowchart of numerical simulation.
temperature field given by Duan et al.8 and experimental data presented by Manabe.20 4.1.1. Validation of Velocity Distribution. The numerical result of the velocity field is compared with Persen’s51 experimental results. The pipeline in the experiment was characterized by an inner diameter of 50 mm, inclination of −0.65°. USG = 3.1 m/s and USL = 0.073 m/s are superficial velocities of experiment 1. USG = 6.4 m/s and USL = 0.073 m/s are superficial velocities of experiment 2. As shown in Figure 8, with a certain superficial velocity, the flow field is obviously asymmetric. The numerical results of velocity distribution are in good agreement with measured experimental data. The vertical centerline velocity shows that the velocity gradient is quite large at vapor−liquid interface or near the pipe wall. 4.1.2. Validation of Temperature Distribution. Manabe20 conducted a loop pipeline experiment to test the characteristics
4. RESULTS AND DISCUSSION The presented model is validated by literature results. Case studies regarding phase change in horizontal pipelines are presented as well. 4.1. Model Validation. Both velocity and temperature are required to calculate the varying liquid level. The experimental data presented by Persen51 is used to validate our hydraulic model in calculating the velocity distribution. Our calculation results of thermal model are compared with the numerical H
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
Therefore, the relative error (RE) should be calculated as follows: relative error =
of heat transfer during vapor−liquid two-phase stratified nonisothermal flow in pipe. The experimental pipeline was with an inner diameter of 52.5 mm and a length of 6.18 m. Natural gas and crude oil were used as experimental fluids. A reverse flow of ethylene glycol in the casing with a flow rate of 0.003680 m3/s and a temperature of 48.88 °C was applied to cool the fluid. The convective heat transfer coefficient of the experimental pipe wall is 810 w/(m2·K). Manabe20 employed thermocouples and RTD (resistance temperature detector) to measure temperature at the inlet, outlet, and the pipe wall of the outlet. Detailed experimental parameters and processes can be referred from Manabe’s work.20 During the process of nonisothermal flow, the average temperature of vapor phase TG and average temperature of liquid phase TL are calculated on the basis of taking the respective velocities as the weighting factor.52 The calculation method is
∫A TG(x ,y) wG(x ,y) dA G
∫A wG(x ,y) dA
(42)
G
TL =
∫A TL(x ,y) wL(x , y) dA L
∫A wL(x ,y) dA
(43)
L
The bulk-temperature drop is defined as
ΔTbulk,measured
(45)
As shown in Table 2, with 0.06 and 0.09 m/s superficial velocities of liquid, the comparison between numerical calculation and experimental results are listed for the outlet under the fully developed two phases stratified flow. The average temperature of the liquid is measured as shown in Table 2. Compared with Duan et al.,8 whose numerical simulations are slightly overpredicted by approximately 8.3% when compared to the measurement, the model in this paper presents improved results with the maximum error of 7.92%, which indicates that the numerical predictions of our model is quite reasonable. 4.2. Example of Nonisothermal Pipe Flow Coupled with Phase Change. Vapor−liquid two-phase flow and heat transfer coupled with phase change have been simulated in this section. The simulated pipeline is with inner diameter of 660 mm and total length of 12000 m. The superficial velocities of vapor and liquid are respectively USG = 4.223 m/s and USL = 0.016 m/s. The pressure gradient is about dp/dz = −17 Pa/m, and the liquid holdup at pipe inlet is HL = 0.032.The mass flow rate of the fluid at pipe inlet is 182.38 kg/s. The fluid of vapor−liquid mixture at the pipe inlet also has a pressure of 10.343 MPa and a temperature of 48 °C. The convective heat transfer coefficient is about 10 w/(m2·K) and the ambient temperature is 15 °C. The simulated pipeline is divided into 12000 segments where the cross-sectional area of each segment has a mesh grid scale of 84 × 122 for numerical calculation. The fluid is a mixture of multicomponent hydrocarbons which are shown in Table 3. The MATLAB 2014a software is used to program and calculate the solution. The computer hardware is the i5-CPU and the memory is 6 GB. Almost 50 h were cost to calculate the 6000 m pipe. The total iteration cycles were nearly 1.2 × 109 times. 4.2.1. Grid Structure. The governing equations are discretized and solved using finite difference schemes on a composite grid with local grid refinements near the interface and near the pipe wall (Figure 9). 4.2.2. Convergence Test. The numerical accuracy is improved as the increasing number of computational grid cells, which consumes more computational resources. Two convergence tests have been performed to determine the convergence rate of the numerical scheme. The average absolute error for variable ϕ with the grid cells number NJ is given by
Figure 8. Comparison of calculated and measured vertical centerline velocity profiles.
TG =
ΔTbulk,calculated − ΔTbulk,measured
NJ
8
ΔTbulk = Tbulk,inlet − Tbulk,outlet
average absolute error =
(44)
∑ j = 1 |ϕcurrent step, j − ϕlast step, j| NJ
(46)
Table 2. Comparison of the Bulk-Temperature Drop of Liquid Phase between Measured by Manabe20 and Predicted by Duan et al.8 as Well as This Modela outlet liquid phase bulk-temperature drop test no. USL (m/s) 1 2 3 4 5 6 a
0.06 0.06 0.06 0.09 0.09 0.09
20
USG (m/s)
20
0.17 0.34 2.33 0.14 0.79 2.19
inlet temp (°C) 66.86 63.09 64.79 66.60 63.53 59.96
20
outlet temp measured (°C) 63.38 60.44 61.93 63.46 60.82 56.92
20
measured (°C) 3.48 2.65 2.86 3.14 2.71 3.04
20
predicted (°C)8 RE (%)8 calculated (°C) RE (%) 3.66 2.87 2.92 3.15 2.76 3.25
5.17 8.30 2.10 0.32 1.85 6.91
3.39 2.86 3.08 2.92 2.75 2.88
−2.59 7.92 7.69 −7.01 1.48 −5.26
Data taken from refs 8 and 20. I
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research Table 3. Chemical Composition of the Light Hydrocarbons Used in Current Study mole percent (%)
CH4
C2H6
C3H8
i-C4H10
n-C4H10
i-C5H12
n-C5H12
n-C6H14
C7+
75.84
6.73
5.98
4.05
3.54
1.85
0.74
0.69
0.59
Figure 9. Composite grid with local refinement and with adaptive grid spacing.
Figure 10. Numerical convergence test for hydrodynamic and thermal calculation.
liquid mass flow rate leads to further rise of liquid level, as shown in Figure 11b. The liquid holdup first increases until it reaches a maximum value and then gradually decreases. The reason behind this is as follows: The increase of liquid holdup results from the liquid precipitation caused by dominant temperature drop. Due to the large difference between fluid temperature and ambient temperature, the amount of liquid precipitation is greater than liquid evaporation. On the contrary, the decrease of liquid holdup is lead by liquid evaporation due to dominant pressure drop. Being same to liquid holdup, the liquid mass flow rate maintains the same trend, that is, gradually increasing to reach a maximum value and then gradually reduced, which is depicted in Figure 11c. As the total mass flow is constant, the mass flow rate of the vapor phase decreases first and then increases. When compared with the process of evaporation, the precipitation process caused by temperature drop is transient and intense, which is related to the temperature difference between the inside and outside of the pipeline and to the convective heat transfer coefficient. The tendency of temperature drop is similar to that in the existing research.33 But there also exists a difference between vapor bulk temperature, liquid bulk temperature, and the total
The average absolute errors of velocity and temperature in both phases are chosen for hydrodynamic and thermal convergence test, respectively. It can be seen in Figure 10, the absolute error in velocity and temperature drops with the increase of iteration steps and the number of computational grid cells. When the 84 × 122 grid structure is used, the calculation has achieved a good converging result, so the subsequent calculation uses the grid with current structure. 4.2.3. Pressure Gradient, Liquid Level, Fluid Mass Flow Rate, and Temperature along the Pipeline. The pressure gradient, liquid level, fluid mass flow rate, and temperature along the pipeline obtained in the condition of condensation production have been compared with that found in literature.33 When the phase change behavior is considered along pipe flow, the vapor in the gas phase starts to condense to liquid, which begins from the pipe inlet due to the significant temperature drop at pipe wall. The pressure drop fit well with the simulated results presented by Sadegh,33 as shown in Figure 11a. During the condensing process, vapor mass flow rate gradually reduces, as shown in Figure 11c, and the liquid mass flow rate increases due to the constant total mass flow rate. The increase of J
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
Figure 11. Pressure gradient, liquid level, fluid mass flow rate, and temperature along the pipeline.
Figure 12. Velocity distribution at a pipe length of 6000 m. K
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
Figure 13. Temperature distribution at a pipe length of 6000 m.
bulk temperature, which cannot be revealed by a one-dimensional model. The liquid bulk temperature is always lower than the vapor bulk temperature whereas the vapor bulk temperature is almost always equal to the total bulk temperature. Latent heat is revealed during the vapor condensing process, which slows down the temperature drop, as shown in Figure 11d. 4.2.4. Velocity Distribution at a Pipe Length of 6000 m. Through solving the model, with phase change happening, it can be found that the pressure gradient is 17.203 Pa/m and the liquid level is 118.72 mm when the axial distance reaches 6000 m. Figure 10 shows the distribution of vapor and liquid axial velocity w in the condition of fully developed vapor−liquid twophase flow. From Figure 12a,b, it can be seen that the velocity of the vapor phase slows down while approaching either the pipe wall or the vapor−liquid interface because of the hindering effects and fluid viscosity. The velocity of the liquid phase keeps increasing from the pipe wall to the interface. The maximum velocity at the pipe cross section occurs within the vapor phase.
Figure 12c reveals the velocity distribution at the centerline of the pipe. By the dragging force of the interface, the velocity of the liquid phase reaches the maximum value at the interface while the velocity of the vapor phase reaches the maximum value at the location between the interface and its bulk center. The liquid phase slows down while approaching the pipe wall because of the hindering of the pipe wall and its high viscosity. Figure 12d represents the velocity distribution at the interface, where the velocity at the interface decreases from the center to the wall of the pipe. In the pipe center area, the velocity decreases so slowly as to be nearly constant, whereas near the pipe wall area, velocity decreases dramatically. 4.2.5. Temperature Distribution at a Pipe Length of 6000 m. From Figure 13a,b, it can be seen that the temperatures of both the vapor and liquid phases drop while approaching the pipe wall because of the lowest ambient temperature and the convective heat transfer effects. The temperature of the liquid phase keeps increasing from the pipe wall to the interface. The maximum L
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research Table 4. Mole Fraction of Each Component at Pipe Inlet in Both Phases mole fraction in vapor phase (%) mole fraction in liquid phase (%)
CH4
C2H6
C3H8
i-C4H10
n-C4H10
i-C5H12
n-C5H12
n-C6H14
C7 +
total
77.99 43.06
6.67 7.57
5.66 10.81
3.65 10.07
3.12 9.97
1.51 7.02
0.59 3.07
0.47 3.97
0.34 4.46
100 100
temperature at the pipe cross section exists within the vapor phase near the interface. In Figure 13c, the temperature of the vapor phase at the pipe wall is lower than that of the liquid phase when convective heat transfer exists due to the smaller heat carried by the vapor phase than the liquid phase. Thus, a lower specific heat capacity results in a bigger temperature drop of the vapor phase at the pipe wall. The thermal conductivity of the liquid phase is greater than that of vapor phase; hence the temperature gradient in the liquid phase is smaller than that in the vapor phase, and the bulk average temperature of the liquid phase is lower than that of the vapor phase. Heat is transferring from the vapor phase to the liquid phase through the interface, which makes the temperature of the liquid phase rise and reduces the temperature difference between the two phases. In Figure 13d, the temperature is almost constant in the center area of pipe cross section and decreases toward the pipe wall, where temperature drops dramatically in approaching the pipe wall. The temperature distribution at the pipe wall in both the vapor and liquid phases is shown in Figure 13f. The points near the wall and at the interface are selected to illustrate the temperature distribution, as shown in Figure 13e. From the top to the interface, temperature at pipe wall of vapor phase increases first and then decreases, where the lowest point take place at the pipe top and the highest point near the interface. From the interface to the bottom, the temperature at pipe wall of the liquid phase increases, where the lowest point appears at the pipe wall near the interface and the highest point at the bottom of the pipe. 4.2.6. Mole Fraction, Density Distribution, and Liquid Level along the Pipeline. The mole fraction of each component, density distribution, temperature distribution, and liquid level along the pipeline are obtained in the condition of condensation production. Mole fractions of each component at pipe inlet in both phases are shown in Table 4. In the vapor phase, the mole fraction of methane is larger than all the other light hydrocarbons (C2+). In the liquid phase, the mole fractions of the other light hydrocarbons are bigger than that in the vapor phase, but the methane is still the main component. Mole fractions of each component in vapor phase are shown in Figure 14. The content of methane in the vapor phase increases when flowing in the pipe while the content of the other light hydrocarbons become less and less in vapor phase. During the condensing process which is dominated by temperature drop, the methane keeps evaporating, whereas during the evaporating process, which is dominated by pressure drop, the other light hydrocarbons keep condensing. Meanwhile, the bigger the molar mass is, the faster the condensing rate is. Six pipe cross sections located at every 1000 m along the pipeline are selected to illustrate the change of density and temperature distribution, as shown in Figures 15 and 16. It is exactly because the pressure on the cross-section is the same, so the uneven distribution of the temperature leads to the uneven distribution of the fluid density. In the two phases, the density distribution is opposite to the temperature distribution. The high temperature means low density. The temperature value
Figure 14. Mole fraction of each component in vapor phase.
Figure 15. Density distribution depicted in 6 pipe cross sections along the pipeline (every 1000 m).
Figure 16. Temperature depicted in 6 pipe cross sections along the pipeline (every 1000 m). M
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
Figure 17. Liquid level depicted in 12 pipe cross sections along the pipeline (every 1000 m).
hydrocarbon transportation process in the pipeline, the temperature drop leads to the reduction of vapor mass flow rate and the rise of liquid level as well as mass flow rate. A larger temperature drop results in a bigger liquid holdup whereas a larger pressure drop causes a smaller liquid holdup due to the change of physical properties and phase equilibrium. After the increase of liquid holdup caused by dominant temperature drop reaching the maximum value, the decrease of liquid holdup maintains its trend until the pipe outlet, which results from liquid evaporation due to dominant pressure drop. The highest values of both velocity and temperature locate in the vapor phase. The liquid bulk temperature is always lower than the vapor bulk temperature whereas the vapor bulk temperature is almost always equal to the total bulk temperature, which cannot be revealed by the one-dimensional model. Latent heat is revealed during the vapor condensing process, which slows down the temperature drop. The average temperature and velocity of the liquid are lower than that of the vapor. But near the pipe wall the temperature of the liquid is higher than that of the vapor. When the fluid flows in the pipeline, the content of methane in the vapor phase increases all the time while the content of the other light hydrocarbons (C2+) become less and less. The condensing rates of the other light hydrocarbons have a positive correlation with their molar mass. The temperature value distributed at every single cross section can be ranked in descending order: the interior of the vapor phase, most parts in the liquid phase, and the vapor phase near the wall. However, the rank of density is the liquid phase, the vapor phase near the wall, and the interior of the vapor phase. Along the pipeline, the temperature in the area referenced above decreases gradually while the density increases gradually. Therefore, the temperature of the sequential cross sections tends to be the same and the density distribution within the two phases gradually becomes uniform. As for the varying trend of the liquid level presented in three-dimensional coordinate system, it has the same trend of the liquid holdup which first increases until it reaches a maximum value and then gradually decreases. Thus, this model can be utilized to accurately predict pressure gradient, velocity, temperature field, liquid holdup, fluid physical properties, and mole fraction, which are essential to the determination of pipe size, design of downstream equipment and guarantee of flow assurance.
distributed at every single cross section can be ranked in descending order: the interior of the vapor phase, most parts in the liquid phase, and the vapor phase near the wall. However, the rank of density is the liquid phase, the vapor phase near the wall, and the interior of the vapor phase. Along the pipeline, the temperature in the area referenced above decreases gradually while the density increases gradually. Therefore, the temperature of the sequential cross sections tends to be the same and the density distribution within the two phases gradually becomes uniform. Figure 17 illustrates the selected 12 pipe cross sections located at every 1000 m along the pipeline, where the varying trend of the liquid level is presented in a three-dimensional coordinate system. The minimum liquid level is 103.65 mm at the pipeline inlet and the maximum one is 121.5 mm, which exists between pipeline lengths of 3120 and 3740 m. The liquid level at the outlet is about 113 mm and keeps the declining trend, which can also be found in Figure 11b.
5. CONCLUSION The vapor−liquid two-phase pipe flow and heat transfer are studied by virtue of numerical simulation in light hydrocarbon transportation pipeline coupled with hydraulics, thermodynamics, and phase change. A three-dimensional nonisothermal vapor−liquid stratified flow model including the phase change model in a bipolar coordinate system has been established, where the LES turbulence model is utilized to simulate the turbulence flow and the wall attenuation function is used to describe the inadequacy performance of the vapor−liquid interface. The vapor phase and the liquid phase are both considered to be compressible and the PR equation of state is chosen for the vapor−liquid equilibrium calculation where the multicomponent hydrocarbon flash calculation is used to evaluate the physical properties, gasification rate, and enthalpy departure of the phases. The P−T flash calculation has been applied to predict the varying liquid level and the multicomponent mass fraction in each phase during the process of vapor/liquid stratified pipe flow. After the convergence of velocity, turbulence, and temperature field, the total mass flow conservation is used to adjust the pressure drop. The predictions of velocity distribution and temperature value have been compared with experimental data and showed acceptable agreement with the experimental results given by Persen51 and Manabe,20 respectively. The axial pressure gradient, liquid holdup, velocity, and temperature fields have been presented. The fluid mass flow rate, mole fraction, density distribution, and liquid level along the pipeline are also given out. The simulation results indicate that the influence of pressure and temperature on liquid holdup is different. During the light
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Tel: +86 15101186627. N
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research ORCID
X = mole fraction of each component in liquid phase, dimensionless Y = mole fraction of each component in vapor phase, dimensionless Z = mole fraction of each component, dimensionless
Guoxi He: 0000-0001-7982-4448 Notes
The author declares no competing financial interest.
■
Greek Symbols
ACKNOWLEDGMENTS
α = convective heat transfer coefficient, W/(m2·K) α(Tr) = parameter in PR equation of state, dimensionless γ = half of the circumferential angle corresponding to the interface, rad Γw = effective diffusion coefficient as viscosity, Pa·s ΓT = effective diffusion coefficient as thermal conductivity, W/ (m2·K) δ = circumferential distance at pipe wall, rad Δ = filter width, m ΔG = relative density of light hydrocarbon vapor, dimensionless η, ξ = bipolar cylindrical coordinate system, rad ηmax = maximum value of η, rad θ = pipe angle from the horizontal, rad λ = conductivity coefficient, W/(m·K) μ = dynamic viscosity, Pa·s μm, μt = molecular, turbulent dynamic viscosity, Pa·s μi = viscosity of component i, Pa·s ρ = density of the fluid, kg/m3 τ = viscous stress, Pa τint,W = shear stress at interface or pipe wall, Pa ω = acentric factor, dimensionless
This work was supported by the National Natural Science Foundation of China (Grant 51474228); and the Beijing Scientific Research and Graduate Joint Training Program (Grant ZX20150440).
■
NOTATIONS a = distance from centerline to pipe wall at interface, m ac = coefficient in PR equation of state, dimensionless a(T) = parameter in PR equation of state, dimensionless A = area, m2 Apipe = cross sectional area, m2 b = parameter in PR equation of state, dimensionless cP = specific heat capacity at constant pressure, kJ/(kg·K) CS = Smagorinsky constant, dimensionless d = distance to interface or pipe wall, m d+ = dimensionless distance to interface or pipe wall, dimensionless D = pipeline diameter, m DS = wall-damping function on d+, dimensionless e = rate of gasification, kmol/kmol E = enthalpy of the fluid, kJ/kmol f int,W = factor of dimensionless distance to interface or pipe wall, dimensionless F = body force, Pa/m g = gravity acceleration, m/s2 G = mass flow rate of the fluid, kg/s ΔG = change of mass flow rate of the fluid, kg/s hL = liquid level, m HL, HG = liquid holdup, vapor holdup, dimensionless K = phase equilibrium constant, dimensionless lη, lξ, lz = Lame coefficient, dimensionless LH = latent heat energy, kJ/kg m = coefficient in PR equation of state, dimensionless M = molar mass of the fluid, kJ/kmol MFC = mole fraction of each component, dimensionless Nc = total number of light hydrocarbon, dimensionless NJ = total number of grid cells, dimensionless p = system pressure, Pa pc = critical pressure, Pa Pr = Prandtl number, dimensionless Prt = turbulent Prandtl number, dimensionless q = heat flux, W/m2 R = general vapor constant, 8.314 kJ/(kmol·K) |S| = module of rate-of-strain tensor, s−1 t = time, s T = temperature of the fluid, K Tc = critical temperature, K Tr = relative temperature, dimensionless T∞ = ambient temperature, K ΔT = change of temperature of the fluid, K u, v = the time-averaged x, y axial velocity m/s u⃗ = velocity vector, m/s V = mole volume, m3/kmol w = the time-averaged z axial velocity of the fluid m/s x, y, z = rectangular coordinate system, m
Superscripts
′ = instantaneous value + = dimensionless value − = average value Subscripts
■
bulk = value of bulk c = critical C7+ = pseudocomponent G = vapor phase i, [i] = component number int = value at interface inlet = value at pipe inlet j = grid cell number k = pipe cross section number L = liquid phase m = molecular MIX = mixture of vapor and liquid phase outlet = value at pipe inlet r = relative t = turbulence W = value at wall x, y, z = axis direction
REFERENCES
(1) Akansu, S. O. Heat transfers and pressure drops for porous-ring turbulators in a circular pipe. Appl. Energy 2006, 83, 280−298. (2) Siavashi, M.; Bahrami, H. R. T.; Saffari, H. Numerical investigation of flow characteristics, heat transfer and entropy generation of nanofluid flow inside an annular pipe partially or completely filled with porous media using two-phase mixture model. Energy 2015, 93, 2451−2466. (3) Wróblewski, W.; Dykas, S. Two-fluid model with droplet size distribution for condensing steam flows. Energy 2016, 106, 112−120. (4) Wang, H.; Zhu, T. A novel model for steam transportation considering drainage loss in pipeline networks. Appl. Energy 2017, 188, 178−189.
O
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research (5) Revellin, R.; Lips, S.; Khandekar, S. Local entropy generation for saturated two-phase flow. Energy 2009, 34 (9), 1113−1121. (6) Zhang, L.; Yang, S.; Xu, H. Experimental study on condensation heat transfer characteristics of steam on horizontal twisted elliptical tubes. Appl. Energy 2012, 97, 881−887. (7) Ferreira, R. B.; Falcão, D. S.; Oliveira, V. B. Numerical simulations of two-phase flow in an anode gas channel of a proton exchange membrane fuel cell. Energy 2015, 82, 619−628. (8) Duan, J.; Liu, H.; Gong, J.; Jiao, G. Heat transfer for fully developed stratified wavy gas−liquid two-phase flow in a circular cross-section receiver. Sol. Energy 2015, 118, 338−349. (9) Lun, I.; Calay, R. K.; Holdo, A. E. Modelling two-phase flows using CFD. Appl. Energy 1996, 53, 299−314. (10) Rakopoulos, C. D.; Kosmadakis, G. M.; Pariotis, E. G. Critical evaluation of current heat transfer models used in CFD in-cylinder engine simulations and establishment of a comprehensive wall-function formulation. Appl. Energy 2010, 87, 1612−1630. (11) Guan, X.; Zhang, D.; Wang, J.; Jin, Y.; Li, Y. Numerical and electrochemical analyses on carbon dioxide corrosion of X80 pipeline steel under different liquid film thicknesses in NACE solution. J. Nat. Gas Sci. Eng. 2017, 37, 199−216. (12) Hassam, N. C.; Ben, R. H. Analysis of the thermal cooling capacity of heat pipes under a low Reynolds number flow. Appl. Therm. Eng. 2014, 71, 559−572. (13) Chen, H.; Xu, J.; Li, Z.; Xing, F.; Xie, J. Stratified two-phase flow pattern modulation in a horizontal tube by the mesh pore cylinder surface. Appl. Energy 2013, 112, 1283−1290. (14) Almanza, R.; Lentz, A.; Jimeenez, G. Receiver behavior in direct steam generation with parabolic troughs. Sol. Energy 1997, 61, 275−278. (15) Dukhovnaya, Y.; Adewumi, M. A. Simulation of non-isothermal transients in gas/condensate pipelines using TVD scheme. Powder Technol. 2000, 112 (1), 163−171. (16) Newton, C. H.; Behnia, M. A numerical model of stratified wavy gas−liquid pipe flow. Chem. Eng. Sci. 2001, 56, 6851−6861. (17) Berthelsen, P. A.; Ytrehus, T. Numerical modeling of stratified turbulent two- and three-phase pipe flow with arbitrary shaped interfaces. The 5th International Conference on Multiphase Flow, Yokohama, Japan, May 30−June 4, 2004; North-Holland: Amsterdam, 2004. (18) Newton, C. H.; Behnia, M. Numerical calculation of turbulent stratified gas−liquid pipe flows. Int. J. Multiphase Flow 2000, 26, 327− 337. (19) Gada, V. H.; Datta, D.; Sharm, A. Analytical and numerical study for two-phase stratified-flow in a plane channel subjected to different thermal boundary conditions. Int. J. Therm. Sci. 2013, 71, 88−102. (20) Manabe, R. A comprehensive mechanic heat transfer model for two-phase flow with high pressure flow pattern validation. Ph.D. thesis, Department of Petroleum Engineering, University of Tulsa, 2001. (21) Fontoura, V. R.; Matos, E. M.; Nunhez, J. R. A three-dimensional two-phase flow model with phase change inside a tube of petrochemical pre-heaters. Fuel 2013, 110, 196−203. (22) Mansoori, Z.; Yoosefabadi, Z. T.; Saffar-Avval, M. Two Dimensional Hydro Dynamic and Thermal Modeling of a Turbulent Two Phase Stratified Gas-Liquid Pipe Flow//ASME 2009 Fluids Engineering Division Summer Meeting. American Society of Mechanical Engineers 2009, 753−758. (23) Singh, A. K.; Goerke, U. J.; Kolditz, O. Numerical simulation of non-isothermal compositional gas flow: application to carbon dioxide injection into gas reservoirs. Energy 2011, 36, 3446−3458. (24) Babus'Haq, B; Probert, D. Single and multiphase convective heat transfer. Appl. Energy 1992, 43, 291−292. (25) Gong, G.; Chen, F.; Su, H.; Zhou, J. Thermodynamic simulation of condensation heat recovery characteristics of a single stage centrifugal chiller in a hotel. Appl. Energy 2012, 91, 326−333. (26) Hu, H.; Zhang, C. A modified k−e turbulence model for the simulation of twophase flow and heat transfer in condensers. Int. J. Heat Mass Transfer 2007, 50, 1641−1648.
(27) Sarica, C.; Panacharoensawad, E. Review of paraffin deposition research under multiphase flow conditions. Energy Fuels 2012, 26, 3968−3978. (28) Duan, J.; Gong, J.; Yao, H. Numerical modeling for stratified gas− liquid flow and heat transfer in pipeline. Appl. Energy 2014, 115, 83−94. (29) Ullmann, A.; Brauner, N. Closure relations for two-fluid models for two-phase stratified smooth and stratified wavy flows. Int. J. Multiphase Flow 2006, 32 (1), 82−105. (30) Schouten, J. A.; Janssen-van, R. R.; Michels, J. P. Condensation in gas transmission pipelines. Int. J. Hydrogen Energy 2005, 30 (6), 661− 668. (31) Zhou, J.; Adewumi, M. A. Transients in gas-condensate natural gas pipelines. Transactions-American Society Of Mechanical Engineers. J. Energy Resour. Technol. 1998, 120, 32−40. (32) Zaghloul, J. S. Multiphase Analysis of Three-phase (gascondensate-water) Flow in Pipes. Ph.D. Thesis, The Pennsylvania State University, 2006. (33) Sadegh, A. A.; Adewumi, M. A. Temperature distribution in natural gas/condensate pipelines using a hydrodynamic model. SPE Eastern Regional Meeting; Society of Petroleum Engineers, 2005. (34) Jin, T.; Ityokumbul, M. T. Network modelling and prediction of retrograde gas behavior in natural gas pipeline systems. Int. J. Eng. Syst. Model. Simul. 2016, 8 (3), 169−182. (35) Mucharam, L.; Adewumi, M. A.; Watson, R. W. Study of gas condensation in transmission pipelines with a hydrodynamic model. SPE Prod. Eng. 1990, 5 (03), 236−242. (36) Vincent, P. A.; Adewumi, M. A. Engineering design of gascondensate pipelines with a compositional hydrodynamic model. SPE Prod. Eng. 1990, 5 (04), 381−386. (37) Abbaspour, M.; Chapman, K. S.; Glasgow, L. A. Transient modeling of non-isothermal, dispersed two-phase flow in natural gas pipelines. Appl. Math. Model. 2010, 34 (2), 495−507. (38) Oliemans, R. V. Modelling of Gas-Condensate Flow in Horizontal and Inclined Pipes. Proceedings of the ASME Pipeline Engineering Symposium-ETCE, Dallas, 1987; p 73. (39) Pitton, E.; Ciandri, P.; Margarone, M. An experimental study of stratified−dispersed flow in horizontal pipes. Int. J. Multiphase Flow 2014, 67, 92−103. (40) Xiao, F.; Dianat, M.; McGuirk, J. J. LES of turbulent liquid jet primary breakup in turbulent coaxial air flow. Int. J. Multiphase Flow 2014, 60, 103−18. (41) Reboux, S.; Sagaut, P.; Lakehal, D. Large-eddy simulation of sheared interfacial flow. Phys. Fluids 2006, 18 (10), 105105. (42) Fulgosi, M.; Lakehal, D.; Banerjee, S. Direct numerical simulation of turbulence in a sheared air−water flow with a deformable interface. J. Fluid Mech. 2003, 482, 319−345. (43) Berthelsen, P. A.; Ytrehus, T. Calculations of stratified wavy twophase flow in pipes. Int. J. Multiphase Flow 2005, 31 (5), 571−592. (44) Berthelsen, P. A.; Ytrehus, T. Stratified smooth two-phase flow using the immersed interface method. Comput. Fluids 2007, 36 (7), 1273−1289. (45) Jones, W. P.; Lauder, B. E. The calculation of low-Reynoldsnumber phenomena with a two-equation model of turbulence. Int. J. Heat Mass Transfer 1973, 16 (6), 1119−1130. (46) Hishida, M.; Nagano, Y.; Tagawa, M. Transport processes of heat and momentum in the wall region of turbulent pipe flow. Proceedings of the 8th International Heat Transfer Conference; Hemisphere Publishing Corp.: Washington, DC, 1986; Vol. 3, pp 925−930. (47) Kays, W. M. Turbulent Prandtl numberwhere are we? J. Heat Transfer 1994, 116 (2), 284−295. (48) Deng, D.; Gong, J. Prediction of Transient Behaviors of GasCondensate Two-Phase Flow in Pipelines with Low Liquid Loading. 2006 International Pipeline Conference; American Society of Mechanical Engineers, 2006; p 175−182. (49) Helgans, B.; Richter, D. H. Turbulent latent and sensible heat flux in the presence of evaporative droplets. Int. J. Multiphase Flow 2016, 78, 1−11. P
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research (50) Roboinson, D. B.; Peng, D. Y. The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs; Gas Processors Association: Tulsa, OK, 1978. (51) Persen, L. N. An experimental study of slug flow in a slightly inclined horizontal pipeline and of disturbance waves appearing in annular flow in vertical risers. Proceedings of the 2nd International Conference on Multiphase Flow; North-Holland: Amsterdam, 1985. (52) Deen, W. M. Analysis of Transport Phenomena. Topics in Chemical Engineering; Oxford University Press: New York, 1998.
Q
DOI: 10.1021/acs.iecr.7b01898 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
